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

chapter2

非线性方程求根

2.1引言

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

定义：对光滑函数$f$，若$f(x^*)={f}'(x^*)=\cdots=f^{m-1}(x^*)=0$，但$f^m(x^*)\neq0$，则称$x^*$为方程的$m$重根。

2.2二分法

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

算法：二分法
输入：a，b，函数$f(x)$；输出：x。
While (b-a)>$\varepsilon $ do

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

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

a:=x;

Else

b:=x;

End

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

假设二分法得到的有根区间序列为$\{ (a_k,b_k),k=0,1,\cdots\}$，若取解$x_k=(a_k+b_k)/2$，则误差

\[\left | x_k-x^* \right | < (b_k-a_k)/2=(b_0-a_0)/2^{k+1},k=0,1,2,\cdots. \]

二分法求解方程$f(x)=x^2-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次便结束了，最终区间的两个端点已经是两个相邻的浮点数。

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

2.3不动点迭代法

基本原理

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

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

算法：基于函数$\varphi(x)$的不动点迭代法
输入：$x_0$，函数$f(x)$，$\varphi(x)$;输出：x
k:=0;
While $\left | f(x_k) \right |>\varepsilon_1$ 或 $\left | x_k-x_{k-1} \right |>\varepsilon_2$ do

$x_{k+1}$:=$\varphi(x_k)$;

k:=k+1;

End
$x=:x_k$.

全局收敛的充分条件

设$\varphi(x)\in C[a,b]$，若满足如下两个条件：

对任意$x\in[a,b]$，有$a\leq\varphi(x)\leq b$，
存在正常数$L\in (0,1)$，使对任意$x_1,x_2\in[a,b]$，

\[\left | \varphi(x_1)-\varphi(x_2) \right |\leq L\left | x_1-x_2 \right | \]

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

定理：设$\varphi(x)\in C[a,b]$满足以上两个条件，则对于任意初值$x_0\in [a,b]$，由不动点迭代法得到的序列$\{x_k\}$收敛到$\varphi(x)$的不动点$x^*$，并有误差估计：

\[\left | x_k-x^* \right |\leq \frac{L^k}{1-L}\left | x_1-x_0 \right | \]

定理：对于不动点迭代法$x_{k+1}=\varphi(x_k)$，若在所求根$x^*$的邻域上函数$\varphi(x)$的$p$阶导数连续，$p\geq 2$，则该迭代法在$x^*$的邻域上$p$阶收敛的充分必要条件是：$\varphi'(x^*)=\varphi"(x^*)=\cdots=\varphi^{p-1}(x^*)=0$，且$\varphi^p(x^*)\neq 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$的根，选取$x_0$作为$r$的初始近似值，过点$(x_0,f(x_0))$做曲线$y=f(x)$的切线$L$，$L$的方程为$y=f(x_0)+f'(x_0)(x-x_0)$，求出$L$与$x$轴的交点横坐标$x_1=x_0-\frac{f(x_0)}{f'(x_0)}$，称$x_1$为$r$的一次近似值。过点$(x_1,f(x_1))$做曲线$y=f(x)$的切线，并求该切线与$x$轴交点的横坐标$x_2=x_1-\frac{f(x_1)}{f'(x_1)}$，称$x_2$为$r$的二次近似值。重复以上过程，得$r$的近似值序列，其中，$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$称为$r$的$n+1$次近似值，上式称为牛顿迭代公式。

算法：解单个非线性方程的牛顿迭代法
输入：$x_0$，函数$f(x)$；输出：$x$
k:=0;
While $\left | f(x_k) \right |>\varepsilon_1$或$\left | x_k-x_{k-1} \right |>\varepsilon_2$ do

$x_{k+1}:=x_k-\frac{f(x_k)}{f'(x_k)}$；

$k:=k+1$；

End
$x:=x_k$.

牛顿法也是一种不动点迭代法，相应的公式中的函数$\varphi(x)=x-\frac{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

判停准则

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

残差判据，即要求$\left | f(x_k) \right |\leq \varepsilon_1$,其中$\varepsilon_1$为某个阈值
误差判据，即要求$\left | x_{k+1}-x_k \right |\leq \varepsilon_2$，其中$\varepsilon_2$为某个阈值。

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

牛顿法的不足

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

2.5割线法和抛物线法

割线法

平行弦法：$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_0)}$。

缺点是收敛较差。

割线法：割线法的基本思路是用差商来近似导数，从而避免复杂的导数计算，利用相邻两次迭代的函数值做差商，得$f'(x_k)\approx \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}$

解单个非线性方程的割线法
输入：$x_0$,$x_1$,函数$f(x)$；输出：x
$k:=1$；
While $\left | f(x_k) \right |>\varepsilon_1$或$\left | x_k-x_{k-1} \right |>\varepsilon_2$ do

$x_{k+1}:=x_k-\frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1})$;

k=k+1;

End
$x:=x_k$.

抛物线法

实用的方程求根技术

阻尼牛顿法

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

\[\left | f(x_{k+1}) \right |< \left | f(x_k) \right |,k=0,1,2,\cdots \]
判断新的解是否可接受。设阻尼因子为$\lambda_1$，则迭代新解为

\[x_{k+1}=x_k-\lambda_1 \frac{f(x_k)}{f'(x_k)} \]
这个方法称为阻尼牛顿法

算法：阻尼牛顿法
输入：$x_0$，函数$f(x)$；输出：$x$
$k:=0$;
While $\left | f(x_k) \right | >\delta _1$或$\left | x_k-x_{k-1} \right |> \delta _2$ do

$s:=\frac{f(x_k)}{f'(x_k)}$;

$x_{k+1}:=x_k-s$;

$i:=0$;

While $\left| f(x_{k+1}) \right |\geq \left | f(x_k) \right |$ do

$x_{k+1}:=x_k-\lambda_i s$;

$i:=i+1$;

End

$k:=k+1$;

End
$x:=x_k$

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

通用求根算法zeroin

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

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

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

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

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

附上实现zeroin算法的MATLAB代码

posted @ 2017-01-23 16:35 main_c 阅读(1078) 评论(0) 收藏举报

刷新页面返回顶部

CZW

be better