数值计算笔记-非线性方程求解

2.1根本问题

1.问题：我们经常遇到f(x) = 0求解的问题．我们想要知道当f(x)为非线性方程时有什么方法求解f(x)的零点，也就是上面方程的根．在这里会介绍二分法，不动点迭代法，牛顿法，割线法，逆二次插值法，zeroin法．

2.2问题的基本概念和敏感性

1.根的重数：对光滑函数f，若
$f(x^*) = f'(x^*) = \cdots = f^{m - 1}(x^*) = 0 \land f^m(x^*) \neq 0$f(x∗)=f′(x∗)=⋯=fm−1(x∗)=0∧fm(x∗)​=0

则$x^*$x∗为方程的m重根

2.问题的敏感性

求解非线性方程的f(x) = y问题是“已知y = 0，求x = ?”

问题条件数：$cond = |\frac{\Delta x}{\Delta y}| \approx \frac{1}{|f'(x^*)|}$cond=∣ΔyΔx​∣≈∣f′(x∗)∣1​

2.3二分法

1.二分法（初中学的二分法）

输入：有根区间[a, b]，函数f(x)

输出：零点$x^*$x∗

（1）$x_0 = a, k = 0, a_0 = a, b_0 = b$x0​=a,k=0,a0​=a,b0​=b

（2）若$b_k - a_k < \epsilon$bk​−ak​<ϵ，则$x^* = x_k$x∗=xk​结束

（3）若$sign(f(x_k)) = sign(f(a_k))$sign(f(xk​))=sign(f(ak​))，则$a_{k + 1} = x_k, b_{k+1} = b_k$ak+1​=xk​,bk+1​=bk​，否则$a_{k + 1} = a_k, b_{k + 1} = x_k$ak+1​=ak​,bk+1​=xk​

（4）k = k + 1转（2）

2.误差分析，机械精度$\epsilon$ϵ

最终停止时：$b_k - a_k = 2^{\lfloor log_2|x^*|\rfloor} \cdot 2\epsilon$bk​−ak​=2⌊log2​∣x∗∣⌋⋅2ϵ

绝对误差：$e(x_k) = |x_k - x^*| \leq 2^{\lfloor log_2|x^*|\rfloor} \cdot 2\epsilon$e(xk​)=∣xk​−x∗∣≤2⌊log2​∣x∗∣⌋⋅2ϵ

相对误差：$e_r(x_k) = \frac{|x_k - x^*|}{|x^*|} \leq 2\epsilon$er​(xk​)=∣x∗∣∣xk​−x∗∣​≤2ϵ

2.4不动点迭代法

1.思路：将f(x) = 0改为$x = \phi(x)$x=ϕ(x)，构造迭代公式$x_{k + 1} = \phi(x_k)$xk+1​=ϕ(xk​)，下图给出了迭代的图形解释

注：确实有可能迭代得结果越来偏离正确解，之后会用特定的算法约束迭代过程．

2.不动点迭代法

输入：函数f(x)，初值$x_0$x0​，迭代公式$x = \phi(x)$x=ϕ(x)

输出：零点$x^*$x∗

（1）初值$x_0 = x_0, k = 0$x0​=x0​,k=0

（2）计算$x_{k + 1} = \phi(x_{k})$xk+1​=ϕ(xk​)

（3）若$|f(x_k)| < \epsilon_1 \land |x_{k + 1} - x_{k}| < \epsilon_2$∣f(xk​)∣<ϵ1​∧∣xk+1​−xk​∣<ϵ2​则$x^* = x_{k + 1}$x∗=xk+1​结束，否则转（2）

3.收敛定理：若$\phi(x) \in C[a, b]$ϕ(x)∈C[a,b]满足，C[a, b]表示[a, b]上的连续函数

（1）$\forall x\in [a, b], a \leq \phi(x) \leq b$∀x∈[a,b],a≤ϕ(x)≤b

（2）$\exist L \in (0, 1), 对 \forall x_1, x_2 \in[a, b]$∃L∈(0,1),对∀x1​,x2​∈[a,b]有$\phi(x_1) - \phi(x_2) \leq L|x_1 - x_2|$ϕ(x1​)−ϕ(x2​)≤L∣x1​−x2​∣

则$\phi(x)$ϕ(x)在[a, b]上有唯一不动点，称为全局收敛

证：若$\phi(a) = a$ϕ(a)=a或$\phi(b) = b$ϕ(b)=b则a或b即为不动点

否则$\phi(a) > a \land \phi(b) < b$ϕ(a)>a∧ϕ(b)<b

令$f(x) = \phi(x) - x$f(x)=ϕ(x)−x则f(x) < 0, f(b) > 0

f(x)为连续函数，因此$\exist x^* \in [a, b]$∃x∗∈[a,b]

定理：若将收敛定理中的（2）改为$|\phi'(x)| < 1$∣ϕ′(x)∣<1，则同样有全局收敛性

4.局部收敛：函数$\phi(x)$ϕ(x)存在不动点$x^*$x∗，若存在其某个邻域$D:[x^* -\delta, x^* + \delta]$D:[x∗−δ,x∗+δ]，对于任意初值$x_0 \in D$x0​∈D，迭代法$x_{k + 1} = \phi(x_k)$xk+1​=ϕ(xk​)产生的解序列$\{x_k\}${xk​}收敛到$x^*$x∗，则迭代法局部收敛

定理：设$x^*$x∗为函数$\phi(x)$ϕ(x)的不动点，若$\phi'(x) < 1, \forall x \in [x^* - \delta, x^* + \delta]$ϕ′(x)<1,∀x∈[x∗−δ,x∗+δ]，则迭代法局部收敛

5.收敛阶：迭代解序列$\{x_0, x_1, \cdots, x_k, \cdots \}${x0​,x1​,⋯,xk​,⋯}收敛于$x^*$x∗，误差$e(x_k) = x_k - x^*$e(xk​)=xk​−x∗
$\lim_{k \rightarrow \infty} \frac{|e(x_{k + 1})|}{|e(x_k)|^p} = C \neq 0$k→∞lim​∣e(xk​)∣p∣e(xk+1​)∣​=C​=0
则收敛阶为p，其中C为常数

p阶收敛定理：迭代公式$x_{k + 1} = \phi(x_k)$xk+1​=ϕ(xk​)，若$\phi(x)$ϕ(x)在$x^*$x∗邻域上p阶导数连续$p \geq 2$p≥2，则该迭代法p阶收敛的充要条件为$\phi'(x^*) = \cdots = \phi^{p - 1}(x^*) = 0$ϕ′(x∗)=⋯=ϕp−1(x∗)=0

证明思路：
$e(x_{k + 1}) = x_{k + 1} - x^* = \phi(x_k) - \phi(x^*)\\ = \phi'(x^*)(x_k - x^*) + \cdots + \frac{1}{p!}\phi^{(p)}(\xi)(x_k - x^*)^p$e(xk+1​)=xk+1​−x∗=ϕ(xk​)−ϕ(x∗)=ϕ′(x∗)(xk​−x∗)+⋯+p!1​ϕ(p)(ξ)(xk​−x∗)p
在此泰勒公式的基础上再带入定义化简

2.5牛顿法

1.思路：$f(x) = 0 \Rightarrow f(x) \approx f(x_k) + (x - x_k)f'(x_k) \approx 0$f(x)=0⇒f(x)≈f(xk​)+(x−xk​)f′(xk​)≈0

迭代公式即为：$x_{k + 1} = x_k - \frac{f(x_k)}{f'(x_k)}$xk+1​=xk​−f′(xk​)f(xk​)​（牛顿法也是不动点迭代法）

2.牛顿法

输入：函数f(x) = 0，初值$x_0$x0​，迭代公式$x_{k + 1} = x_k - \frac{f(x_k)}{f'(x_k)}$xk+1​=xk​−f′(xk​)f(xk​)​

输出：输出最优值$x^*$x∗

（1）初值$x_0 = x_0$x0​=x0​，k = 0

（2）计算$x_{k + 1} = x_{k} - \frac{f(x_k)}{f'(x_k)}$xk+1​=xk​−f′(xk​)f(xk​)​

（3）若$|f(x_k)| < \epsilon_1 \land |x_{k + 1} - x_{k}| < \epsilon_2$∣f(xk​)∣<ϵ1​∧∣xk+1​−xk​∣<ϵ2​，则$x^* = x_{k + 1}$x∗=xk+1​结束

（4）k = k + 1转（2）

3.单根收敛定理：设$x^*$x∗是f(x) = 0的单跟，且f(x)在$x^*$x∗附近有2阶连续导数，则牛顿法产生的解序列2阶收敛

4.牛顿法m重根收敛定理：$\lim_{k\rightarrow \infty} |\frac{e(x_{k + 1})}{e(x_k)}| = 1 - \frac{1}{m}$limk→∞​∣e(xk​)e(xk+1​)​∣=1−m1​，其中m为非线性方程根的重数．

2.6割线法与抛物线法

1.割线法：在牛顿法迭代时，用差商代替导数$f'(x) \approx \frac{f(x_k) - f(x_{k - 1})}{x_k - x_{k - 1}}$f′(x)≈xk​−xk−1​f(xk​)−f(xk−1​)​，即$(x_k, y_k),(x_{k - 1}, y_{k - 1})$(xk​,yk​),(xk−1​,yk−1​)的斜率

2.抛物线法（二次插值法）：带入$x_k, x_{k - 1}, x_{k - 2}$xk​,xk−1​,xk−2​算出抛物线，用抛物线与x轴交点作为$x_{k + 1}$xk+1​

注：抛物线有可能不与x轴相交，可以使用“侧向抛物线”把xy轴较换求抛物线，这样一定会与x轴相交（逆二次插值法）

2.7使用的方程求根基数

1.阻尼牛顿法：将牛顿法的迭代公式变为$x_{k + 1} = x_{k} - \lambda_i\frac{f(x_{k})}{f'(x_{k})}$xk+1​=xk​−λi​f′(xk​)f(xk​)​

其中$\lambda_i \in \{\lambda\} \subset (0, 1]$λi​∈{λ}⊂(0,1]，每次选择的$\lambda_i$λi​保证$|f(x_{k + 1})| \leq|f(x_k)|$∣f(xk+1​)∣≤∣f(xk​)∣

2.zeroin通法

思路：综合使用逆二次插值法，割线法，二分法得到收敛速度快，稳定的算法

输入：方程f(x) = 0，有根区间[a, b]

输出：解$x^*$x∗

（1）选取初始值a和b，使得f(a)和f(b)异号

（2）将a赋值给c

（3）重复循环直到满足精度$|f(b) \leq \epsilon_1 \lor |a - b| \leq \epsilon_2|b|$∣f(b)≤ϵ1​∨∣a−b∣≤ϵ2​∣b∣

a．若f(b)的正负号与f(a)相同，a = c；

b. 若|f(a)| < |f(b)|，则c = b，然后对调a，b；

c. 如果$c \neq a$c​=a，使用a，b，c对应的点作逆二次插值法，得到与x轴交点赋值给b；否则执行割线法，得到与x轴交点赋值给b

d. 若得到的近似解“比较满意”，则赋值给b不变，否则使用a，c重新用二分法得到新的b；将3.b中得到的b赋值给c

注：3.a和3.b都是为了调整a，b，c的位置，让|f(b)| < |f(a)|；3.c和3.d则是选择最优的迭代方法．

3.zeroin例子

下面列出zeroin实际执行的每一步，每一步分别用（1），（2），（3），3.a，3.b，3.c，3.d表示

（1）初始化$a_0, b_0$a0​,b0​如上图所示

（2）初始$c_0 = a_0$c0​=a0​

（3）判断精度未达到

3.a f(b)与f(a)符号不同，无操作

3.b |f(a)| < |f(b)|，$c_1 = b_0, b_1 = a_0, a_1 = b_0$c1​=b0​,b1​=a0​,a1​=b0​（到这里已经调整好a，b，c的位置）

3.c $c_1 = a_1$c1​=a1​得到做割线法，得到上图中$b_2$b2​

3.d 得到的结果满意，$c_2 = b_1$c2​=b1​，为了说明方便，$a_2 = a_1$a2​=a1​

（3）判断精度未达到

3.a 无操作

3.b 无操作

3.c $c_2 \neq a_2$c2​​=a2​，则根据逆二次插值法得到$b_3$b3​如上图所示

3.d 得到的结果满意，$c_3 = b_2, a_3 = a_2$c3​=b2​,a3​=a2​

……

注：此处只是为了说明方便才列出所有的角标，实际程序运行过程中都是同一个变量，而不是数组