数值分析 非线性方程与方程组的数值解法

引言

在自然科学和工程技术中，非线性问题是经常出现的，如求解曲线与直线交点这样简单的数学问题，即便我们给出的曲线为简单的 $y=sin(x)$，直线方程为 $y=kx+b(k,b均为常数)$，且交点存在，在大多数情况下，仍然无法解析求出非线性方程$sin(x)=kx+b$的解 x . 很多实际应用的模型都是非线性问题，以往采用的方法是，简化模型形成线性问题进行求解，解的精确性常常得不到保证. 随着计算技术的进步和计算机的普及，为了得到更符合实际的解，人们已经开始研究非线性科学. 非线性问题要使用计算机进行求解，往往需要先转换成非线性方程或非线性方程组，因此对非线性方程（组）的求解问题一直是人们研究的课题之一.
一般，非线性方程分为代数方程和超越方程两种. 如果 $f(x)$ 是n(n>1)次多项式，则称方程为代数方程或 n 次多项式方程，否则称为超越方程.

概述

定理 1（根的存在性） 若 $f(x)$ 在$[a,b]$上连续，且 $f(a)f(b)<0$ ，则方程 $f(x)=0在[a,b]$上至少存在一个实根.
用数值计算方法求解方程（2.1）的根一般分成两个步骤进行.

1. 根的隔离. 即确定根所在的区间，使方程在这个区间内有且仅有一个根，所得区间称为方程根的隔离区间.
2. 根的精确化. 用一种方法将近似值精确化，使其满足精度要求.

显然，所求隔离区间越小越好. 通常，隔离区间的确定方法为

1. 作 $y=f(x)$ 的草图，由 $y=f(x)$与横轴交点的大致位置来确定；或者将$f(x)=0$ 改写成$f_1(x)=f_2(x)$，根据 $y=f_1(x)$和$y=f_2(x)$ 交点横坐标来确定根的隔离区间.
2. 逐步搜索：从 a 点出发，选取适当的步长 h ，根据定理 1，比较$f(a+ih)$与f的符号来确定根的隔离区间.

在具体运用上述方法时，步长的选择是关键. 如果步长 h 过小，且区间长度过大，则会使计算量增大；如果 h 过大，则有可能产生漏根现象. 因此，这种根的隔离法，只适用于求根的初始近似值.
根的逐步精确化的方法，包括二分法、迭代法、牛顿法和割线法等，我们将在以下几节中介绍上述方法，并重点学习迭代法的思想.

二分法

若非线性方程$f(x)=0$中的 $f(x)$ 在$[a,b]$ 上连续，且严格单调,$f(a)f(b)$，则非线性方程在[a,b] 上有且仅有一个根. 此时可以使用二分法求出该单根.
二分法的基本思想是，逐步将含根区间二等分，通过判别区间端点的函数值符号，进一步搜索含根区间，使含根区间长度缩小到充分小，从而求出满足给定精度的根的近似值.

利用二分法求解方程单根的具体步骤如下.

1. 计算$f(x)$在区间$[a,b]$端点处的值$f(a),f(b)$

2. 计算$f(x)$在区间中点$x_1=\frac{a+b}{2}$处的值$f(x_1)$

3. 判断，若$|f(x_1)|<\eta$, $\eta$为事先给定的精度，则$x_1=\frac{a+b}{2}$就是所求近似根，否则检验；
   若$f(x_1)$与$f(a)$同号，则根位于区间$[x_1,b]$,则取$a=x_1$
   若$f(x_1)$与$f(a)$异号，则根位于区间$[a,x_1]$,则取$b=x_1$

4. 若$|b-a|<\epsilon$($\epsilon$为预先给定精度)，则方程的近似根为($\frac{a+b}{2}$);否则，则会步骤2继续计算。

例题

使用用二分法求方程$f(x)=x^3-x-1=0$在区间$[1,2]$内的实根，要求精确到$10^{-2}$

解 已知$f(1)=-1<0,f(2)=5>0$，且在区间$[1,2]$内，$f^{\prime }(x)=3x^2-1>0$,则方程$f(x)=x^3-x-1=0$在区间$(1,2)$上只有一个实根。

1. 当$k=1,x_1=1.5$,计算$f(1.5)=0.875>0$,则区间在$(1,1.5)$内有根，$|b_k-a_k|>10^{-2}$,继续二分。
2. 当$k=2,x_2=1.25$,计算$f(1.25)=-0.296875<0$在区间$(1.25,1.5)$内有根，$|b_k-a_k|>10^{-2}$,继续二分
3. 若此反复二分，结果如下表。

$k$	$a_k$	$b_k$	$x_k$	$f(x_k)的符号$
0	1	2	1.5	+
1	1	1.5	1.5	-
2	1.25	1.5	1.5	+
3		1.25	1.375	-
4	1	1.3125	1.375	+
5	1	1.3125	1.3428	+
6	1	1.3125	1.3281	-
7	1	1.3203	1.3281	-

优点

二分法的优点是，算法简单，收敛性总能保证，对$f(x)$ 要求不高（只要连续即可），程序容易实现. 它的缺点是，只能求单实根，不能求重根和复根，也不能推广到方程组情形，且收敛速度仅与比值为$\frac{1}{2}$的几何级数相同，不算太快. 因此二分法常用于求根的初始近似值，然后再使用其他方法求根.

python 代码

import sympy
def Dichotomy(a, b, f, R):#a,b是区间端点，f是求根的函数，R是精确度
    x = sympy.symbols('x')
    f_a, f_b = f.subs(x, a), f.subs(x, b)
    if abs(b - a) < R:
        return (a + b) / 2
    if (f.subs(x, (a + b) / 2) > 0 and f_a > f_b) or (f.subs(x, (a + b) / 2) < 0 and f_a < f_b):
        return Dichotomy((a + b) / 2, b, f, R)
    elif (f.subs(x, (a + b) / 2) > 0 and f_a < f_b) or (f.subs(x, (a + b) / 2) < 0 and f_a > f_b):
        return Dichotomy(a, (a + b) / 2, f, R)

简单迭代法

构造原理

迭代算法是数值计算方法中一种逐次逼近的方法，其基本思想是，首先将方程 $f(x)=0$改写成某种等价的形式，由等价形式构造相应的迭代公式，然后选取方程的某个初始近似值 $x_0$代入公式，反复校正根的近似值，直到满足精度要求为止。

   将非线性方程$f(x)=0$改写成等价形式$$x=\phi (x)$$式中，$\phi (x)$连续，称为迭代函数，给定初始近似值$x_0$,构造如下迭代公式$$x_{k+1}=\phi (x_k),k=0,1,2,\cdots$$由上式得到的解的近似序列 { x k } \{x_k\} {xk​}的过程，称为简单迭代法，如果近似序列 { x k } \{x_k\} {xk​}有极限$$\underset{k\to \infty }{\lim }{x}_k=x^{\ast }$$则$x^{\ast }$为方程的根,此时称迭代法收敛；否则称迭代法发散.由于$x^{\ast }=\phi (x^{\ast })$,故$x^{\ast }$是迭代函数的不动点，因此简单迭代法又称为不动点迭代法

例题

  使用迭代法求方程$f(x)=x-10^x+2=0$在[0,1]内的根，计算结果保留 4 位有效数字.

  解 因为$f(0)=1>0,f(1)=-7<0且\forall x\in [0,1],f^{\prime }(x)<0$,所以方程在$[0,1]$之间必有一个实根，现将该方程改写为迭代方程
$$x=lg(x+2)$$
迭代格式
$$x_{k+1}=lg(x_k+2),k=0,1,2,\cdots$$
取初始值$x_0=1$,逐次算得$$x_1=0.47712,x_2=0.39395,x_3=0.37911,x_4=0.37642$$
$$x_5=0.37592,x_6=0.37583,x_7=0.37582$$

  因为$x_6$,$x_7$已趋于一致，所以$x_7$就是放我们方程的数值解。

一个方程的迭代公式并不是一致的，本题的迭代公式还有下面一种
$$x=10^x-2$$$$x_{k+1}=10^{x_k}-2,k=0,1,2,\cdots$$
  仍然取$x_0=1$算的$$x_1=10-2=8,x_2=10^8-2,x_3=10^{10^8-2}-2$$
显然，该迭代序列不趋向于某个定值，这种不收敛的迭代过程即为发散. 由例 题可以看出，迭代函数$f(x)=0$ 选得不同，相应的迭代序列收敛情况也不一样.

特点

要解决的问题

1. 如何选择初始近似解$x_0$ 及迭代函数$f(x)$才能保证迭代公式产生的序列$\phi (x)$ 收敛？
2. 当序列 { x k } \{x_k\} {xk​}收敛时，如何估计 k 次近似解的误差，即何时可以使迭代终止？
3. 若两个迭代均收敛，则如何取舍? 即如何衡量迭代收敛速度的快慢?

收敛性定理

   定理 2（收敛性定理） 设迭代函数 $\phi (x)\in C[a,b]$满足以下两个条件：

1. 当 $x\in [a,b]$ 时，有 $a\le \phi (x)\le b$；
2. 对任意 $x\in [a,b]$，有${\phi }^{\prime }(x)\le L<1$ .
      则 $x=\phi (x)在[a,b]$上有唯一根 $x^{\ast }$且对任意初值 $x_0\in [a,b]$,由迭代过程 $x_{k+1}=\phi (x_k)$ ， 得到的迭代序列$x_k$均收敛于$x^{\ast }$并且误差估计为$$|x_k-x^{\ast }|\le \frac{L}{1-L}|x_k-x_{k-1}|$$$$|x_k-x^{\ast }|\le \frac{L^k}{1-L}|x_1-x_0|$$

该定理的证明过程在此不详述

例题代码

接下来我们看一下该例题的代码

import math

phi_x = [1] #构造序列phi(x)且为将x0赋值为1

while True:
    phi_x.append(round(math.log10(phi_x[-1] + 2),5))#添加x_k+1 到序列中 
    if abs(phi_x[-1] - phi_x[-2]) < 1 / 2 * 10 ** -4:#检验精度
        break

for i in phi_x:
    print(i)

迭代结果

1
0.47712
0.39395
0.37912
0.37642
0.37592
0.37583
0.37582