计算方法-[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]$范围。

  2. 在判断子区间${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$
$$
\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