非线性方程的解

From 2022-12-2 To 2022-12-
Learning from 物理学中的非线性方程的逐步搜索法和二分法求解非线性方程的迭代法弦截法

带有空气阻力的抛射体运动模型，最后解得的是非线性、超越方程，因此如何解非线性方程是所需要知道的

逐步搜索法

二分法

简单迭代法(不动点迭代法)

设$f(x)=0$，$f(x)$在有根区间$[a,b]$上是连续函数。
设计一个迭代公式，将$f(x)=0$写成等价形式:

\[ x = \varphi(x) \]

在$[a,b]$上任取$x_0$，带入上式右端，记所得的值为$x_1$

\[ x_1 = \varphi(x_0) \]

同理得$x_2 = \varphi(x_1), x_3 = \varphi(x_2)$\
于是迭代公式为:

\[ x_{n + 1} = \varphi(x_n), \, n = 0,1,2,3,... \tag{1} \]

若

\[\lim_{n \to \infty} x_{n+1} = \lim_{n \to \infty} \varphi(x_n) = x^* = \varphi(x^*) \]

称$x_n$为第$n$次近似

注：

当$x_{n + 1}与x_n$的差值的绝对值小于等于一个数时，认为是收敛了
迭代公式可能收敛,可能发散
迭代公式的设计方法可以有多种，但存在某种迭代公式是发散的，即无法求得解，因此迭代公式的设计方法要合适

收敛:

发散:

例：求$e^{2x} + x - 4 = 0$的根

将原方程改写为：

\[ e^{2x} = 4 - x \\ \Rightarrow 2x = \ln(4 - x) \\ \Rightarrow x = \cfrac{1}{2}\ln(4 - x) \]

迭代公式即为：

\[ x_{n + 1} = \cfrac{1}{2}\ln(4 - x_n) \, , n = 0,1,2,3,... \]

程序求解：

#include <iostream>
#include <cmath>
#include <vector>

#define EPS 1e-6

using namespace std;

int main(){
	cout << "process:" << endl; //迭代过程： 
	
	vector<double> x;
	x.push_back(0.0); // 等价于 x[0] = 0
	
	cout << x[0] << endl;
	int n;
	for(n = 0; ; ++n){
		x.push_back(0.5 * log(4 - x[n])); // 等价于 x[n + 1] = 0.5 * log(4 - x[n])
		cout << x[n + 1] << endl;
		if(abs(x[n + 1] - x[n]) <= EPS) // x[n + 1]与x[n]的差值小于等于EPS时，认为是收敛了
			break;
	}

	cout << "times = " << n + 1 << endl;
	cout << "result = " << x[n + 1] << endl; //输出结果 
	
	return 0;
}

结果：

process:
0
0.693147
0.597998
0.612182
0.610093
0.610401
0.610356
0.610362
0.610361
times = 8
result = 0.610361

另一种形式：

\[ x = 4 - e^{2x} \]

迭代公式即为：

\[ x_{n + 1} = 4 - e^{2x_n} \, , n = 0,1,2,3,... \]

程序求解(主要代码)：
将x.push_back(0.5 * log(4 - x[n]));改为x.push_back(4 - exp(2*x[n]));

结果：

process:
0
3
-399.429
4
-2976.96
4
-2976.96
4
-2976.96
4
-2976.96
4
-2976.96
4
-2976.96
4
-2976.96
4
...

程序陷入死循环，显然是发散，不收敛的，无法求得解

迭代函数收敛判定

若$\varphi(x)$在$x^*$的某个邻域内有一阶连续导数，且对该邻域内的$x$有$|\varphi^\prime(x)|\leq q < 1$,则由微分中值定理得：

\[ |x_{n+1} - x^*| = |\varphi^\prime(\xi)(x_n - x^*)| \leq q|x_n - x^*| \]

反复递推得：

\[ |x_n - x^*| \leq q|x_{n-1} - x^*| \leq q^2|x_{n-2} - x^*| \leq ... \leq q^n|x_0 - x^*| \]

同时

\[ \begin{aligned} &\because q < 1 \\ &\therefore 当\, n \to \infty \, , q^n \to 0\\ &即 |x_n - x^*| \to 0 \, , x_n \to x^* \end{aligned} \]

因此，递推公式若是收敛的，则需满足其在某一领域内有一阶连续导数，且导数的绝对值小于1

迭代函数收敛定理

综上，有以下定理：
设$x^*$是方程$x=\varphi(x)$的根，若$\varphi(x)$在$x^*$的某个邻域内有一阶导，且对该邻域内的一切$x$有

\[ |\varphi^\prime(x)|\leq q < 1 \]

则迭代公式$x_{n + 1} = \varphi(x_n)$对该邻域内任一初值$x_0$均收敛
迭代常用条件 $|x_n - x_{n+1}| < \varepsilon$

在上面的例子中：
第一个迭代公式：

\[ x_{n + 1} = \cfrac{1}{2}\ln(4 - x_n) \\ |\varphi^\prime(x_n)| = \cfrac{1}{2}|\cfrac{1}{4 - x_n}| < 1 \,(x^*的领域内) \]

因此是收敛的

第二个迭代公式：

\[ x_{n + 1} = 4 - e^{2x_n} \\ |\varphi^\prime(x_n)| = 2e^{2x_n} \,不满足小于1 \,(x^*的领域内) \]

因此是发散的

迭代法中有些迭代虽收敛，但速度慢!

加速算法 - 埃特金(Aitken)法

方程$x = \varphi(x)$如图所示：

根据初值$x_0$，得到$A(x_0,\bar{x_1})，其中\bar{x_1} = \varphi(x_0)$
再由$\bar{x_1}$，得到$B(\bar{x_1},\bar{x_2})，其中\bar{x_2} = \varphi(\bar{x_1})$
点$A$与$B$的连线与$y = x$得到一交点$C(x_1,x_1)$
根据$x_1$得到$(x_1,\bar{x_3})$、$(\bar{x_3},\bar{x_4})$、两者的连线与$y=x$的交点$(x_2,x_2)$
根据$x_2$ ... 得到$x_3$
根据$x_3$ ... 得到$x_4$
...

即：
第一次，计算：

\[ \bar{x_1} = \varphi(x_0) \\ \bar{x_2} = \varphi(\bar{x_1}) \]

连接$A(x_0,\bar{x_1}),B(\bar{x_1},\bar{x_2})$，与$y=x$交于点$C(x_1,x_1)$，则有$AC$斜率 $=$ $CB$斜率：

\[ \cfrac{x_1 - \bar{x_1}}{x_1 - x_0} = \cfrac{\bar{x_2} - \bar{x_1}}{\bar{x_1} - x_0} \\ \Rightarrow x_1 = \cfrac{x_0\bar{x_2} - \bar{x_1}^2}{x_0 - 2\bar{x_1} + \bar{x_2}} \]

第二次，计算：

\[ \bar{x_1} = \varphi(x_1) \\ \bar{x_2} = \varphi(\bar{x_1}) \]

同理可得：

\[ x_2 = \cfrac{x_1\bar{x_2} - \bar{x_1}^2}{x_1 - 2\bar{x_1} + \bar{x_2}} \]

第三次：

\[ \bar{x_1} = \varphi(x_2) \\ \bar{x_2} = \varphi(\bar{x_1}) \\ x_3 = \cfrac{x_2\bar{x_2} - \bar{x_1}^2}{x_2 - 2\bar{x_1} + \bar{x_2}} \]

...

综上，
第$n + 1$次，$n = 0,1,2,...$：

\[ \bar{x_1} = \varphi(x_n) \tag{2.1} \]

\[ \bar{x_2} = \varphi(\bar{x_1}) \tag{2.2} \]

\[ x_{n + 1} = \cfrac{x_n\bar{x_2} - \bar{x_1}^2}{x_n - 2\bar{x_1} + \bar{x_2}} \tag{2.3} \]

收敛定理

与函数迭代法一致

例：求$e^{2x} + x - 4 = 0$的根

将原方程改写为：

\[ x = \cfrac{1}{2}\ln(4 - x) \]

程序求解(主要代码)：

	vector<double> x;
	x.push_back(0.0); // 等价于x[0] = 0
	
	double bar_x1, bar_x2; // 对应\bar{x_1}和\bar{x_2}
	
	cout << x[0] << endl;
	int n;
	for(n = 0; ; ++n){
		bar_x1 = 0.5 * log(4 - x[n]); //计算\bar{x_1}
		bar_x2 =  0.5 * log(4 - bar_x1); //计算\bar{x_2}
		x.push_back((x[n] * bar_x2 - pow(bar_x1, 2))/(x[n] - 2 * bar_x1 + bar_x2)); //计算x[n + 1] 
		cout << x[n + 1] << endl;
		if(abs(x[n + 1] - x[n]) <= EPS) // x[n + 1]与x[n]的差值小于等于EPS时，认为是收敛了
			break;
	}

结果：

process:
0
0.609483
0.610362
0.610362
times = 3
result = 0.610362

直观的可见，相比于普通函数迭代法(见上)，普通函数迭代次数为$8$，而经过加速算法，迭代次数缩减至$3$，收敛速度要快，精度也在$1e-6$

注：对某些发散过程，埃特金法也适用

牛顿迭代法(切线法)

设方程$f(x)=0$
首先在根$x^*$的附近取一点$x_0$作为初始近似值，过曲线$y=f(x)$上的点$(x_0,f(x_0))$作切线

切线方程：

\[ y = f(x_0) + f^\prime(x_0)(x - x_0) \]

第一次：
若$f^\prime(x_0) \neq 0$，该切线与$x$轴的交点为

\[ x_1 = x_0 - \cfrac{f(x_0)}{f^\prime(x_0)} \]

$x_1$作为根$x^*$的第一次近似值

第二次：
以$x_1$过曲线$y=f(x)$上的点$(x_1,f(x_1))$作切线，若$f(x_1) \neq 0$，该切线与$x$轴交点：

\[ x_2 = x_1 - \cfrac{f(x_1)}{f^\prime(x_1)} \]

$x_2$作为根$x^*$的第二次近似值

...

综上，
第$n + 1$次，$n = 0,1,2,...$：

\[ x_{n + 1} = x_n - \cfrac{f(x_n)}{f^\prime(x_n)} \tag{3.1} \]

当近似值越趋近于方程的解$x^*$时，所得斜率会越来越大，即趋近方程的解$x^*$的迭代速度会越来越快

相比于埃特金(Aitken)法，埃特金(Aitken)法计算量大

收敛判断

牛顿法仍是迭代法，牛顿法迭代函数为：

\[ \varphi(x) = x - \cfrac{f(x)}{f^\prime(x)} \]

若$f^\prime(x)$在根$x$的某个邻域内不为零，且$f^{\prime\prime}(x)$存在，

牛顿迭代法中，要求$f^\prime(x_n) \neq 0$，因此$f^\prime(x)$在根$x$的某个邻域内不为零
判断迭代函数的收敛性需要对\varphi(x)求导，即求$|\varphi^\prime(x)|$，因此$f^{\prime\prime}(x)$存在

由迭代函数的收敛性判断定理可知，若：

\[ |\varphi^\prime(x)| = |\cfrac{f(x)f^{\prime\prime}(x)}{[f^{\prime}(x)]^2}|\leq q_1 < 1 \]

$q_1$为常数，则牛顿迭代法收敛

保证牛顿迭代法收敛：

用逐步搜索法，找到方程的有根区间$[a,b]$(尽量小)
使得在$[a,b]$上，$f^{\prime}(x),f^{\prime\prime}(x)$都不变号，即$f(x)$在$[a,b]$上保持严格单调(由$f^{\prime}(x)$确定)和凹向不变(由$f^{\prime\prime}(x)$确定)
有以下四种情况:

可以看出，用牛顿迭代法求得的序列${x_n}$均是单调地趋于$x^*$，故牛顿迭代法是收敛的。
凡满足关系式

\[ f(x_0)f^{\prime\prime}(x_0) > 0 \]

$x_0$均可作为初始值
例如图(1)，(4)取$x_0 = b$，图(2)，(3)取$x_0 = a$

牛顿迭代法收敛定理

设函数$f(x)$在$[a,b]$上存在二阶导数，且满足下列条件:

$f(a)f(b) < 0$
$f^{\prime}(x_0)$在$[a,b]$上不为零
$f^{\prime\prime}(x_0)$在$[a,b]$上不变号
对$x_0 \in [a,b]$，有$f(x_0)f^{\prime\prime}(x_0)>0$

则牛顿法迭代序列${x_n}$收敛于方程$f(x)=0$在$(a,b)$内的唯一根$x^*$

使用牛顿法，初始值$x_0$的选取很重要，若在$[a,b]$上，$f^{\prime}(x),f^{\prime\prime}(x)$的符号不易判别，如何选取初值$x_0$？

如何选取初值$x_0$

由：

\[ x_1 = x_0 - \cfrac{f(x_0)}{f^\prime(x_0)} \]

得：

\[ x_1 - x^* = (x_0 - x^*) - \cfrac{f(x_0)}{f^\prime(x_0)} \]

设第$n$次的近似值与方程的根$x^*$的误差为$e_n$：

\[ e_n = x_n - x^* \]

即$e_0 = x_0 - x^*$,$e_1 = x_1 - x^*$
那么：

\[ e_1 = e_0 - \cfrac{f(x_0)}{f^\prime(x_0)} \\ \Rightarrow \cfrac{e_1}{e_0} = 1 - \cfrac{f(x_0)}{f^\prime(x_0)e_0} \\ \Rightarrow \cfrac{e_1}{e_0} = 1 - \cfrac{f(x_0)}{f^\prime(x_0)(x_0 - x^*)} \\ \Rightarrow \cfrac{e_1}{e_0} = \cfrac{f(x_0) + f^\prime(x_0)(x^* - x_0)}{f^\prime(x_0)(x^* - x_0)} \\ \]

利用一阶及零阶泰勒公式泰勒公式：

\[ 0 = f(x^*) = f(x_0) + f^\prime(x_0)(x^* - x_0) + \cfrac{1}{2}f^{\prime\prime}(\xi)(x^* - x_0)^2 \\ 0 = f(x^*) = f(x_0) + f^\prime(\eta)(x^* - x_0) \]

那么：

\[ \cfrac{e_1}{e_0} \approx \cfrac{-\cfrac{1}{2}f^{\prime\prime}(\xi)(x^* - x_0)^2}{f^\prime(x_0)(x^* - x_0)} = -\cfrac{1}{2}\cfrac{f^{\prime\prime}(\xi)(x^* - x_0)}{f^\prime(x_0)} = -\cfrac{1}{2}\cfrac{f^{\prime\prime}(\xi)-\cfrac{f(x_0)}{f^\prime(\eta)}}{f^\prime(x_0)} = \cfrac{f^{\prime\prime}(\xi)f(x_0)}{2f^\prime(x_0)f^\prime(\eta)} \]

在$x_0$处可计算$f(x),f^{\prime}(x_0),f^{\prime\prime}(x_0)$的值，但无法计算$f^{\prime\prime}(\xi)$与$f^\prime(\eta)$
若$f^{\prime}(x),f^{\prime\prime}(x)$在$x_0$附近相对变化不大，即只要$f^{\prime\prime}(x_0) \neq 0$，则有近似公式:

\[ \cfrac{e_1}{e_0} \approx \cfrac{f^{\prime\prime}(x_0)f(x_0)}{2[f^\prime(x_0)]^2} \]

要使牛顿迭代法收敛，误差必须减少，即$|e_1|<|e_0|$

\[ [f^\prime(x_0)]^2 > |\cfrac{f^{\prime\prime}(x_0)}{2}|\cdot|f(x_0)| \tag{3.2} \]

只要$x_0$满足上式，且$f^{\prime\prime}(x_0) \neq 0$，就可将$x$作为牛顿迭代法初值

使用牛顿迭代法步骤

选初值$x_0$(为保证迭代收敛，可用收敛定理或式$(3.2)$选定初值)，并计算$f(x_0), f^{\prime}(x_0)$
迭代，$x_1 = x_0 - \cfrac{f(x_0)}{f^\prime(x_0)}$，计算$f(x_1), f^{\prime}(x_1)$
若$|x_1 - x_0| < \varepsilon$(或$|f(x_1)| < \varepsilon$)，迭代终止，$x_1$即为根，否则到4.
以$x_1, f(x_1), f^{\prime}(x_1)$替代$x_0, f(x_0), f^{\prime}(x_0)$回到2.继续，直到$|x_1 - x_0| < \varepsilon$(或$|f(x_1)| < \varepsilon$)为止

例：求方程$f(x) = x^3 - 2x^2 + 4x + 1 = 0$在$x = 0$附近的一个实根

对$f(x)$求导

\[ f^{\prime}(x) = 3x^2 - 4x + 4 \\ f^{\prime\prime}(x) = 6x - 4 \]

程序求解(主要代码)：

	vector<double> x;
	x.push_back(1.0); // 初值x[0]
	
	double x1, f_x, f_1_x; // 对应x_1, f(x_1), f^{\prime}(x_1)
	
	cout << x[0] << endl;
	int n;
	for(n = 0; ; ++n){
		f_x = pow(x[n], 3) - 2 * pow(x[n], 2) + 4 * x[n] + 1;
		if(abs(f_x) <= EPS) //f(x)的绝对值小于等于EPS时，认为是收敛了
			break;
		f_1_x = 3 *  pow(x[n], 2) - 4 * x[n] + 4;
		x1 = x[n] - f_x/f_1_x;
		x.push_back(x1); 
		cout << x[n + 1] << endl;
		if(abs(x[n + 1] - x[n]) <= EPS) // x[n + 1]与x[n]的差值的绝对值小于等于EPS时，认为是收敛了
			break;
	}

结果：

process:
1
-0.333333
-0.228758
-0.222515
-0.222495
times = 5
result = -0.222495

初值取0.0结果：

process:
0
-0.25
-0.222892
-0.222495
times = 4
result = -0.222495

例：求$e^{2x} + x - 4 = 0$在区间$[0.5, 1]$内的根的近似值

设：

\[ f(x) = e^{2x} + x - 4 \]

则：

\[ f^{\prime}(x) = 2e^{2x} + 1 > 0 \\ f^{\prime\prime}(x) = 4e^{2x} > 0 \]

且：

\[ f(1) = e^2 + 1 - 4 > 0 \]

因此：取初值为$x_0 = 1$

程序求解(主要代码)：
修改f_x：f_x = exp(2 * x[n]) + x[n] - 4;
修改f_1_x：f_1_x = 2 * exp(2 * x[n]) + 1;

结果：

process:
1
0.721826
0.620693
0.610454
0.610362
times = 5
result = 0.610362

可见，使用牛顿迭代法也比一般的函数迭代法($8$次)收敛速度要快

例：计算$f(x) = x^{41} + x^3 + 1 = 0$在$x_0 = -1$附近的实根

\[ f^{\prime}(x) = 41x^{40} + 3x^2 \\ \cfrac{1}{2}f^{\prime\prime}(x) = 820x^{39} + 3x \]

可知：

\[ f(-1) = -1 \\ f^{\prime}(-1) = 44 \\ \frac{1}{2}f^{\prime\prime}(-1) = -823 \]

\[ [f^{\prime}(-1)]^2 = 1936 > |\frac{f^{\prime\prime}(-1)}{2}|\cdot|f(-1)| = 823 \]

因此，可取$x_0 = -1$为初值

...

例：计算平方根

\[ x = \sqrt{c} \]

求$f(x) = x^2 - c = 0$的正根

\[ f^{\prime}(x) = 2x\\ f^{\prime\prime}(x) = 2 \\ x_{n + 1} = x_n - \cfrac{x_n^2 - c}{2x_n} = \cfrac{1}{2}(x_n + \cfrac{c}{x_n}) \]

初值的选取：符合定理或条件

用于求复根

牛顿法可求方程的实根，也可以求方程的复根
令$z_0 = x_0 + y_0i$，迭代按：

\[ z_{n + 1} = z_n - \cfrac{f(z_n)}{f^{\prime}(z_n)} \]

例：计算$f(z) = z^3 - z - 1 = 0$的复根

取初值$z_0 = -0.5 + 0.5i$

...

弦截法(割线法)

牛顿法计算时要用到函数的导数，很多情况下难以使用，更何况选取初值还需使用到二阶导数

弦截法(割线法)是它的一种改进，在牛顿迭代公式中，用差商代替导数，即

\[ f^{\prime}(x_n) \approx \cfrac{f(x_n) - f(x_{n-1})}{x_n - x_{n-1}} \]

代入牛顿迭代公式可得：

\[ x_{n+1} = x_n - \cfrac{f(x_n) }{f(x_n) - f(x_{n-1})}(x_n - x_{n - 1}) \, , n = 0,1,2,... \tag{4} \]

迭代格式中没有用到函数的导数，计算方便
但收敛速度较牛顿迭代法要慢，开始时要用到两个不同的根的近似值作初值

例：求方程$e^{2x} + x - 4 = 0$在区间$[0.5,1]$内根的近似值

设$f(x) = e^{2x} + x - 4$
两个初值取区间端点，$x_0 = 0.5, x_1 = 1$

使用弦截法步骤

选定初始值$x_0, x_1$，计算$f(x_0),f(x_1)$
按迭代公式 $x_{n+1} = x_n - \cfrac{f(x_n) }{f(x_n) - f(x_{n-1})}(x_n - x_{n - 1}) \, , n = 0,1,2,...$
计算$x_2$，再计算$f(x_2)$
判定若$|f(x_2)| < \varepsilon$，则迭代终止
否则，用$(x_2, f(x_2))和(x_1, f(x_1)),$分别代替$(x_1, f(x_1)),(x_0, f(x_0))$
重复步骤2，3，直至$|f(x_2)| < \varepsilon$