Mathematics & Programme

为讨论方便设计算过程中不会刚好取到 f(m)=0

1)  m=(a+b)/2 2)  f(m)f(a)>0吗? 是则令a=m, 否则令 b=m 3)  达到精度要求了吗? 是则输出m, 否则返回(1)

函数曲线的一种情况.

f(m)f(a)>0 这一步不需真的做乘法.

可以想象的出, 对于多数函数 a,b中点是一个不差的选择, 每次可以缩小根区间的一半. 

很容易的, 马上会想到多数常见函数在a-b较小时函数会比较"直", 取过(a,f(a)) (b,f(b))的直线与x轴的交点可能会更好:

直线方程: y=((f(b)-f(a))/(b-a))(x-a)+f(a), y=0 时 x=-f(a)*(b-a)/(f(b)-f(a))+a.

于是有以下算法:

1)  m=-f(a)*(b-a)/(f(b)-f(a))+a 2)  f(m)f(a)>0吗? 是则令a=m, 否则令 b=m 3)  达到精度要求了吗? 是则输出m, 否则返回(1)

不过这一算法有个问题, 由于函数在这一区间很可能只是凹的或凸的, 于是f(m)的符号总是不变, a或b也总是保持原有的值, 影响收敛速度, 根区间基本不再缩小(难以确定根的精度).与二分法交替使用也许会有些改进.

至于牛顿迭代法则是取一个端点的切线, 把切线与x轴的交点作为更好的根的逼近.

切线方程(假设作(b,f(b))的切线): y=f'(b)(x-b)+f(b), y=0 时 x=-f(b)/f'(b)+b. f'(x)是f(x)的导函数, 见微积分.

于是有以下算法:

1)  x=x-f(x)/f'(x) 2)  x不属于[a,b]吗? 是则停止计算(计算失败). 3)  达到精度要求了吗? 是则输出x, 否则返回(1)

举一个实例, 解 sin(3x)-x=0, 已知在(0.5,1)内有根.

y=sin(3x)-x 的函数图像

再举一例, 解 x/ln(x)-20=0, 已知在(60,110)内有根. 图像几乎是一根直线.

可以看到, 二分法收敛速度慢而稳, 割线法对于较直的函数也比较快. 而牛顿迭代法只需几步就完成计算, 非常快. 

 { 生这些数据的程序 root01.pas , 用Turbo Pascal运行时可能要打开编译开关 {$N+} 以便能使用浮点类型double }

定义 设迭代序列{x_k}收敛到x^*,若存在实数p≥1及a≠0, 且当k≥k₀时x_k=x^*, 有

则称序列{x_k}是p阶收敛的. α称为收敛因子.

Newton法  (2.01)

    局部收敛定理: 设f(x^*)=0, f'(x^*)≠0, 且f(x)在x^*附近的二阶导数f''(x)连续, 则(2.01) 具有二阶收敛, 且

 二阶收敛表现在计算得数上就是每算一次会有多几乎一倍的精度.

函数的计算一般总是可以转化为方程求根, 问题是怎样转换才适合高精度计算.

2.1 求倒数

方程 x-1/A=0 显然不行. 应该是 f(x)=1/x-A=0

代入公式 (2.01), f'(x)=-1/x^2. 得:

x=x+(1/x-A)*x^2=x+(1-Ax)x ①  = 2x-Ax^2 ②

可以看到 x 递推式的结构是 x = x + (x的修正值). 并且x的修正值是x的二阶小量: 1-Ax≈0,

计算 1-Ax 时可以省一些计算量, 所以计算的时候会用①式而不是②式.

收敛因子α=-A.

递推式 x = (x+a/x)/2,  α=0.5/√A.

涉及到一次除法.

若是方程 x^-2-A=0

递推式 x = x+(1-Ax^2)x/2,  α=-1.5*√A.

有3次乘法, 并且由于结果是 x=1/√A. 所以计算结束前还要计算一次倒数. 尽管如此, 在高精度计算时仍会用后者而计算器可能用的是前者, 因为3n次乘法和1次除法所花的时间在n和位数较大时比n次除法更短.