计算方法-[2]非线性方程的数值解法

发表于 更新于 分类于

计算方法第二章非线性方程的数值解法课程内容整理。

以下Matlab代码于Notes/计算方法/programs/Solution of Nolinear Equations at master · fish-404/Notes 中，使用时创建一个f.m文件写入待求解函数，如：

1
2
function y=f(x)
y=sqrt(x)-cos(x);

使用牛顿法，需要另外创建一个fd.m文件，写入待求解函数一阶导函数。

二分法 (Bisection or Binary-search method)

二分法

原理

函数 \(f\) 在\([a, b]\) 上连续（\(f \in C [a, b]\))，\(f\)在区间\([a,b]\)的两个端点处异号(\(f(a)f(b) < 0)\), 方程\(f(x) = 0\)在 \((a, b)\)内至少有一个实根\(p\)，使得\(f(p) = 0\)，区间\([a,b]\)称为方程的有根区间。

采用二分法求出\(p\)的逼近过程

取\(a_0 = a,b_0 = b\)，中点\(x_0 = \frac{a_0 + b_0}{2}\)，若\(f(x_0) = 0\)，则 \(p = x_0\)，算法停止，否则如果\(f(a_0)f(x_0)<0\)，取\(a_1 = a_0\)，\(b_1 = x_0\)；如果\(f(x_0)f(b_0) <0\)，取\(a_1 = x_0\)，\(b_1 = b_0\)，方程的有根区间缩小为\([a_1,b_1]\)。

重复上述步骤，不断缩小有根区间，直至找到满足精度的解。

算法的停止条件可以是： * \(|x\_n-x\_{n-1}|<\varepsilon\) * \(\frac{|x\_n-x\_{n-1}}{|x\_ n|}<\varepsilon, \quad x\_n \not= 0\) * \(|f(x_n)|<\varepsilon\)

注意： 1. 在使用计算机进行数值逼近时，需注意计算精度问题，计算区间中心值时，应选用公式 \[ x_n = a_n + \frac{b_n-a_n}{2} \] 因为当\(a_n\)和\(b_n\)都接近机器精度时，\(x_n = \frac{a_n + b_n}{2}\)得到的值可能超出\([a_n，b_n]\)范围。

  1. 在判断子区间\({a_n,b_n}\)中是否存在\(f\)的根时，应使用符号函数： \[ sgn(x) = \begin{cases} -1, \quad x < 0 \\\\ 0, \quad x = 0 \\ 1,\quad x > 0 \end{cases} \] 检查\(f(a_n)\)与\(f(b_n)\)是否异号，使用\(sgn(f(a_n)\cdot sgn(f(b_n))<0\)，而不是\(f(a_n)\cdot f(b_n) <0\)，避免乘积越界而导致错误。

定理2.1 设\(f\in C[a,b]\)，\(f(a)f(b)<0\)，则二分法产生的数列\({x_n}\)满足 \[ |x_n - p| \le \frac{b-a}{2^{n+1}} \] 其中\(p\in[a,b]\)是\(f(x) = 0\)的根，\(n = 0,1,\cdots\)。

  • 收敛速率 由定理2.1，当\(n\to \infty\)时，数列\({x_n}\)以\(O(\frac{1}{2^n})\)的速度收敛到实际解\(p\) \[ x_n = p + O(\frac{1}{2^n}) \] 收敛速度与比值为\(\frac{1}{2}\)的等比级数相同。

  • 估计求解步数 对于给定的精度\(\varepsilon\)，二分法找到满足某个精度的解所需步数\(n+1\) 满足 \[ \frac{b-a}{2^{n+1}} < \varepsilon \Rightarrow n+1 > \frac{ln(b-a)-lnc}{ln 2} \] 即： \[ n + 1\ge \lceil \frac{ln(b-a) - ln\varepsilon}{ln2} \rceil \] 实际误差可能比这个值小得多。

优缺点

  • 优点:
    • 简单，总会收敛到解
    • 只要求函数连续即可
  • 缺点：
    • 收敛速度慢，只利用函数值的正负
    • 只能用于求实根，无法用于求复根和偶重根

二分法Matlab代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
% Bisection: get the solution of function f in [a,b]
% a : the interval lower bound 
% b : the interval upper bound
% varphi : the presicion
function Bisection(a, b, varphi)
clc
i=1;
 fprintf('n\t \t a_n \t\t b_n\t\t x_n\t |f(x_n)|\t \n')
 while i<100
    p=a+(b-a)/2;
    w=f(p);
    if w==0 || ((b-a)/2)<varphi
        fprintf('%d\t%8.5g\t%8g\t%8.5g\t%8.5g\t%8.5g\t%8.5g\n',i - 1,a,b,p,w,(b-a)/2,(b-a)/2/abs(p));
        return;
    end
    fprintf('%d\t%8.5g\t%8.5g\t%8.5g\t%8.5g\t%8.5g\t%8.5g\n',i - 1,a,b,p,w,(b-a)/2,(b-a)/2/abs(p));
    i=i+1;  
    if w*f(a)<0
        b=p;
    else a=p;
    end  
end

线性插值二分法

原理

在二分法的基础上，利用函数在区间两端点\(a\)和\(b\)上的连线与\(x\)轴的交点，缩小有根区间。 \[ x_n = b_n - f(b_n) \cdot \frac{b_n-a_n}{f(b_n)-f(a_n)} \]

在一定条件下，线性插值二分法比二分法的收敛速度快些。

线性插值二分法Matlab代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
% Linear interpolation dichotomy
% a : the interval lower bound 
% b : the interval upper bound
% varphi : the presicion

function LinearInterpolation(a, b, varphi)
clc
fprintf('One point secant Iteration\n');
fprintf('n\t x_n\n');
i=1;
fa=f(a);
fb=f(b);
while i<100
    p=b-fb*(b-a)/(fb-fa);
    if abs(p-b)<varphi
        fprintf('%d\t%g\n',i,p);
        return;
    end
    fprintf('%d\t%g\n',i,p);
    i=i+1;
    q=f(p);
    if q*fb<0
        a=b;
        fa=fb;
    end
    b=p;
    fb=q;
end

迭代法 (Iteration Method)

将区间\([a,b]\)上的方程\(f(x) = 0\)改写成等价形式： \[ x = \varphi(x), a\le x\le b \] 其中\(\varphi \in C[a,b]\)。 * 迭代数列\(\{x_n\}\) \[ x\_n = \varphi(x\_{n - 1}), n = 1, 2, \cdots \] 上式称为\(f(x) = 0\)的一个迭代格式， \(x_0\) 称为迭代初值

一般迭代法(Fixed-Point Iteration)(Functional Iteration)

如果数列\(\{x\_n\}\) , 设\(\lim\limits\_{n\to \infty} x\_n = p\)，则

\[ p = \lim\limits\_{n\to \infty} x\_{n+1} = \lim \limits\_{n\to \infty} \varphi(x\_n) =\varphi(\lim \limits\_{n\to \infty } x\_n ) = \varphi(p) \]

此时\(p\)为方程 \(f(x) = 0\) 的根。 #### 不动点理论 * 定义2.1 对于函数\(\varphi\), 若点\(p\) 使得\(\varphi(p) = p\)，称\(p\)为\(\varphi\)的不动点。 ##### 不动点存在的充分条件 * 定理2.2 若函数\(\varphi \in C[a,b]\)且对于\(\forall x\in [a, b]\)有\(\varphi(x) \in [a,b]\) ，那么\(\varphi\)在\([a,b]\)上有不动点。 ##### 区间内不动点存在且唯一的充分条件 * 定理 2.3 满足定理2.2的前提下，如果函数\(\varphi\)满足李普希茨(Lipschitz)条件 \[ |\varphi(x) - \varphi(y)| \le L|x- y|, \forall x, y \in [a,b] \] 并且常数\(0 < L <1\)，那么\(\varphi\)在\([a,b]\)上有唯一的不动点。 * 定理2.4 满足定理2.2的前提下，如果\(\varphi\)在\((a,b)\)上有一阶连续导数，且存在常数\(0<L<1\)使得 \[ |\varphi'(x) |\le L, \forall x \in[a,b] \] 那么\(\varphi\)在\([a,b]\)上有唯一的不动点。

注意 ，定理2.3与定理2.4为充分条件，也即即使\(\varphi\)不满足这两个定理，也可能存在不动点唯一的情况。

迭代收敛的条件

  • 定义2.2 满足定理2.3或定理2.4的函数\(\varphi\)称为区间\([a,b]\)上的压缩映射，式中的\(L\)称为\(\varphi\)的压缩系数。

  • 定理2.5 **迭代收敛定理* 若\(\varphi\)是区间\([a,b]\)上的压缩映射，那么对于\(\forall x\_0 \in [a,b]\)，由迭代公式\(x\_n = \varphi(x\_{n-1})\)产生的数列\(\{x\_n\}\)收敛到方程\(x = \varphi(x)\)在\([a,b]\)上的唯一解。

不动点迭代法Matlab代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% Functional Iteration / Fixed-Point Iteration
% varphi : the presicion
% x0 : the initial value
function FunctionalIteration(varphi, x0)
clc
i=1;
fprintf('n\t x_n\n');
while i<100
    x=f(x0);
    if abs(x-x0)<varphi
        fprintf('%d\t%g\n',i,x);
        return;
    end
    fprintf('%d\t%g\n',i,x);
    i=i+1;
    x0=x;
end

迭代算法理论

  • 迭代法的几何意义 相当于把曲线\(y = \varphi(x)\)与直线\(y = x\)的求交问题，近似地转化为数列\(\{x\_n\}\)向不动点逼近的过程，其中\(x\_n = \varphi(x\_{n-1}), n = 1, 2, \cdots\)。
  • 定理 2.6 压缩系数与迭代收敛速度的关系 若\(\varphi\)是区间\([a,b]\)上的压缩映射，迭代过程中解的误差满足 \[ |x\_n - p| <L^n \max(x\_0 - a, b-x\_0) \] \[ |x\_n - p| \le \frac{L}{1-L}|x\_n - x\_{n - 1}|\le \frac{L^n}{1-L}|x\_1 - x\_0| \] 其中\(p\)为区间\([a,b]\)上的不动点，\((0<L<1)\)为\(\varphi\)在区间\([a,b]\)的压缩系数，\(n = 1, 2, \cdots\)。

压缩系数\(L\)的值越小，\(x\_n\)趋近\(p\)的速度越快，而当\(L\)的值接近于1时，迭代法的收敛速度就变得很慢。

  • 定理2.7 局部收敛性 设\(p\)时函数\(\varphi\)的不动点，且\(\varphi'(p) \ne 1\)，如果存在常数\(\delta<0\)使得\(\varphi\)在区间\([p-\delta, p+\delta]\)上满足定理2.4， 那么一般迭代法得到的数列对任意初值\(x\_0 \in [p-\delta, p+ \delta]\)收敛。

收敛阶

  • 定义2.3 对于一个收敛到\(p\)的数列\({x\_n}\)， \(x\_n \ne p\)， 如果存在实数\(\lambda\)和\(\alpha\)使得 \[ \lim \limits\_{n\to \infty} \frac{|x\_{n+1} - p|}{|x\_n - p|^\alpha} = \lambda \]

  • 当\(\lambda \ne 0\)时数列\({x_n}\)的收敛阶是\(\alpha\)，\(\lambda\)为渐近误差常数。

    1. 若\(\alpha = 1\), 称数列线性收敛或1阶收敛。
    2. 若\(\alpha = 2\)，称数列平方收敛或2阶收敛。
    3. 若\(\alpha >1\)，称数列超线性收敛或\(\alpha\)阶收敛。

  • 当\(\lambda = 0\)时称数列\(\{x_n\}\)的收敛阶高于\(\alpha\)。

  • 迭代法的普遍性定理 满足压缩映射的函数\(\varphi\)得到的数列\(\{x_n\}\):

    • 当\(\varphi'(p) \ne 0\)时，对\(\forall x\_0 \in [a,b]\)，数列\(\{x\_n\}\)都收敛到区间\([a,b]\)上唯一的不动点\(p\)。
    • 当\(\varphi'(p) = 0\)，\(\varphi\)在\((a,b)\)有二阶连续导数，且\(\varphi''(p) \ne 0\)时数列的收敛阶是2.

加速收敛迭代法

埃特金方法 （Aitken)

对于一个收敛到 \(p\) 的数列\(\{x\_n\}\)，如果其收敛阶为1， 那么 \[ \lim \limits\_{n \to \infty}\frac{|x\_{n+1} - p|}{|x\_n - p|} = C， C 为常数，0< C < 1 \] 可以认为当\(n \to \infty\)时有 \[ \frac{x\_n -p}{x\_{n - 1} - p } \approx \frac{x\_{n - 1} - p}{x\_{n - 2} - p} \] 得到Aitken数列\(\hat{\{x\_n\}}\) \[ \hat{x\_n} = x\_n - \frac{(x\_n - x\_{n - 1})^2}{x\_n - 2 x\_{n - 1} + x\_{n - 2}} \]

  • 定理2.9 若数列\(\{x\_ n\}\)以1阶收敛到\(p\)，记 \[ e\_n = x\_n - p \ne 0 \] 对一切\(n\ge 0\)成立，且 \[ \lim \limits\_{n\to \infty} \frac{e\_{n+1}}{e\_n} = \lambda, \quad |\lambda| < 1 \] 则Aitken数列\(\hat\{x\_n\}\)是完全确定的，且 \[ \lim \limits\_{n\to \infty} \frac{\hat{x\_n} -p}{x\_n - p} = 0 \] 即Aitken数列\(\{\hat{x\_n}\}\)比数列\(\{x\_n\}\)更快收敛到\(p\)。 ##### Aitken 加速迭代法Matlab代码 
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    % AitkenIteration
    % x0 : the initial value
    % varphi : the presicion
    function AitkenIteration(x0, varphi)
    clc
    i=1;
    x1=f(x0);
    p=x0;
    fprintf('Aitken Iteration\n');
    fprintf('n\t x_n\n');
    fprintf('0\t%g\n',x0);
    while i<100   
        x2=f(x1);
        x=x0-(x1-x0).^2./(x2-2*x1+x0);
        if abs(x-p)<varphi
            fprintf('%d\t%g\n',i,x);
            break;
        end
        fprintf('%d\t%g\n',i,x);
        i=i+1;
        x0=x1;
        x1=x2;
        p=x;
    end

    #### Steffensen 方法 Aitken 数列可以对任意线性收敛的数列进行处理。若将Aitken方法与不动点结合则可以得到Steffensen方法
  • 定理2.10 设\(p\)是函数\(\varphi\)的不动点，那么具有下列迭代格式的函数\(\psi\)与\(\varphi\)具有相同的不动点。 \[ x\_n = \psi(x\_{n-1}) = \frac{x\_{n-1}\varphi(\varphi(x\_{n-1})) - [\varphi(x\_{n-1})]^2}{\varphi(\varphi(x\_{n-1})) - 2\varphi(x\_{n - 1})+ x\_{n-1}} =x\_{n-1} - \frac{[\varphi(x\_{n-1} - x\_{n-1}]^2}{\varphi(\varphi(x\_{n-1})) - 2\varphi(x\_{n-1}) + x\_{n-1}} \]
  • 定理2.11 设\(p\)是函数\(\varphi\)的不动点，且\(\varphi'(p)\ne1\)，如果存在常数\(\delta>0\)使得\(\varphi\)在区间\([p - \delta, p + \delta]\)上有三阶连续导数，那么 Steffensen 方法得到的数列对任意初值\(x\_0 \in [p - \delta, p+ \delta]\)二阶收敛。
Steffensen 迭代的特点
  • 对于一阶收敛到\(p\)的数列\(x\_{n+1} = \varphi(x_n)\)，在\(\varphi '(p)\ne 1\)时，Steffensen 方法得到的数列二阶收敛
  • 对于原本不收敛的数列，Steffensen 方法也可能把它改进为二阶收敛
  • 对于超线性收敛（大于1阶）数列，改用 Steffensen 方法的意义不大
  • Steffensen 方法多用于改进一阶方法
Steffensen 迭代法Matlab代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
% SteffensenIteration
% x0 : the initial value
% varphi : the presicion
function SteffensenIteration(x0, varphi)
i=1;
fprintf('Steffensen Iteration\n');
fprintf('n\t x_n\n');
fprintf('0\t%g\n',x0);
while i<100   
    x1=f(x0);
    x2=f(x1);
    x=x0-(x1-x0).^2./(x2-2*x1+x0);
    if abs(x-x0)<varphi
        fprintf('%d\t%g\n',i,x);
        return;
    end
    fprintf('%d\t%g\n',i,x);
    i=i+1;
    x0=x;
end

牛顿法/牛顿切线法/牛顿迭代法(Newton-Raphson Methon)

将非线性方程线性化，有 \[ f(x) = f(x\_n) + f'(x\_n)(x-x\_0) \] 当\(f'(x\_n)\ne0\)，用\(x\_{n+1}\)代替\(x\)，由\(f(x) = 0\)得到迭代公式： \[ x\_{n+1} = x\_n - \frac{f(x\_n)}{f'(x\_n)} \] 即 \[ \varphi(x) = x - \frac{f(x)}{f'(x)} \] 从几何角度来解释，即过\((x\_n, f(x\_n))\)作曲线\(y = f(x)\)的切线，切线与\(x\)轴交点的横坐标即为\(x\_{n+1}\)。

  • 定理2.12 牛顿法收敛的充分条件 设\(f\in C^2[a, b]\)，即\(f\)在\([a,b]\)上有二阶连续导数，若
  1. \(f(a)f(b)<0\)
  2. \(f'\)与\(f''\)在\([a,b]\)上符号保持不变
  3. \(f(x\_0)f''(x)>0, x\_0, x\in [a,b]\)（\(x\_0\)为迭代初值） 那么由牛顿迭代格式生成的数列\(\{x_n\}\)收敛于方程\(f(x) = 0\)在\({a,b]\)上的唯一解，且收敛阶为2。
  • 推论2.1 设\(f\in C^2[a,b]\)且\(f(a)f(b)<0\)，并在\([a,b]\)上\(f'f''>0\)，那么当迭代初值\(x_0 = b\)时，由牛顿迭代格式生成的数列\(\{x_n\}\)是一个严格递减有下界的数列，它收敛于方程\(f(x) = 0\)在\([a,b]\)上的唯一解，且收敛阶为2.
  • 推论2.2 设\(f\in C^2[a,b]\)且\(f(a)f(b)<0\)，并在\([a,b]\)上\(f'f''<0\)，那么当迭代初值\(x_0 = a\)时，由牛顿迭代格式生成的数列\(\{x_n\}\)是一个严格递增有上届的数列，它收敛于方程\(f(x) = 0\)在\([a,b]\)上的唯一解，且收敛阶为2.
  • 定理2.13 牛顿法的局部收敛性 设\(f\in C^2[a,b]\)，如果存在\(p \in [a,b]\)使得\(f(p) = 0\)和\(f'(p)\ne 0\)，那么存在一个正数\(\delta >0\)使得牛顿法生成的数列\(\{x_n\}\)对\(\forall x_0 \in [p - \delta, p+ \delta]\)都收敛于\(p\)，即数列\(\{x_n\}\)局部收敛于\(p\)，并且二阶收敛。

牛顿法的致命弱点在于每次迭代除了计算函数值外还要计算导数值。

牛顿法Matlab代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
% Newton-Raphson Method
% x0 : the initial value
% varphi : the presicion
function NewtonRaphsonMethod(x0, varphi)
clc
i=1;
fprintf('Newton Iteration\n');
fprintf('n\t x_n\n');
fprintf('0\t%g\n',x0);
while i<100   
    x=x0-f(x0)/fd(x0);
    if abs(x-x0)<varphi
        fprintf('%d\t%g\n',i,x);
        return;
    end
    fprintf('%d\t%g\n',i,x);
    i=i+1;
    x0=x;
end

弦截法/割线法/双点弦截法(Secant method)

使用差商代替牛顿法中的\(f'(x\_n)\)： \[ x\_{n+ 1} = x\_n - \frac{f(x\_0)(x\_0 - x\_{n-1})}{f(x\_n)-f(x\_{n-1}）} \] 几何意义： 用过点\((x\_{n-1}, f(x\_{n-1}))\)和\((x\_n, f(x\_n))\)的割线与\(x\)轴的交点逼近方程\(f(x) = 0\)的解。 * 定理2.14 设\(f(p) = 0\)，在\(p\)的一个邻域\(\Delta = [p - \delta, p+ \delta]\)内\(f \in C^2(\Delta)\)， \(f'(\Delta)\ne 0\), \(M \delta<1\)，其中 \[ M = \frac{\max \limits\_{x\in \Delta}|f''(x)|}{2 \min \limits\_{x\in \Delta}|f'(x)|} \] 则当\(x\_0x_1\in \Delta\)时，由弦截法迭代格式产生的数列\(\{x_n\} \subset \Delta\)收敛于\(p\)，且收敛阶为 \(\frac{1+\sqrt{5}}{2} \approx 1.618\) ### 弦截法Matlab代码 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
% Secant Method: get the solution of function f in [a,b]
% x0 : the initial value
% x1 : the initial value
% varphi : the presicion

function SecantMethod(x0, x1, varphi)
clc
i=1;
fprintf('Secant Iteration\n');
fprintf('n\t x_n\n');
fprintf('0\t%.5f\n',x0);
fprintf('1\t%.5f\n',x1);
w0=f(x0);
while i<100
    w1=f(x1);
    x=x1-w1*(x1-x0)/(w1-w0);
    if abs(x-x1)<varphi
        fprintf('%d\t%g\n',i+1,x);
        return;
    end
    fprintf('%d\t%g\n',i+1,x);
    i=i+1;
    x0=x1;
    x1=x;
    w0=w1;
end
## 抛物线法/Müller法（Parabola Iteration) 弦截法是通过\((x\_{n-1}, f(x\_{n-1}))\)和\((x\_n, f(x\_n))\)两点做直线与\(x\)轴的焦点得出\(x\_{n+1}\)，而抛物线法则是过\((x\_{n-2}, f(x\_{n-2})), (x\_{n-1}, f(x\_{n-1}))\)和\((x\_n, f(x\_n))\)三点作抛物线与\(x\)轴的交点作为\(x\_{n+1}\)。 迭代公式为 \[ x\_{n+1} = x\_n - \frac{2c}{b + sgn(b)\sqrt{b^2 - 4ac}} \] 抛物线法可用于求解非线性方程的实根和复根，特别适用于多项式的求根问题。收敛阶为1.839。

抛物线法Matlab代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
% Parabolalteration Method
% x0 : the initial value
% x1 : the initial value
% x2 : the initial value
% varphi : the presicionclc
function ParabolaIteration(x0, x1, x2, varphi)
i=3;
fprintf('Parabola Iteration\n');
fprintf('n\t x_n\n');
fprintf('0\t%.5f\n',x0);
fprintf('1\t%.5f\n',x1);
fprintf('2\t%.5f\n',x2);
h1=x1-x0;
h2=x2-x1;
w1=(f(x1)-f(x0))/h1;
w2=(f(x2)-f(x1))/h2;
d=(w2-w1)/(h2+h1);
while i<100
    b=w2+h2*d;
    D=sqrt(b*b-4*f(x2)*d);
    if abs(b-D)<abs(b+D)
        E=b+D;
    else
        E=b-D;
    end
    h=-2*f(x2)/E;
    p=x2+h;
    if abs(h)<varphi
        fprintf('%d\t%g\n',i, p);
        return;
    end
    fprintf('%d\t%g\n',i,p);
    i=i+1;
    x0=x1;
    x1=x2;
    x2=p;
    h1=x1-x0;
    h2=x2-x1;
    w1=(f(x1)-f(x0))/h1;
    w2=(f(x2)-f(x1))/h2;
    d=(w2-w1)/(h2+h1);
end