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

最新推荐文章于 2024-07-07 23:58:36 发布

OopspoO

最新推荐文章于 2024-07-07 23:58:36 发布

阅读量2.7k

点赞数 2

分类专栏： Read_Notes 文章标签：数值分析与算法

本文链接：https://blog.csdn.net/cclethe/article/details/54694430

版权

Read_Notes 专栏收录该内容

72 篇文章 1 订阅

订阅专栏

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$ ，则误差

$| x k - x * | < (b k - a k) / 2 = (b 0 - a 0) / 2 k + 1, k = 0, 1, 2, \dots .$ $\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]$ ，

$| φ (x 1) - φ (x 2) | \leq L | x 1 - x 2 |$ $\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^*$ ，并有误差估计：

| x k - x * | \leq L k 1 - L | x 1 - x 0 |

$\left | x_k-x^* \right |\leq \frac{L^k}{1-L}\left | x_1-x_0 \right |$
定理：对于不动点迭代法

xk+1=φ(xk) $x_{k+1}=\varphi(x_k)$ ，若在所求根

x∗ $x^*$ 的邻域上函数

φ(x) $\varphi(x)$ 的

p $p$ 阶导数连续，

p≥2 $p\geq 2$ ，则该迭代法在

x∗ $x^*$ 的邻域上

p $p$ 阶收敛的充分必要条件是：

φ′(x∗)=φ"(x∗)=⋯=φp−1(x∗)=0 $\varphi'(x^*)=\varphi"(x^*)=\cdots=\varphi^{p-1}(x^*)=0$ ，且

φp(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^*$ 较远时，牛顿法可能发散。为了防止这种情况，在得到牛顿法的下一步解后，可引入一个阻尼因子缩小解的该变量。然后通过单调性要求

$| f (x k + 1) | < | f (x k) |, k = 0, 1, 2, \dots$ $\left | f(x_{k+1}) \right |< \left | f(x_k) \right |,k=0,1,2,\cdots$
判断新的解是否可接受。设阻尼因子为 $\lambda_1$ ，则迭代新解为
$x k + 1 = x k - λ 1 f ( x k ) f ' ( x k )$ $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$
- 重复下面的步骤，直到 |f(b)|≤ε1 或 |a−b|≤ε2|b| ， ε1、ε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 算法将方程的根困在不断缩小的区间中，很稳定，也兼顾了割线法、逆二次插值法收敛快的特点。它的主要优点如下：
  1. 本身不要求函数 $f(x)$ 具有光滑性
  2. 不需要计算导数 $f'(x_k)$ ，只需要有办法算出任一 $x_k$ 对应的 $f(x_k)$
  3. 初始解只是包含准确解的区间，不需要和准确解很接近
  4. 算法简单、稳定，每步迭代都使有根区间缩小
  附上实现zeroin算法的MATLAB代码