MATLAB和Python求解非线性方程组

解非线性方程

考虑以下形式的非线性方程：

$$\mathbf{f}\left(\mathbf{x}\right)=0$$

其中$\mathbf{f}:\mathbb{R}\to \mathbb{R}$是非线性函数。 问题是找到一个点$\mathbf{x}={\mathbf{x}}^{\ast }$使得$\mathbf{f}\left({\mathbf{x}}^{\ast }\right)=0$。然后，$\mathbf{x}={\mathbf{x}}^{\ast }$被称为函数$f$的根。 例如，如果$\mathbf{f}\left(\mathbf{x}\right)={\mathbf{e}}^{-\mathbf{x}}-\cos \mathbf{x}+0.5$，则$\mathbf{f}\left(\mathbf{x}\right)={\mathbf{e}}^{-\mathbf{x}}-\cos \mathbf{x}+0.5=0$的根位于曲线${\mathbf{e}}^{-\mathbf{x}}+0.5$与$\cos \left(\mathbf{x}\right)$ 的交点 如xia图4.1所示。

在本中，将讨论四种求解给定非线性方程的方法，即二分法，牛顿-拉夫森法，割线和迭代法，并在MATLAB和Python中实现。

二等分法

假设函数$\mathbf{f}\left(\mathbf{x}\right)$ 在间隔$\left[\mathbf{a},\mathbf{b}\right]$ 中是连续的，并且随着函数从a移至b，其符号从正变为负，反之亦然。 中间值定理的直接结果是，连续函数无法在不经过零的情况下在一个间隔内改变其符号。

二等分方法是一种迭代方法，该方法在每次迭代中将包含根的间隔分为两个相等的子间隔，将不包含根的一半除去，然后在另一半中寻找根。

如果间隔$\left[\mathbf{a},\mathbf{b}\right]$ 分为两个相等的子间隔$\left[\mathbf{a},\mathbf{c}\right]$ 和$\left[\mathbf{c},\mathbf{b}\right]$ ，其中$\mathbf{c}=\left(\mathbf{a}+\mathbf{b}\right)/2$是间隔$\left[\mathbf{a},\mathbf{b}\right]$ 的中点。函数$f$的根在$\left[\mathbf{a},\mathbf{c}\right]$ 或$\left[\mathbf{c},\mathbf{b}\right]$ 中。 如果它位于$\left[\mathbf{a},\mathbf{c}\right]$ 中，则$\mathbf{f}\left(\mathbf{a}\right)\cdot \mathbf{f}\left(\mathbf{c}\right)<0$。在这种情况下，我们知道根不在$\left[\mathbf{a},\mathbf{c}\right]$ 区间内，因此我们在$\left[\mathbf{c},\mathbf{b}\right]$ 中寻找它 ，二等分法继续将包含根的间隔分为两个子间隔，这样一个子间隔包含根，而另一个不包含根。 因此，该方法仅考虑包含根的间隔并丢弃不包含。

实现二分法的MATLAB代码部分如下，

 function x = Bisection(f, a, b, Epsilon)
	 while b-a ≥ Epsilon
		 c = (a+b)/2 ;
		 if f(a)*f(c) < 0
			 b = c ;
		 elseif f(b)*f(c) < 0
			 a = c ;
		 else
			 x = c ;
		 end
	end
	x=c ;

现在，我们可以从命令提示符处调用函数二分法，

>> format long
>> Epsilon = 1e-8 ;
>> f = @(x) xˆ2 - 3 ;
>> r = Bisection(f, 1, 2, Epsilon)
r =
1.732050813734531
>> s = Bisection(f, -2, -1, Epsilon)
s =
-1

实现二分法的Python代码部分如下，

In [1]: f = lambda x: x**2 - 3
In [2]: a, b = 1., 2.
In [3]: Eps = 1e-8
In [4]: x, Iters = Bisection(f, a, b, Eps)
In [5]: print(’Approximate root is:’, x, ’\nIterations:’, Iters)
Approximate root is: 1.7320508062839508
Iterations: 25

牛顿-拉普森法

如果函数$\mathbf{f}\left(\mathbf{x}\right)$ 在点${\mathbf{x}}_0$附近是连续的，则可以用以下形式编写：

$$\mathbf{f}\left(\mathbf{x}\right)=\mathbf{f}\left({\mathbf{x}}_0\right)+{\mathbf{f}}^{\prime }\left({\mathbf{x}}_0\right)\left(\mathbf{x}-{\mathbf{x}}_0\right)+\mathcal{O}\left({\left(\mathbf{x}-{\mathbf{x}}_0\right)}^2\right)$$

现在，如果${\mathbf{x}}_1$是函数$\mathbf{f}\left(\mathbf{x}\right)$ 的根，则$\mathbf{f}\left({\mathbf{x}}_1\right)=0$。即：

$$\mathbf{f}\left({\mathbf{x}}_0\right)+{\mathbf{f}}^{\prime }\left({\mathbf{x}}_0\right)\left(\mathbf{x}-{\mathbf{x}}_0\right)\approx \mathbf{f}\left({\mathbf{x}}_1\right)=0$$

根据以上公式：

$${\mathbf{x}}_1={\mathbf{x}}_0-\frac{\mathbf{f}\left({\mathbf{x}}_{\mathbf{n}}\right)}{{\mathbf{f}}^{\prime }\left({\mathbf{x}}_0\right)}$$

假设${\mathbf{f}}^{\prime }\left({\mathbf{x}}_0\right)\ne 0$。

牛顿-拉夫森法是一种迭代方法，用于找到最接近函数$\mathbf{f}\left(\mathbf{x}\right)$ 到初始猜测值$x_0$的近似值。 从$x_0$开始，它会生成一个数字序列${\mathbf{x}}_0,{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,$其中，

$${\mathbf{x}}_{\mathbf{n}+1}={\mathbf{x}}_{\mathbf{n}}-\frac{\mathbf{f}\left({\mathbf{x}}_{\mathbf{n}}\right)}{{\mathbf{f}}^{\prime }\left({\mathbf{x}}_{\mathbf{n}}\right)}$$

这个数字序列收敛到$f$的最接近${\mathbf{x}}_0$的根。

要使用牛顿-拉夫森方法为函数$\mathbf{f}\left(\mathbf{x}\right)={\mathbf{x}}^2-3$求根，请注意${\mathbf{f}}^{\prime }\left(\mathbf{x}\right)=2\mathbf{x}$，因此，迭代方法的形式为：

$${\mathbf{x}}^{\left(\mathbf{n}+1\right)}={\mathbf{x}}^{\left(\mathbf{n}\right)}-\frac{\mathbf{f}\left({\mathbf{x}}^{\left(\mathbf{n}\right)}\right)}{{\mathbf{f}}^{\prime }\left({\mathbf{x}}^{\left(\mathbf{n}\right)}\right)}={\mathbf{x}}^{\left(\mathbf{n}\right)}-\frac{\left({\mathbf{x}}^{{\left(\mathbf{n}\right)}^2}-3\right)}{2{\mathbf{x}}^{\left(\mathbf{n}\right)}}$$

MATLAB函数NewtonRaphson.m实现了NewtonRaphson方法：

function [x, Iter] = NewtonRaphson(f, fp, x0, Epsilon)
	 Iter = 0 ;
	 x1 = x0 - f(x0)/fp(x0) ;
	 while abs(f(x1)) ≥ Epsilon
		 x0 = x1 ;
		 x1 = x0 - f(x0)/fp(x0) ;
		 Iter = Iter + 1 ;
	 end
	 x = x1 ;

从命令提示符处调用NewtonRaphson函数：

>> format long
>> f = @(x) xˆ2 - 3 ;
>> fp = @(x) 2*x ;
>> Epsilon = 1e-8 ;
>> x0 = 1 ;
>>[x, Iterations] = NewtonRaphson(f, fp, x0, Epsilon)
x =
1.732050810014728
Iterations =
3
>> [x, Iterations] = NewtonRaphson(f, fp, -x0, Epsilon)
x =
-1.732050810014728
Iterations =
3

Python函数NewtonRaphson的代码为：

一次运行上述代码${\mathbf{x}}_0=1.0$，另一次运行${\mathbf{x}}_0=-1.0$：

In [6]: x, Iters = NewtonRaphson(f, fp, x0, Eps)
In [7]: print(’Approximate root is:’, x, ’\nIterations:’, Iters)
Approximate root is: 1.7320508100147276
Iterations: 4
In [8]: x, Iters = NewtonRaphson(f, fp, -x0, Eps)
In [9]: print(’Approximate root is:’, x, ’\nIterations:’, Iters)
Approximate root is: -1.7320508100147276
Iterations: 4

割线法

割线法与牛顿-拉夫森法具有近似的形式，但是它不需要解析形式为${\mathbf{x}}_{\mathbf{n}}\left({\mathbf{f}}^{\prime }\left({\mathbf{x}}_{\mathbf{n}}\right)\right)$ 的$\mathbf{f}\left(\mathbf{x}\right)$ 的导数。 它将${\mathbf{f}}^{\prime }\left({\mathbf{x}}_{\mathbf{n}}\right)$ 替换为有限差分公式：

$${\mathbf{f}}^{\prime }\left({\mathbf{x}}_{\mathbf{n}}\right)\approx \frac{\mathbf{f}\left({\mathbf{x}}_{\mathbf{n}}\right)-\mathbf{f}\left({\mathbf{x}}_{\mathbf{n}}-1\right)}{{\mathbf{x}}_{\mathbf{n}}-{\mathbf{x}}_{\mathbf{n}-1}}$$

因此，割线方法的形式为：

$${\mathbf{x}}_{\mathbf{n}+1}={\mathbf{x}}_{\mathbf{n}}-\frac{{\mathbf{x}}_{\mathbf{n}}-{\mathbf{x}}_{\mathbf{n}-1}}{\mathbf{f}\left({\mathbf{x}}_{\mathbf{n}}\right)-\mathbf{f}\left({\mathbf{x}}_{\mathbf{n}-1}\right)}\mathbf{f}\left({\mathbf{x}}_{\mathbf{n}}\right)$$

从某个包含$\mathbf{f}\left(\mathbf{x}\right)$ 根的间隔$\left[\mathbf{a},\mathbf{b}\right]$ 开始，割线方法迭代逼近$\left[\mathbf{a},\mathbf{b}\right]$ 中$\mathbf{f}\left(\mathbf{x}\right)$ 的零。

MATLAB函数Secant实现割线方法。 它接收一个函数$f$，其间隔$\left[\mathbf{a},\mathbf{b}\right]$ 的界限包含$f$的根，并且公差$\mathbf{\epsilon }>0$。它应用割线方法返回近似解$x$和迭代次数。

function [x, Iterations] = Secant(f, a, b, Eps)
	 x = b - ((b-a)*f(b))/(f(b)-f(a)) ;
	 Iterations = 1 ;
	 while fabs(f(x)) ≥ Eps
		 a = b ;
		 b = x ;
		 x = b - ((b-a)*f(b))/(f(b)-f(a)) ;
		 Iterations = Iterations + 1 ;
	 end

调用MATLAB函数以找到${\mathbf{x}}^2-3=0$的近似根，如下所示。

>> f = @(x) xˆ2-3 ;
>> a = 1; b = 2 ;
>> Eps = 1e-8 ;
>> [x, Iterations] = Secant(f, a, b, Eps)
x =
1.732050807565499
Iterations =
5

函数Secant的Python代码如下。

为了找到等式${\mathbf{x}}^2-3=0$的根，使用了以下Python指令：

固定点的迭代方法

点${\mathbf{x}}^{\ast }$是函数$\mathbf{g}\left(\mathbf{x}\right)$ 的不动点，如果

$$\mathbf{g}\left({\mathbf{x}}^{\ast }\right)={\mathbf{x}}^{\ast }$$

假设函数$\mathbf{f}\left(\mathbf{x}\right)$ 是可微的，可以写成

$$\mathbf{f}\left(\mathbf{x}\right)=\mathbf{g}\left(\mathbf{x}\right)-\mathbf{x}$$

如果$x_1$是$\mathbf{f}\left(\mathbf{x}\right)$ 的根，则

$$\mathbf{f}\left({\mathbf{x}}_1\right)=\mathbf{g}\left({\mathbf{x}}_1\right)-{\mathbf{x}}_1=0⟹{\mathbf{x}}_1=\mathbf{g}\left({\mathbf{x}}_1\right)$$

这意味着$x_1$是函数$\mathbf{g}\left(\mathbf{x}\right)$ 的不动点。 朝固定点进行迭代的方法的思想是将函数$\mathbf{f}\left(\mathbf{x}\right)$ 编写为$\mathbf{g}\left(\mathbf{x}\right)-\mathbf{x}$的形状，并从一些初始猜测$x_0$开始，该方法生成数字${\mathbf{x}}_0,{\mathbf{x}}_1,{\mathbf{x}}_2,\cdots$ 的序列，使用迭代规则收敛到函数$\mathbf{g}\left(\mathbf{x}\right)$ 的固定点：

$${\mathbf{x}}_{\mathbf{n}+1}=\mathbf{g}\left({\mathbf{x}}_{\mathbf{n}}\right)$$

为了说明迭代方法的工作原理，将使用它来找到函数的根

$$\mathbf{f}\left(\mathbf{x}\right)={\mathbf{x}}^2-3$$

首先，函数${\mathbf{x}}^2-3$以$\mathbf{g}\left(\mathbf{x}\right)-\mathbf{x}$的形式编写。 一种可能的选择是编写：

$$\mathbf{f}\left(\mathbf{x}\right)={\mathbf{x}}^2-3={\mathbf{x}}^2+2\mathbf{x}+1-2\mathbf{x}-4={\left(\mathbf{x}+1\right)}^2-2\mathbf{x}-4$$

如果${\mathbf{x}}^{\ast }$是$\mathbf{f}\left(\mathbf{x}\right)$ 的根，则：

$$\mathbf{f}\left({\mathbf{x}}^{\ast }\right)={\left({\mathbf{x}}^{\ast }+1\right)}^2-2{\mathbf{x}}^{\ast }-4=0$$

从中得出，

$${\mathbf{x}}^{\ast }=\frac{{\left({\mathbf{x}}^{\ast }+1\right)}^2-4}{2}$$

从初始点$x_0$开始，为$\mathbf{f}\left(\mathbf{x}\right)$ 求根的迭代方法是：

使用MATLAB和Python solve函数

使用非线性方程组

详情请参考 - 亚图跨际