第 3 章：非线性方程求解（二）

3.4Newton迭代法

对非线性方程\(f(x)=0\)可构造多种迭代格式\(x_{k+1}=φ(x_k)\)，\(Newton\)迭代法是对函数\(f(x)\)做\(Taylor\)展开取其线性部分构造的一种迭代格式

作\(f(x)\)在\(x_0\)处的泰勒展开：

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

取展开式的线性部分作为\(f(x_0)=0\)的近似方程，则有

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

设\(f'(x_0) ≠ 0\)，则

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

令

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

类似的，再作\(f(x)\)在\(x_1\)的泰勒展开，并取其线性部分得到

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

一直做下去得到Newton迭代格式：

\[x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}, \quad k = 0, 1, \ldots \]

牛顿迭代格式对应于\(f(x)=0\)的等价方程是

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

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

若\(x^*\)是\(f(x)\)的单根时，\(f(x^*)=0,f'(x^*)≠0\)，则有\(|φ'(x^*) = 0|\)，只要初值\(x_0\)充分接近\(x^*\)，牛顿迭代收敛，且对于充分大的\(k\),

\[|x_{k+1} - x^*| \approx |x_k - x^*|^2 \frac{|\varphi''(x^*)|}{2} \]

因此牛顿迭代是平方收敛的

牛顿迭代法的几何意义

以\(f'(x_0)\)为斜率作过\((x_0,f(x_0))\)点的直线，即作\(f(x)\)在点\(x_0\)的切线

\[y - f(x_0) = f'(x_0)(x - x_0) \]

令\(y=0\),则得此切线与\(x\)轴的交点\(x_1\)，即

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

再作\(f(x)\)在\(x_1\)处的切线，得交点\(x_2\)，逐步逼近方程的根\(x^*\)，如图所示。因此牛顿法也称为切线法

例题

\(用牛顿法解方程xe^x-1=0\)

牛顿迭代法的局限性

在牛顿迭代法中收敛与否与迭代值\(x_0\)密切相关，当迭代初始值在某根附近时，迭代才能收敛到这个根，有时会发生从一个根附近跳跃到另一个根附近的情况，尤其在导数\(f'(x)\)很小时

如果\(f(x)=0\)没有实根，初始值\(x_0\)是实数，则迭代序列不收敛。下图给出迭代函数\(f(x)=2+x^2\),初始值\(x_0=2\)的发散的迭代序列

3.4弦截法和抛物线法

弦截法

设\(x_k,x_{k+1}\)是\(f(x)=0\)的近似根，我们利用\(f(x_k),f(x_{k+1})\)构造一次插值多项式\(p_1(x)\)，并用\(p_1(x)=0\)的根作为\(f(x)=0\)的新的近似根\(x_{k+1}\)
由于

\[p_1(x) = f(x_k) + \frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}}(x - x_k) \]

因此有

\[x_{k+1} = x_k - \frac{f(x_k)}{f(x_k) - f(x_{k-1})}(x_k - x_{k-1}) \\ = x_k - \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})}f(x_k) \]

这样导出的迭代公式可以看做牛顿公式

\[x_{k+1} = x_k + \frac{f(x_k)}{f'(x_k)} \]

中导数\(f'(x)\)用差商

\[\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} \]

代替的结果
作过\((x_{k-1},f(x_{k-1}))\)和\((x_{k},f(x_{k}))\)的一条直线（弦），该直线的斜率等于差商值\(\frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}}\)，该直线与\(x\)轴的交点的横坐标就是生成的迭代点\(x_{k+1}\)。这种算法因此被称为弦截法

例题

用弦截法解方程

\[f(x)=xe^x-1=0 \]

弦截法的收敛速度

抛物线法

设已知方程\(f(x)=0\)的三个近似根\(x_k,x_{k-1},x_{k-2}\),我们以这三个点为节点，构造二次插值多项式\(p_2(x)\)，并适当地选取\(p_2(x)\)的一个零点\(x_{k+1}\)作为新的近似根，这样确定的迭代过程称为抛物线法，亦称为密勒法。从几何图形上看，这种方法是用抛物线\(y=p_2(x)\)与\(x\)轴的交点\(x{k+1}\)作为所求根\(x^*\)的近似位置

插值多项式

\[p_2(x) = f(x_k) + f[x_k, x_{k-1}](x - x_k) + f[x_k, x_{k-1}, x_{k-2}](x - x_k)(x - x_{k-1}) \]

有两个零点

\[x_{k+1} = x_k - \frac{2f(x_k)}{\omega \pm \sqrt{\omega^2 - 4f(x_k)f[x_k, x_{k-1}, x_{k-2}]}} \]

其中

\[\omega = f[x_k, x_{k-1}] + f[x_k, x_{k-1}, x_{k-2}](x_k - x_{k-1}) \]

为了确定哪一个值作为\(x_{k+1}\)，需要讨论根式前正负号的取值

计算表明，抛物线法比弦截法收敛更快
在一定条件下可以证明，对于抛物线法，迭代误差有下列渐进关系式：

\[\frac{|e_{k+1}|}{|e_k|^{1.840}} \rightarrow \left| \frac{f'''(x^*)}{6 f'(x^*)} \right|^{0.42} \]

课件抛物线法是超线性收敛的，其收敛的阶
\(p=1.840\)（是方程\(\lambda^3 - \lambda^2 - \lambda - 1 = 0\)的根），收敛速度比弦截法更接近牛顿法

3.5解线性方程组的牛顿法

为叙述简单，我们以解二阶非线性方程组为例，演示牛顿法迭代求解的方法和步骤。类似地可以得到求高阶非线性方程组的方法和步骤
设二阶方程组

\[\begin{cases} f_1(x, y) = 0 \\ f_2(x, y) = 0 \end{cases}\]

\(x,y\)为自变量。为了方便起见，将方程组写成向量形式\(F(\omega)=0\)，其中

\[F(\omega) = \begin{pmatrix} f_1(x, y) \\ f_2(x, y) \end{pmatrix}, \quad \omega = \begin{pmatrix} x \\ y \end{pmatrix} \]

作\(f_1(x,y),f_2(x,y)在(x_0,y_0)\)的二元泰勒展开，并取其线性部分，得到方程组

\[\begin{cases} f_1(x_0, y_0) + (x - x_0)\frac{\partial f_1(x_0, y_0)}{\partial x} + (y - y_0)\frac{\partial f_1(x_0, y_0)}{\partial y} = 0 \\ f_2(x_0, y_0) + (x - x_0)\frac{\partial f_2(x_0, y_0)}{\partial x} + (y - y_0)\frac{\partial f_2(x_0, y_0)}{\partial y} = 0 \end{cases}\]

令\(x - x_0 = \Delta x, \quad y - y_0 = \Delta y\),则有

\[\begin{cases} \Delta x \frac{\partial f_1(x_0, y_0)}{\partial x} + \Delta y \frac{\partial f_1(x_0, y_0)}{\partial y} = -f_1(x_0, y_0) \\ \Delta x \frac{\partial f_2(x_0, y_0)}{\partial x} + \Delta y \frac{\partial f_2(x_0, y_0)}{\partial y} = -f_2(x_0, y_0) \end{cases}\]

如果

\[\text{det}(J(x_0, y_0)) = \left| \begin{matrix} \frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y} \\ \frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y} \end{matrix} \right|_{(x_0, y_0)} \neq 0 \]

解出\(\Delta x,\Delta y\),得

\[\omega_1 = \omega_0 + \begin{pmatrix} \Delta x \\ \Delta y \end{pmatrix} = \begin{pmatrix} x_0 + \Delta x \\ y_0 + \Delta y \end{pmatrix} = \begin{pmatrix} x_1 \\ y_1 \end{pmatrix} \]

继续做下去，每一次迭代都是解一个方程组

\[J(x_k, y_k) \begin{pmatrix} \Delta x \\ \Delta y \end{pmatrix} = \begin{pmatrix} -f_1(x_k, y_k) \\ -f_1(x_k, y_k) \end{pmatrix} \]

\[\Delta x = x_{k+1} - x_k, \quad \Delta y = y_{k+1} - y_k \]

即

\[x_{k+1} = x_k + \Delta x, \quad y_{k+1} = y_k + \Delta y \]

直到

\[\max(|\Delta x|, |\Delta y|) < \varepsilon \]

为止