非线性方程的解

非线性方程的解

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\)的差值的绝对值小于等于一个数时,认为是收敛了
  • 迭代公式可能收敛,可能发散
  • 迭代公式的设计方法可以有多种,但存在某种迭代公式是发散的,即无法求得解,因此迭代公式的设计方法要合适

收敛:
image-202212040713216

发散:
image-202212040712457

例:求\(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)\)如图所示:

image-20221204084405214

根据初值\(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))\)作切线
image-20221204092412257

切线方程:

\[ 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)\)确定)
有以下四种情况:

image-202212041004321

  • 可以看出,用牛顿迭代法求得的序列\({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]\)上存在二阶导数,且满足下列条件:

  1. \(f(a)f(b) < 0\)
  2. \(f^{\prime}(x_0)\)\([a,b]\)上不为零
  3. \(f^{\prime\prime}(x_0)\)\([a,b]\)上不变号
  4. \(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\)作为牛顿迭代法初值

使用牛顿迭代法步骤

  1. 选初值\(x_0\)(为保证迭代收敛,可用收敛定理或式\((3.2)\)选定初值),并计算\(f(x_0), f^{\prime}(x_0)\)
  2. 迭代,\(x_1 = x_0 - \cfrac{f(x_0)}{f^\prime(x_0)}\),计算\(f(x_1), f^{\prime}(x_1)\)
  3. \(|x_1 - x_0| < \varepsilon\)(或\(|f(x_1)| < \varepsilon\)),迭代终止,\(x_1\)即为根,否则到4.
  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_xf_x = exp(2 * x[n]) + x[n] - 4;
修改f_1_xf_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\)

使用弦截法步骤

  1. 选定初始值\(x_0, x_1\),计算\(f(x_0),f(x_1)\)
  2. 按迭代公式 \(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)\)
  3. 判定若\(|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\)
posted @ 2023-08-19 20:31  Champrin  阅读(23)  评论(0编辑  收藏  举报