迭代格式

非线性方程的数值方法

张瑞
中国科学技术大学数学科学学院

迭代格式

二分法

若函数$f(x)$在$[a,b]$上连续，且$f(a)f(b)<0$，则$f(x)$在$(a,b)$内至少有一个零点。

利用这个定理，可以给出二分法求函数$f(x)$的根。

输入区间a,b，函数f(x)，及误差控制e

while (|a-b|>e) do
  x=(a+b)/2
  if (|f(x)|<e)  return x; // 找到根为x
  else if (f(a)*f(x)<0) b=x; // 根在区间[a,(a+b)/2]
  else if (f(b)*f(x)<0) a=x; // 根在区间[(a+b)/2,b]
end while
return (a+b)/2 // 还没有找到根，但b-a已经足够小了

例 1. 求函数$f(x)=xe^x-1$在$[0.5,0.8]$之间的根

解. 用Python计算，得到 

root =  0.5671432822942734

ex-iterator-bisect

Total  25  steps
   1 0.65 0.15000000000000002
   2 0.575 0.07500000000000001
   3 0.5375 0.03749999999999998
   4 0.5562499999999999 0.01874999999999999
   5 0.5656249999999999 0.009375000000000022
   6 0.5703125 0.004687500000000011
   7 0.5679687499999999 0.0023437500000000333
   8 0.5667968749999999 0.001171874999999989
   9 0.5673828124999999 0.0005859375000000222
  10 0.5670898437499998 0.0002929687500000111
  11 0.5672363281249999 0.0001464843750000333
  12 0.5671630859374999 7.324218750004441e-05
  13 0.5671264648437498 3.6621093750022204e-05
  14 0.5671447753906249 1.8310546875011102e-05
  15 0.5671356201171873 9.155273437533307e-06
  16 0.5671401977539061 4.577636718794409e-06
  17 0.5671424865722655 2.2888183593972045e-06
  18 0.5671436309814453 1.1444091796986022e-06
  19 0.5671430587768553 5.722045898770567e-07
  20 0.5671433448791503 2.861022949662839e-07
  21 0.5671432018280028 1.4305114748314196e-07
  22 0.5671432733535766 7.152557374157098e-08
  23 0.5671433091163635 3.5762786843029915e-08
  24 0.5671432912349701 1.7881393421514957e-08
  25 0.5671432822942734 8.940696738513054e-09

~~

  • 二分法的算法简单
  • 二分法需要保证初始时$f(a)$, $f(b)$异号
  • 二分法的收敛速度较慢，需要的计算步数$n$满足
    \[(b-a)(\frac12)^n\leq\epsilon \, \Rightarrow \, n\geq\frac{|\ln\frac{\epsilon}{b-a}|}{\ln 2} \]

简单迭代

函数$f(x)=0$的等价形态$x=\phi(x)$有很多种，如

\[x=\phi(x)=x+af(x) \]

其中$a$是一个非零常数。

可以看到

\[\phi'(x)=1+a f'(x) \]

因而，取合适的$a$使得$\phi'(x)$在$f(x)$的根$\alpha$附近满足

\[|\phi'(x)|=|1+a f'(x)|\leq L<1 \]

则迭代格式收敛。

特别，若有$a=-\frac1{f'(\alpha)}$，则

\[\phi'(\alpha)=0 \]

此时，格式达到2阶收敛。

但这种情况很难发生，在实际应用中，可以取

\[a=-\frac1{f'(x_0)} \]

Newton 迭代

函数$f(x)$在$x_n$处做Taylor展开，

\[f(x)=f'(x_n)+f'(x_n)(x-x_n)+\frac{f''(x_n)}{2!}(x-x_n)^2+\cdots \]

取线性部分近似为$f(x)$，然后求根

\[f(x_n)+f'(x_n)(x-x_n)\approx 0 \]

得到

\[x=x_n-\frac{f(x_n)}{f'(x_n)} \]

取这个值为下一个迭代值，得到

定义 1.
Newton迭代格式为

\[x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)} \]

iterator-newton

可以看到，直线

\[y=f(x_n)+f'(x_n)(x-x_n) \]

是$f(x)$在$x_n$处的切线，与$x$轴的交点就是下一步迭代的值。

因而，Newton迭代又称为切线法。

Newton 迭代的收敛性

由Newton迭代格式，可以得到迭代格式对应的等价方程是

\[x=x-\frac{f(x)}{f'(x)} \]

记

\[\phi(x)=x-\frac{f(x)}{f'(x)} \]

则有

\[\phi'(x)=\frac{f(x)f''(x)}{(f'(x))^2} \]

设$f(\alpha)=0$。若$f'(\alpha)\neq 0$，则有

\[\phi'(\alpha)=\frac{f(\alpha)f''(\alpha)}{(f'(\alpha))^2}=0 \]

因而有

定理 1.
Newton迭代格式在$f(x)$的单根附近收敛，且至少2阶收敛。

利用前面的定理即可

或者

证明. $f(\alpha)$在$x_n$处Taylor展开

\[0=f(\alpha)=f(x_n)+f'(x_n)(\alpha-x_n)+\frac{f''(\xi_n)}{2!}(\alpha-x_n)^2 \]

其中$\xi_n$在$\alpha$与$x_n$之间。移项，可得

\[-f'(x_n)(\alpha-x_n)-f(x_n)=\frac{f''(\xi_n)}{2!}(\alpha-x_n)^2 \]

进而

\[-\alpha+x_n-\frac{f(x_n)}{f'(x_n)}=\frac{f''(\xi_n)}{2f'(x_n)}(\alpha-x_n)^2 \]
\[-\alpha+x_n-\frac{f(x_n)}{f'(x_n)}=\frac{f''(\xi_n)}{2f'(x_n)}(\alpha-x_n)^2 \]

注意到$\displaystyle x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$

\[x_{n+1}-\alpha=\frac{f''(\xi_n)}{2f'(x_n)}(\alpha-x_n)^2 \]

因此，若$\alpha$是$f(x)$的单根，

\[\frac{|e_{n+1}|}{|e^n|^2}=\frac{f''(\xi_n)}{2f'(x_n)} \to \frac{f''(\alpha)}{2f'(\alpha)}, n\to\infty \]

可以看到，若$\alpha$是$f(x)$的单根，则格式至少2阶收敛。

定理 2.
若$\alpha$是$f(x)$的$m\geq2$重根，则Newton迭代

\[x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)} \]

是1阶收敛的。而格式

\[x_{n+1}=x_n-m\frac{f(x_n)}{f'(x_n)} \]

至少2阶收敛。

此时，$f(x)$可以写成

\[f(x)=(x-\alpha)^m h(x) \]

其中$h(\alpha)\neq 0$。

则Newton迭代表达为

\[\begin{aligned} x_{n+1}=&x_n-\frac{f(x_n)}{f'(x_n)} \\ =&x_n-\frac{(x_n-\alpha)^mh(x_n)}{m(x_n-\alpha)^{m-1}h(x_n)+(x_n-\alpha)^mh'(x_n)} \\ =&x_n-\frac{(x_n-\alpha)h(x_n)}{mh(x_n)+(x_n-\alpha)h'(x_n)} \\ \end{aligned} \]

因而

\[\begin{aligned} x_{n+1}-\alpha =&x_n-\alpha-\frac{(x_n-\alpha)h(x_n)}{mh(x_n)+(x_n-\alpha)h'(x_n)} \\ =&(x_n-\alpha)\left(1-\frac{h(x_n)}{mh(x_n)+(x_n-\alpha)h'(x_n)}\right) \end{aligned} \]

得到

\[\begin{aligned} \frac{e_{n+1}}{e_n}=&1-\frac{h(x_n)}{mh(x_n)+(x_n-\alpha)h'(x_n)} \\ \to & 1-\frac{h(\alpha)}{mh(\alpha)+0}=1-\frac1m\neq 0, n\to\infty \end{aligned} \]

即，此时Newton迭代1阶收敛。

若函数$f(x)$的根$\alpha$的重数未知，不防设为$p$，则

\[f(x)=(x-\alpha)^p h(x) \]

令

\[\begin{aligned} F(x)=&\frac{f(x)}{f'(x)} =\frac{(x-\alpha)^ph(x)}{p(x-\alpha)^{p-1}h(x)+(x-\alpha)^ph'(x)} \\ =&(x-\alpha)\frac{h(x)}{ph(x)+(x-\alpha)h'(x)} \end{aligned} \]

则$\alpha$是$F(x)$的一重根。

可以构造$F(x)=0$的Newton迭代格式。

若$\alpha$是函数$f(x)$重根，则解方程$f(x)=0$的格式

\[x_{n+1}=x_n-\frac{F(x_n)}{F'(x_n)} =x_n-\frac{f(x_n)f'(x_n)}{(f'(x_n))^2-f(x_n)f''(x_n)} \]

具有至少2阶精度。

Newton迭代的伪代码

设置f(x)计算函数值，df(x)计算导数值，tol为精度控制，maxiter为最大迭代步数， m为预设的重根数，缺省为1， 

    while k <= maxiter and abs( x - x0 ) >= tol:
        x2, x1, x0 = ( x1, x0, x )
        fx = f(x)
        dfx = df(x)
        x = x - m * fx / dfx
        k = k + 1

        if k > maxiter:
            return (x,float('inf'),stepinfo,k)

例 2. Newton迭代求$f(x)=e^x-1-x=0$的根。

解. 可以看到$f(0)=0$且$f'(0)=1-e^0=0$。因而，$0$是一个二重根。

以初值$x_0=1$开始计算，误差控制$10^{-5}$。

使用Newton迭代计算 

Root is  5.424952541628956e-06
Total steps: 18
Rate: 1.000010240574523

以同样的初值，设置重根阶$m=2$ 

Root is  1.0864531688955798e-11
Total steps: 4
Rate: 2.0147483986450294

未知根重数的Newton迭代 

取初值 1
Root is  -4.218590698935789e-11
Total steps: 5
Rate: math domain error
  1 -0.23421061355351425
  2 -0.00845827991076109
  3 -1.1890183808588653e-05
  4 -4.218590698935789e-11
  5 -4.218590698935789e-11

Newton迭代的中间结果 

  1 0.5819767068693265
  2 0.31905504091081843
  3 0.16799617288577048
  4 0.08634887374778137
  5 0.04379570367371408
  6 0.022057685365768236
  7 0.0110693874777393
  8 0.005544904662931229
  9 0.0027750144941372577
 10 0.0013881489723892668
 11 0.0006942350659796617
 12 0.0003471576966651309
 13 0.00017358889159729097
 14 8.679695657422037e-05
 15 4.339910703392076e-05
 16 2.1699709854160184e-05
 17 1.0849887297322925e-05
 18 5.424952541628956e-06

重根Newton迭代中间结果 

  1 0.1639534137386529
  2 0.0044781144487033575
  3 3.342250383920123e-06
  4 1.0864531688955798e-11

迭代算法的稳定性

  • 可以证明Newton迭代格式在根附近是收敛的。
  • 由于非线性问题的复杂性，在离根稍远处的特性未知， 因而迭代序列的收敛性是与初值的选取相关的。
    • 如果初值合适，经过若干步后，迭代的值落在收敛范围内，则序列收敛。
    • 如果初值不当，迭代的值始终不能落在收敛范围，则序列发散。
    • 甚至可能，在一个根附近的初值，经过迭代后，收敛到方程的另一个根。

例 3. 用Newton迭代求根$\arctan(x)=0$

解. 用初值$x_0=1$，结果为 

Root is  0.0
Total steps: 5
Rate: 2.9936674514109285
  1 -0.5707963267948966
  2 0.1168599039989131
  3 -0.001061022117044716
  4 7.963096044106416e-10
  5 0.0

作业： 结果中，可以看到3阶收敛，试分析原因。

证明： 求解方程$\arctan(x)=0$的Newton迭代格式具有至少3阶精度。

用初值$x_0=2$，结果为 

Root is  nan
Total steps: 11
Rate: nan
  1 -3.535743588970452
  2 13.95095908692749
  3 -279.3440665336173
  4 122016.99891795448
  5 -23386004197.933853
  6 8.590766671950354e+20
  7 -1.1592676698907246e+42
  8 2.110995587610979e+84
  9 -6.9999433953175654e+168
 10 inf
 11 nan

是一个发散到$\infty$的序列。

例 4. 用Newton迭代法解方程

\[x^3-2x+2=0 \]

解. 构造Newton迭代格式为

\[x_{n+1}=x_n-\frac{x_n^3-2x_n+2}{3x_n^2-2} =\frac{2x_n^3-2}{3x^2_n-2} \]

取初值$x_0=0$，则有

\[x_0=0 \, \Rightarrow\, x_1=1 \, \Rightarrow\, x_2=0 \, \Rightarrow\, x_3=1 \]

迭代序列在$0,1$中循环取值，永远不会收敛。

例 5. Newton迭代求根$x^3-2x^2-11x+12=0$

解. 用不同的初值，得到结果为 

取初值 2.35283735
Root is  4.000000000000001
Total steps: 25
Rate: 2.001075621801667
取初值 2.352836327
Root is  -3.0000000000000004
Total steps: 25
Rate: 2.000886478517893
取初值 2.352836323
Root is  0.9999999999999925
Total steps: 16
Rate: 1.89270344136726
取初值 2.35283735
Root is  4.000000000000001
Total steps: 25
Rate: 2.001075621801667
  1 -0.7829479432088857
  2 2.3528750599018755
  3 -0.7832624903449901
  4 2.3542988581341184
  5 -0.7951930373389495
  6 2.4096178613487207
  7 -1.357025522098792
  8 436.8325543155845
  9 291.450201808276
 10 194.53176652201245
 11 129.92416933967377
 12 86.85946275830207
 13 58.160163230315774
 14 39.04298283697368
 15 26.32157018154881
 16 17.8753619404899
 17 12.295986835648872
 18 8.652254263160414
 19 6.334490245199713
 20 4.951264807478875
 21 4.252004042721794
 22 4.025430907342878
 23 4.0003021866887325
 24 4.000000043474304
 25 4.000000000000001
取初值 2.352836327
Root is  -3.0000000000000004
Total steps: 25
Rate: 2.000886478517893
  1 -0.7829394111557035
  2 2.3528364638374564
  3 -0.7829405524081285
  4 2.352841626395325
  5 -0.7829836099095435
  6 2.3530364176331253
  7 -0.78460924829194
  8 2.3604147342051243
  9 -0.8476711759310573
 10 2.6872287516836755
 11 -144.9559225012978
 12 -96.43394799584374
 13 -64.09545863227268
 14 -42.550769282273606
 15 -28.209243369124543
 16 -18.680961503769232
 17 -12.37865102649303
 18 -8.253684045483189
 19 -5.62221452244607
 20 -4.050604387537867
 21 -3.2657018020731914
 22 -3.0239035181209393
 23 -3.0002212761583498
 24 -3.0000000192329486
 25 -3.0000000000000004
取初值 2.352836323
Root is  0.9999999999999925
Total steps: 16
Rate: 1.89270344136726
  1 -0.7829393777949023
  2 2.3528363129272205
  3 -0.7829392937858954
  4 2.352835932905887
  5 -0.7829361243354871
  6 2.352821595739119
  7 -0.7828165551125976
  8 2.3522808476279304
  9 -0.7783146062709512
 10 2.3321026642860514
 11 -0.6205467552357247
 12 1.7993798401408654
 13 0.8042684705782203
 14 0.9981009654022841
 15 0.9999997007081824
 16 0.9999999999999925

弦截法

在Newton迭代中，$f'(x_n)$的计算需要解析求解，有时候不是很方便。

使用差商来近似$f'(x_n)$，得到弦截法

\[x_{n+1}=x_n-(x_n-x_{n-1})\frac{f(x_n)}{f(x_n)-f(x_{n-1})} \]

节点$x_{n+1}$正好是连接$(x_{n-1}, f(x_{n-1}))$和$(x_{n}, f(x_{n}))$的弦 与$x$轴的交点。

num-iterator-scant

  • 相对于Newton迭代来说，弦截法不需要导数的显示表达式，计算上相对简单
  • 弦截法需要2个初值才可以开始计算。
  • 弦截法下一步迭代$x_{n+1}$的值除了跟$x_n$有关外，还跟前一步的$x_{n-1}$相关。 称这种格式为多步格式（或2步格式）。 而Newton迭代中，下一步迭代$x_{n+1}$的值只与$x_n$相关的格式称为单步格式。
  • 弦截法的收敛阶要低于Newton迭代，为1.618。

从插值近似的角度来看，Newton迭代和弦截法都是用一个近似多项式的根来作为函数$f(x)$的根的估计

  • Newton迭代是用满足 $f(x_n)$, $f'(x_n)$的Hermite型插值多项式
  • 弦截法是用满足 $f(x_{n-1})$, $f(x_n)$的线性插值多项式

问题. 可以用3个点，如 $f(x_{n-2})$, $f(x_{n-1})$, $f(x_n)$做2次插值多项式来近似吗？

可以。这样的格式称为抛物线法。

例 6. 求根$f(x)=1-xe^x=0$

解. 取误差控制$\epsilon=10^{-8}$,

二分法， 取初始区间为$[0,2]$ 

root =  0.5671432837843895
Estimated Convergence Rate =  1.00
Total  28  steps

Newton 迭代，取初值$x_0=1$ 

root =  0.567143290409784
Estimated Convergence Rate =  2.00
Total  5  steps

弦截法 

root =  0.5671432904097838
Estimated Convergence Rate =  1.60
Total  7  steps

简单迭代，初值$x_0=1$ 

root =  0.5671432904097838
Estimated Convergence Rate =  2.00
Total  6  steps

程序作业

用Newton迭代和弦截法求根，

\[\frac{x^3}3-x=0 \]

误差控制为$\epsilon=10^{-8}$。

  1. 取初值为$0.1$, $0.2$, $0.9$, $9.0$，给出Newton迭代法的根和迭代步数
  2. 取初值对$(0.1,0.2)$, $(0.2,0.9)$, $(8.0,9.0)$ 给出弦截法的根和迭代步数

输出示例： 

Newton迭代，初值、根和迭代步数：
初值 0.1 , XXXXXXXXXXXXXXXXX , xxxxx
初值 0.2 , XXXXXXXXXXXXXXXXX , xxxxx
初值 0.9 , XXXXXXXXXXXXXXXXX , xxxxx
初值 9.0 , XXXXXXXXXXXXXXXXX , xxxxx

弦截法，初值、根和迭代步数：
初值 0.1,0.2 , XXXXXXXXXXXXXXXXX , xxxxx
初值 0.2,0.9 , XXXXXXXXXXXXXXXXX , xxxxx
初值 8.0,9.0 , XXXXXXXXXXXXXXXXX , xxxxx

方程组的Newton迭代

简记线性方程组

\[\begin{cases} f_1(x_1,x_2,\cdots,x_n)=0 \\ f_2(x_1,x_2,\cdots,x_n)=0 \\ \vdots \\ f_n(x_1,x_2,\cdots,x_n)=0 \\ \end{cases} \]

为$\vec F(\vec x)=0$，其中

\[\begin{aligned} \vec x=&(x_1,x_2, \cdots, x_n)^T, \\ \vec F(\vec x)=&(f_1(x), f_2(x),\cdots, f_n(x))^T \end{aligned} \]

则有

\[\frac{d\vec F(\vec x)}{d \vec x} =\begin{pmatrix} \frac{\partial f_1(\vec x)}{\partial x_1} & \frac{\partial f_1(\vec x)}{\partial x_2} & \cdots & \frac{\partial f_1(\vec x)}{\partial x_n} \\ \frac{\partial f_2(\vec x)}{\partial x_1} & \frac{\partial f_2(\vec x)}{\partial x_2} & \cdots & \frac{\partial f_2(\vec x)}{\partial x_n} \\ \vdots & \vdots & & \vdots \\ \frac{\partial f_n(\vec x)}{\partial x_1} & \frac{\partial f_n(\vec x)}{\partial x_2} & \cdots & \frac{\partial f_n(\vec x)}{\partial x_n} \\ \end{pmatrix} =J_F(\vec x) \]

其中的矩阵为Jacobi矩阵。

因而，函数$\vec F(\vec x)$的在点$\vec \alpha$处的线性近似为

\[\vec F(\vec x)\approx \vec F(\vec \alpha)+J_F(\vec x_n)(\vec x-\vec \alpha) \]

把线性近似的根作为下一个迭代值，就得到了

定义 2.
记解非线性方程组$\vec F(\vec x)=0$的迭代序列为$\{\vec X^{(n)}\}$， 则Newton迭代格式为

\[\vec X^{(n+1)}=\vec X^{(n)}-(J_F(\vec X^{(n)}))^{-1}\vec F(\vec X^{(n)}) \]

为避免求解Jacobi矩阵的逆矩阵，通常将格式改为

\[J_F(\vec X^{(n)})\Delta \vec X=-F(\vec X^{(n)}), \vec X^{(n+1)}=\vec X^{(n)}+\Delta \vec X \]

例 7. (例3.6) 解非线性方程组

\[\begin{cases} 4-x^2-y^2 =0 \\ 1-e^x-y =0 \end{cases} \]

解. Jacobi矩阵为

\[\begin{pmatrix} -2x & -2y \\ -e^x & -1 \end{pmatrix} \]

取初值$\vec X^{(0)}=(1,-1.7)$，则

\[J_F(\vec X^{(0)})=\begin{pmatrix} -2 & 3.4 \\ -2.71828 & -1 \end{pmatrix} \]
\[F(\vec X^{(0)})=\begin{pmatrix} 0.11 \\ -0.01828 \end{pmatrix} \]

解线性方程组

\[\begin{pmatrix} -2 & 3.4 \\ -2.71828 & -1 \end{pmatrix} \begin{pmatrix} dx_1 \\ dx_2 \end{pmatrix} =\begin{pmatrix} 0.11 \\ -0.01828 \end{pmatrix} \]

得到

\[\begin{pmatrix} dx_1 \\ dx_2 \end{pmatrix} =\begin{pmatrix} 0.004256 \\ -0.029849 \end{pmatrix} \]

因而，下一个迭代值为

\[\vec X^{(1)}=\begin{pmatrix} 1 \\ -1.7 \end{pmatrix} +\begin{pmatrix} 0.004256 \\ -0.029849 \end{pmatrix} =\begin{pmatrix} 1.004256 \\ -1.729849 \end{pmatrix} \]

目录

本节读完

例 8.

8. 

************************************************************
Finding root of f(x) = 1 - x * exp(x) on [0.000000,2.000000]
************************************************************

-------- Bisection Method ---------------------------

root =  0.5671432837843895
Estimated Convergence Rate =  1.00
Total  28  steps

-------- Newton's Method ----------------------------

root =  0.567143290409784
Estimated Convergence Rate =  2.00
Total  5  steps

-------- Secant Method ------------------------------

root =  0.5671432904097838
Estimated Convergence Rate =  1.60
Total  7  steps

-------- Fixed-Point Method -------------------------

root =  0.5671432904097838
Estimated Convergence Rate =  2.00
Total  6  steps

左边

右边