求数值近似开平方算法.doc

求数值近似开平方算法
200721101008 陈林

摘要：本文分别用牛顿迭代法、蒙特卡洛实验法和样条插值方法求数值的近似开方值，并用数值3开平方为算例，简单比较了三种方法的优劣。

一． 牛顿迭代法

1.原理分析：牛顿迭代法实际上是由不动点迭代法的原理结合泰勒展开式构造迭代公式，所以迭代初值必须满足迭代收敛条件，设，则，由初值条件解得，不妨设，，则牛顿迭代公式为：
。

2.实验程序及注释：
x0=input('input x0='); %设定迭代初值，本例为1
er=1;
n=0; %定义一个计数器记录循环次数
while er>0.001
x0.5*(x0+3/x0); %牛顿迭代公式
er=abs(x-x0);
x0=x; %数值解
n=n+1;
end
x
sqrt3=sqrt(3) %真实数值
n
er
clear %清除内存变量

3.实验结果及分析：
x0
x
sqrt3
n
er
1
1.7321
1.7321
4
9.2047e-005

因为计算机的截断误差和matlab数值显示格式的原因，显示出来的数值解和真实数值是一样的，但从误差变量er的值我们可以看到其计算的误差。误差在小数点后的第5位，保证了5位的有效数字，即1.7321的迭代结果是准确的。迭代4次就达到了要求的精度，说明了牛顿迭代法的收敛速度很快，可以证明它的平方收敛的。

4.注释及花絮（参考文献 Chris Lomont, fast inverse square root, 2000, Mathematics Subject Classification）
牛顿迭代法初值的选取影响收敛的速度，而且有初值条件，不满足初值条件，迭代公式可能发散，这也是牛顿迭代法的缺陷之一；其次迭代公式中要求解方程的一阶导数，而很多函数无法求导或是求其解析导函数很困难，无法用初等函数表示。所以虽然牛顿迭代法能达到平方收敛，但应用却不是想象那么好。
关于迭代初值的研究，很多数学家对此很感兴趣。如卡马克就用0x5f3759df作为初值（前两位0x表示16进制，后面才是数据，换算成十进制为1597463007），达到了很好的收敛效果。但我用这个数作为迭代初值却没有看出可以减少迭代次数，感觉还是取解附近的值作为迭代初值好些，也可能是我所选的数值太小了吧。而对于这个神奇的初值，大家都不知道卡马克当时是怎么算出来的，因为效果很好又具有神秘色彩，被大家誉为Magic Number。
普渡大学的数学家Chris Lomont看了以后觉得有趣，决定要研究一下卡马克弄出来的这个猜测值有什么奥秘。Lomont也是个牛人，在精心研究之后从理论上也推导出一个最佳猜测值，和卡马克的数字非常接近, 0x5f37642f。卡马克真牛，他是外星人吗？
传奇并没有在这里结束。Lomont计算出结果以后非常满意，于是拿自己计算出的起始值和卡马克的神秘数字做比赛，看看谁的数字能够更快更精确的求得平方根。结果是卡马克赢了... 谁也不知道卡马克是怎么找到这个数字的。
最后Lomont怒了，采用暴力方法一个数字一个数字试过来，终于找到一个比卡马克数字要好上那么一丁点的数字，虽然实际上这两个数字所产生的结果非常近似，这个暴力得出的数字是0x5f375a86。
Lomont为此写下一篇论文，"Fast Inverse Square Root"。（有链接）

二． 蒙特卡洛实验法
蒙特卡洛(Monte Carlo)方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战进研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙上了一层神秘色彩。
　　Monte Carlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针试验的方法来决定圆周率π。本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。
　　考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？Monte Carlo方法是这样一种“随机化”的方法：向该正方形“随机地”投掷N个点落于“图形”内，则该“图形”的面积近似为M/N。
　　可用民意测验来作一个不严格的比喻。民意测验的人不是征询每一个登记选民的意见，而是通过对选民进行小规模的抽样调查来确定可能的优胜者。其基本思想是一样的。
　　科技计算中的问题比这要复杂得多。比如金融衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”Course Dimensionality)，传统的数值方法难以对付（即使使用速度最快的计算机）。Monte Carlo方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率，科学家们提出了许多所谓的“方差缩减”技巧。
另一类形式与Monte Carlo方法相似，但理论基础不同的方法—“拟蒙特卡罗方法”(Quasi-Monte Carlo方法)—近年来也获得迅速发展。我国数学家华罗庚、王元提出的“华—王”方法即是其中的一例。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为Low Discrepancy Sequences)代替Monte Carlo方法中的随机数序列。对某些问题该方法的实际速度一般可比Monte Carlo方法提出高数百倍，并可计算精确度。

1. 原理分析：实际原理就是随机数的统计特性，得到两个事件集合的发生的概率之比等于集合对应图形区域面积之比又约等于集合对应的随机点数之比。在第一象限内构造C*C的矩形域和由函数在[0，C]与轴围成的区域，并在矩形域中产生随机点，设在矩形域中的点数为（也为总的随机点数），在函数在[0，C]与轴围成的区域中的点数为。显然，矩形面积为9，函数在[0，C]与轴围成的区域面积为积分,所以，变形得，还是以数值3为算例，则，代入上式得。

2 实验程序及注释：
n=input('input n='); %总随机点数
x=3*rand(n,1);y=3*rand(n,1); %在3*3正方形中产生n个随机点
m=sum(x.^0.5-y>0); %求落在y=x^0.5，x=3与y=0所围成区域点数
sqrt=9*m/(2*n) %利用motecarlo思想所推导的等式
sqrt3=sqrt(3) %真实值
er=abs(sqrt-sqrt3)/sqrt3 %求相对误差
clear %清除内存变量

3. 实验结果及分析：
n
sqrt
sqrt3
er
2000
1.6718
1.7321
0.0348
5000
1.7496
1.7321
0.0101
10000
1.7037
1.7321
0.0164
由实验结果可以看出蒙特卡洛实验方法并不理想，增加随机点数的方法并不能提高精度，甚至可能降低精度。尽管相对误差很小，但有效数字始终只有1至2位。对于无法获取准确数据的模型，蒙特卡洛方法还是提供了一种有效的方法。

三． 样条插值法
样条插值是在分段三次埃尔米特插值法改进基础上得到的，具有很好的光滑性。
1. 原理分析：以数值3开平方为算例，在3附近选取容易计算开平方值的数，如：2.25开平方为1.5，4开平方为2。这样选取4个点，利用这四个点进行样条插值，最后计算插值多项式在3处的值，进而简化计算（就是不用开平方了，可以做做乘法，计算多项式的值就可以了）。

2.实验程序及注释：
x=[1 2.25 4 6.25];y=[1 1.5 2 2.5]; %设置插值结点
xi=[1:0.01:6.25];f=xi.^0.5; %原函数
yi=spline(x,y,xi); %样条插值
y3=spline(x,y,3) %求在3出的样条插值多项式的值
sqrt3=sqrt(3) %真实值
error=abs(y3-sqrt3) %显示误差
plot(x,y,'o',xi,f,xi,yi,'r')
clear %清除内存变量

3.实验结果及分析
y3
sqrt3
error
1.7365
1.7321
0.0045

注：绿线为原函数，红线为样条插值函数，蓝色圆圈为插值结点。
可以看出样条插值在图形上表现得很好，所以广泛应用于精度较高的图片放大。由图可以看出，样条函数光滑性很好，但没有原函数弯曲幅度大。
样条插值是在分段三次埃尔米特插值法改进基础上得到的，因为只有在被插值函数在插值结点处的函数值和导数值已知时才能使用，但在实际问题中不现实，而样条插值克服了这一点，只需知道插值结点和函数值就可以了，而且这类插值函数在插值结点处具有二阶导数连续，从而具有很好的光滑性。不过缺点也很明显，就是插值条件不够求出所有的未知数，要人工的添加边界条件。