数值分析与算法——读书笔记（二）

chapter2

非线性方程求根

2.1引言

线性方程是方程式中仅包含未知量的一次方项和常数项的方程，除此之外的方程都是非线性方程(nonlinear equation)。

定义：对光滑函数f，若f(x∗)=f′(x∗)=⋯=fm−1(x∗)=0，但fm(x∗)≠0，则称x∗为方程的m重根。

2.2二分法

二分法的思想很简单，就是每次将有根区间一分为二，得到长度逐次减半的区间序列{(ak,bk)},则区间中点xk=(ak+bk)/2就是第k步迭代的近似解，具体算法如下：

1. 算法：二分法

2. 输入：a，b，函数f(x)；输出：x。

3. While (b-a)>ε do

   ​ x:=a+(b-a)/2;

   ​ If sign(f(x))=sign(f(a)) then

   ​ a:=x;

   ​ Else

   ​ b:=x;

   ​ End

   End

4. x:=a+(b-a)/2.

假设二分法得到的有根区间序列为{(ak,bk),k=0,1,⋯}，若取解xk=(ak+bk)/2，则误差

|xk−x∗|<(bk−ak)/2=(b0−a0)/2k+1,k=0,1,2,⋯.

二分法求解方程f(x)=x2−2=0：

%%二分法求解f(x)=x^2-2=0
M=2;a=1;b=2;k=0;
while b-a>eps             %%MATLAB中的eps为两倍的机器精度
    x=a+(b-a)/2;
    if x^2>M
        b=x
    else
        a=x
    end
    k=k+1;
end

该程序执行了52次便结束了，最终区间的两个端点已经是两个相邻的浮点数。

1. 二分法是求解单变量方程f(x)=0的实根的一种可靠算法，一定能收敛
2. 二分法解的误差不一定随迭代次数的增加一直减小，在实际的有限精度算术体系中，误差限存在最小值
3. 二分法的缺点是又是不容易确定合适的初始有根区间（含两个初始值）、收敛较慢，且无法求解偶数重的根。因此，实际应用中常将二分法与其他方法结合起来。

2.3不动点迭代法

基本原理

通过某种变换，可将非线性方程f(x)=0改写为x=φ(x)，其中φ(x)为连续函数，给定初始值x0后，可构造迭代计算公式xk+1=φ(xk),(k=0,1,⋯)，从而得到近似解序列{xk}。

由于解x∗满足x∗=φ(x∗)，称它为函数φ(x)的不动点(fixed point)，此方法为求解非线性方程的不动点迭代法(fixed-point iterative method)。

1. 算法：基于函数φ(x)的不动点迭代法

2. 输入：x0，函数f(x)，φ(x);输出：x

3. k:=0;

4. While |f(xk)|>ε1 或 |xk−xk−1|>ε2 do

   ​ xk+1:=φ(xk);

   ​ k:=k+1;

   End

5. x=:xk.

全局收敛的充分条件

设φ(x)∈C[a,b]，若满足如下两个条件：

1. 对任意x∈[a,b]，有a≤φ(x)≤b，

2. 存在正常数L∈(0,1)，使对任意x1,x2∈[a,b]，
   

   |φ(x1)−φ(x2)|≤L|x1−x2|

则φ(x)在[a,b]上存在不动点，且不动点唯一。

定理：设φ(x)∈C[a,b]满足以上两个条件，则对于任意初值x0∈[a,b]，由不动点迭代法得到的序列{xk}收敛到φ(x)的不动点x∗，并有误差估计：

|xk−x∗|≤Lk1−L|x1−x0|

定理：对于不动点迭代法xk+1=φ(xk)，若在所求根x∗的邻域上函数φ(x)的p阶导数连续，p≥2，则该迭代法在x∗的邻域上p阶收敛的充分必要条件是：φ′(x∗)=φ"(x∗)=⋯=φp−1(x∗)=0，且φp(x∗)≠0。

%%不动点迭代法求解f(x)=x^4-x-2=0,x_0=1.5
clear
clc
k=0;xk=1.5;
while abs(xk^4-xk-2)>10*eps
    xk=(xk+2)^(1/4);
    k=k+1;
end
x=xk;

x = 1.353209964199325

k =15

2.4牛顿迭代法

方法原理

设r是f(x)=0的根，选取x0作为r的初始近似值，过点(x0,f(x0))做曲线y=f(x)的切线L，L的方程为y=f(x0)+f′(x0)(x−x0)，求出L与x轴的交点横坐标x1=x0−f(x0)f′(x0)，称x1为r的一次近似值。过点(x1,f(x1))做曲线y=f(x)的切线，并求该切线与x轴交点的横坐标x2=x1−f(x1)f′(x1)，称x2为r的二次近似值。重复以上过程，得r的近似值序列，其中，xn+1=xn−f(xn)f′(xn)称为r的n+1次近似值，上式称为牛顿迭代公式。

1. 算法：解单个非线性方程的牛顿迭代法

2. 输入：x0，函数f(x)；输出：x

3. k:=0;

4. While |f(xk)|>ε1或|xk−xk−1|>ε2 do

   ​ xk+1:=xk−f(xk)f′(xk)；

   ​ k:=k+1；

   End

5. x:=xk.

牛顿法也是一种不动点迭代法，相应的公式中的函数φ(x)=x−f(x)f′(x)。

定理

​ 设x∗是方程f(x)=0的单根，且f(x)在x∗附近有连续的2阶导数，则牛顿法产生的解序列至少是局部2阶收敛的。

%%牛顿迭代法求解f(x)=x^4-x-2=0,x_0=1.5
k=0;xk=1.5;
while abs(xk^4-xk-2)>5*eps
    xk=(3*xk^4+2)/(4*xk^3-1)
    k=k+1;
end
x=xk;

k =5

x = 1.353209964199325

判停准则

迭代过程的判停准则一般有两个：

1. 残差判据，即要求|f(xk)|≤ε1,其中ε1为某个阈值
2. 误差判据，即要求|xk+1−xk|≤ε2，其中ε2为某个阈值。

在实际应用时往往要将这两种判据组合起来使用，有时也需要根据问题的特点和经验额外设置条件。

牛顿法的不足

1. 无法保证全局收敛性，也就是说，如果初始解x0不在局部收敛的范围内，迭代过程可能发散
2. 对函数的连续性要求较高，需要f(x)在x∗附近有连续的2阶导数
3. 每步迭代都要计算1阶导数，其计算量可能较大

2.5割线法和抛物线法

割线法

​ 平行弦法：xk+1=xk−f(xk)f′(x0)。

​ 缺点是收敛较差。

​ 割线法：割线法的基本思路是用差商来近似导数，从而避免复杂的导数计算，利用相邻两次迭代的函数值做差商，得f′(xk)≈f(xk)−f(xk−1)xk−xk−1

1. 解单个非线性方程的割线法

2. 输入：x0,x1,函数f(x)；输出：x

3. k:=1；

4. While |f(xk)|>ε1或|xk−xk−1|>ε2 do

   ​ xk+1:=xk−f(xk)f(xk)−f(xk−1)(xk−xk−1);

   ​ k=k+1;

   End

5. x:=xk.

抛物线法

实用的方程求根技术

阻尼牛顿法

当初始值x0偏离准确解x∗较远时，牛顿法可能发散。为了防止这种情况，在得到牛顿法的下一步解后，可引入一个阻尼因子缩小解的该变量。然后通过单调性要求

|f(xk+1)|<|f(xk)|,k=0,1,2,⋯

判断新的解是否可接受。设阻尼因子为λ1，则迭代新解为

xk+1=xk−λ1f(xk)f′(xk)

这个方法称为阻尼牛顿法

1. 算法：阻尼牛顿法

2. 输入：x0，函数f(x)；输出：x

3. k:=0;

4. While |f(xk)|>δ1或|xk−xk−1|>δ2 do

   ​ s:=f(xk)f′(xk);

   ​ xk+1:=xk−s;

   ​ i:=0;

   ​ While |f(xk+1)|≥|f(xk)| do

   ​ xk+1:=xk−λis;

   ​ i:=i+1;

   ​ End

   ​ k:=k+1;

   End

5. x:=xk

算法使用了一个阻尼因子序列{λi},其中每个值都在(0,1)之间，并按照递减顺序排列。当迭代解充分靠近准确解时，则不需要阻尼因子的调节。

通用求根算法zeroin

zeroin算法由Richard Brent 发表于1973年。该算法将二分法的稳定性和抛物线法、割线法的快速收敛性结合，是一种稳定、高效的通用求根算法。

zeroin算法一般用变量b表示当前迭代步的近似解，变量c为上一步的b，而变量a的作用则是与b构成有根区间。算法主要包括以下步骤：

1. 选取初始值a和b，使得f(a)和f(b)的正负号正好相反；
2. 将a的值赋给c
3. 重复下面的步骤，直到|f(b)|≤ε1或|a−b|≤ε2|b|，ε1、ε2为误差控制阈值
   
   1. 若f(b)的正负号与f(a)的相同，将c赋值给a;
   2. 若|f(a)|<|f(b)|，则将b的值赋给c，然后对调a、b的值;
   3. 如果c≠a，利用a、b、c以及他们的函数作逆二次插值法的一次迭代，否则执行割线法中的一步
   4. 如果执行一步逆二次插值法或割线法得到的近似解“比较满意”，将他赋值给b，否则执行一步二分法得到b，然后将上一步的b赋值给c.

zeroin 算法将方程的根困在不断缩小的区间中，很稳定，也兼顾了割线法、逆二次插值法收敛快的特点。它的主要优点如下：

1. 本身不要求函数f(x)具有光滑性
2. 不需要计算导数f′(xk)，只需要有办法算出任一xk对应的f(xk)
3. 初始解只是包含准确解的区间，不需要和准确解很接近
4. 算法简单、稳定，每步迭代都使有根区间缩小

附上实现zeroin算法的MATLAB代码