任意方程求根
简介
方程和函数是代数数学中最为重要的内容之一，从初中直到大学，我们都在研究着方程与函数，甚至我们将图形代数化，从而发展出了代数几何、解析几何的内容。而在方程与函数中，我们研究其性质最多的，往往就是方程的根（零点），即使是研究方程的极值点、鞍点等，我们无非也只是研究其微商的零点。
我们在初等数学中已经学过许多简单初等函数、线性方程的求解方法，在本文中，我们重点讨论任意方程，尤其是计算困难的非线性方程的求根方法。
方程
分类和介绍
方程就是指含有未知数的等式。是表示两个数学式（如两个数、函数、量、运算）之间相等关系的一种等式，使等式成立的未知数的值称为“解”或“根”。在这里，根据一些性质的不同，我们将方程分成以下几类：
- 单个方程
- 线性方程：本质是等式两边乘以任何相同的非零数，方程的本质都不受影响。通常认为只含有一次项的方程。
- 非线性方程：是因变量与自变量之间的关系不是线性的关系的方程。
- 多项式方程
- 超越方程：指含有未知量的超越式（指数、对数、三角函数、反三角函数等）的方程。换言之，超越方程中都有无法用含有未知数的多项式、分式或开方表示的式子。
- 多个方程
- 线性方程组
- 非线性方程组
方程的零点（根、解）
若有一个值或一些值能够使得方程 \(f(x)=0\) 成立，那么这个值就被成为方程的解，也常常被叫做零点和根。
若方程有且只有一个解\(x^*\)，那么我们称方程有单根\(x^*\)。
若对于方程\(f(x)=0\)，有\(f(x^*) = 0,f^{'}(x^*)=f^{''}(x^*)=\cdots=f^{(k)}(x^*)=0,f^{(k+1)}(x^*)\neq0\),那么称\(x^*\)为方程的k+1重根
PS：若方程是简单幂函数多项式组成，那么方程的解的数量应和最高此项的数值一致，因为存在虚根。
求根方法
求根的方法基本上大同小异，都是通过区间去逼近方程的根的点。
首先我们说一个定理1：对于实函数方程\(f(x)=0\)，当\(x\in(a,b)\)，且\(f(x)\)在\(x\in(a,b)\)时单调且连续，若\(f(a)\cdot f(b)<0\)，则方程在\(x\in(a,b)\)有且只有一个根。
二分法
二分法和算法中的二分搜索法非常的类似，取定有根区间\([a,b]\)进行对分，求得\(mid = \frac{a+b}{2}\)进行判断含根区间，如果\(f(a)\cdot f(mid)<0\),则令\(b=mid\)；反之若\(f(b)\cdot f(mid)<0\),则令\(a=mid\)。当\(|b_n-a_n|<\epsilon\)停止计算返回结果。
产生的截断误差为\(|e_{n-1}| = x_{n+1} - x^*\leq[b_n - a_n] = \frac{b_0 - a_0}{2^n}\)。
可以计算出最小迭代次数为\(n = \frac{lg(c_0-a_0)-lg\epsilon}{lg2}\)
代码实现：
private static double epsilon = 0.001;
// func为函数，写法如x=>x*x+2*x-1，a，b必须为有效的含根开区间
public static double Binary(Func<double, double> func, double a, double b)
{
var f1 = func.Invoke(a);
var f2 = func.Invoke(b);
if (f1 * f2 > 0)
throw new ArgumentException("此区间无根或根不唯一");
double mid = (a + b) / (double)2;
var fm = func.Invoke(mid);
if (fm == 0)
return fm;
if (f1 * fm < 0)
b = mid;
else if (f2 * fm < 0)
a = mid;
if (Math.Abs(b - a) <= epsilon)
return (a + b) / (double)2;
return Binary(func, a, b);
}
一般迭代法
迭代法是将方程求根问题转换成了函数求交点问题再转换成数列求极限的问题，这个过程中，对方程进行一个巧妙的变换之后，方程就可以在若干次迭代之后解出一个近似解。
操作方法如下：首先将方程\(f(x)=0\)改写成\(x = \varphi(x)\)的形式，这样就可以将方程的解看成是函数\(y=x\)和\(y=\varphi(x)\)的交点。给定初始值\(x_0\)后，则\(x_1 = \varphi(x_0)\)，不断重复这个过程，若\(\displaystyle \lim_{k \to \infty}x_k\)存在，则迭代便可以达到使得\(x_k\)趋于交点。
通常，我们需要保证迭代函数在指定的含根区间内的导数值\(|\varphi^{'}(x)|<1\)，否则迭代函数将会不收敛。
迭代过程如下图所示：
迭代法的收敛性证明
在这里，我们将证明迭代法求根的合理性和可行性。
前提条件：设函数\(\varphi(x), x\in[a,b]\)有连续的一阶导数，并且满足以下条件：
-
\(\forall x\in [a,b],\varphi(x)\in[a,b]\)
-
\(\exist L \in (0,1),\forall x\in[a,b],|\varphi^{'}(x)|\leq L\)
要证明和解决以下命题和问题：
-
在\(x\in[a,b],\exist x^*,\varphi(x^*) = 0\)
-
\(\forall x_0\in[a,b]\)，迭代过程\(x_{k+1} = \varphi(x_k)\)均收敛与\(x^*\)
-
求解误差分析式
现在开始证明
1：证明在区间内有且只有一个根存在：
证：在\(x \in [a,b]\)时，\(\varphi^{'}(x)\)存在，所以有\(\varphi(x)\)连续，于是可以作\(g(x) = x - \varphi(x)\)，易知\(g(x)\)连续。
又因为\(\varphi(x) \in [a,b]\)，且\(g(a)*g(b)<0\),故存在实根，使得\(x=\varphi(x)\)。
利用反证法：若在[a,b]上还有一实根\(\bar{x}\)，那么通过拉格朗日中值定理必定有：
\[x^*-\bar{x} = \varphi(x^*)-\varphi(\bar{x}) = \varphi^{'}(\xi)(x^*-\bar{x})\Longrightarrow \varphi^{'}(\xi) = 1 \]
显然与假设不符合。
2：证明这个根收敛于\(x^*\)
根据拉格朗日中值定理，容易知道
\[x^*-x_{k+1} = \varphi(x^*)-\varphi(x_k) = \varphi^{'}(\xi)(x^*-x_k) \]
又因\(|\varphi^{'}(x)|\leq L\)，故易得：
\[|x^*-x_{k+1}|\leq|L(x^*-x_k)|=|L^2(x^*-x_{k-1})=\cdots=|L^k(x^*-x_0)| \]
因为\(\displaystyle \lim_{k \to \infty}L^k=0\)，故\(\displaystyle \lim_{k \to \infty}|x^*-x_{k+1}|=0(绝对值必定非负)\)，得\(x^*=x_{k+1}(k\to\infty)\)
3:迭代法的误差式：
设某一次迭代后误差为\(\epsilon\)，则可以推出：
\[|x_{k+1}-x_{k}| = |(x^*-x_k)-(x^*-x_{k+1})\geq|x^*-x_k|-|x^*-x_{k+1}|\geq|x^*-x_k|-L|x^*-x_k| \]
其中从\(|x^*-x_{k+1}| \leq L|x^*-x_k|\)则是因为每一次迭代后，误差总是减少的。
故可以轻松的计算出误差估计式为：
\[|x^*-x_k|\leq\frac{1}{1-L}|x_{k+1}-x_k| \]
通过中值定理，又可以推出：
\[|x_{k+1}-x_k| = |\varphi^{'}(\xi)(x_k-x_{k-1})| \]
因为\(|\varphi^{'}(\xi)|\leq L\)，将上式递推后放缩成第二种误差估计式：
\[|x^*-x_k|\leq\frac{L^k}{1-L}|x_1-x_0| \]
// func是迭代函数而不是实际函数
public static double Iterative(Func<double, double> func, double x)
{
double temp = func.Invoke(x);
if (Math.Abs(temp - x) <= epsilon)
{
return temp;
}
x = temp;
return Iterative(func, x);
}
牛顿法
牛爵爷在整个微积分、数值分析中都有着举足轻重的地位，这里阐述的牛顿法，就是Taylor展开式的一部分。牛顿迭代法的核心思想就是：设法将一个非线性方程\(f(x)=0\)转化为某种线性方程去求解。在数学意义上，我们知道泰勒公式可以将任意函数以简单幂函数的形式表示出来，而在几何意义上，我们是利用切线与X轴的交点去进行迭代，在根处，切线必过零点。若设\(f(x)=0\)的近似解为\(x_k\)，则方程可以用一阶Taylor展开进行近似表示，牛顿迭代法的核心公式如下：
\[p_1(x) = f(x_k)+f^{'}(x_k)(x-x_k) \]
普通牛顿法
从上述公式中我们知道，其中\(p_1(x)\)就是泰勒多项式的表达，将其看成是方程\(f(x)=0\)，通过迭代的思想，从而得到了一个线性方程：
\[x_{k+1} = x_k - \frac{f(x_k)}{f^{'}(x_k)} \]
这就是普通牛顿法的迭代公式。在几何意义上，他就是一个切线与x轴的交点去逼近零点，如图：
牛顿下山法
所有的迭代法都有一个无法逃过的缺点，如果选取的初始值离根太远，则会导致迭代过程次数过多，并且有可能导致迭代发散，因为迭代法只具有局部收敛性。
为了避免迭代失败或时间过长，我们加上这样一个条件用于保证数值稳定下降
\[|f(x_{k+1})|<|f(x_k)| \]
将这个条件和牛顿法结合再一起，即再保证数值稳定下降的过程中，利用牛顿法加快迭代，这就是牛顿下山法。具体的操作如下：
将牛顿法的结果\(\bar{x}_{k+1}\)与前一步的迭代值\(x_k\)进行加权平均作为新的迭代值\(x_{k+1}\)，
则有：
\[x_{k+1} = \lambda\bar{x}_{k+1} + (1-\lambda)x_k \]
或
\[x_{k+1} = x_k - \lambda\frac{f(x_k)}{f^{'}(x_k)} \]
其中\(\lambda(0\leq\lambda<1)\)称为下山因子，它的值是一个逐步试探的过程，可以从1开始取值，一旦满足\(|f(x_{k+1})|<|f(x_k)|\)则称为下山成功，否则需要另选初始值\(x_0\)进行试算。
简单牛顿法
牛顿法可以说是一个非常有效的计算方法，准确度和迭代次数上都要比普通的迭代法要好得多，但是牛顿法最大的问题是我们需要求方程的导数，对于某些极其复杂的函数而言，导数是无法通过人工的方式计算，假如我们使用微积分的方式去求解导数，这对整个程序的性能会有比较大的影响，因此我们可以利用一个常数值\(\lambda\)来代替方程中的导数项。此时迭代公式为：
\[x_{k+1} =x_k - \frac{f(x_k)}{\lambda} \]
不过对于常数\(\lambda\)的取值是有限制的，因为我们需要保证迭代函数的收敛性，如果函数不收敛于\(x^*\)，那么一切都没有意义。于是有：
\[\varphi(x) = x - \frac{f(x)}{\lambda}\Longrightarrow\varphi^{'}(x) = 1-\frac{f^{'}(x)}{\lambda} \]
牛顿迭代法的收敛性遵循前文中普通迭代法的收敛性，于是可以得到：
\[|\varphi^{'}(x)| = |1-\frac{f^{'}(x)}{\lambda}|\Longrightarrow0<\frac{f^{'}(x)}{c}<2 \]
这样我们就可以很轻松的确定下c的值了。
简单牛顿法的几何意义就简单许多了，和我们之前讨论的普通迭代法一致，只不过普通迭代法是将函数值和\(y=x\)进行处理变换，而简单牛顿法则是和\(y= \lambda(x-x_k)+f(x_k)\)进行变换，本质是一致的。
///这里只写普通牛顿法，另外的函数由读者自己思考
// 其中f1x为导数
public static double Newtown(Func<double, double> fx, Func<double, double> f1x, double x)
{
var temp = f1x.Invoke(x);
if (temp == 0)
{
throw new ArgumentException();
}
x = x - fx.Invoke(x) / temp;
if (Math.Abs(fx.Invoke(x)) <= epsilon)
{
return x;
}
return Newtown(fx, f1x, x);
}
弦截法
弦截法是牛顿法的一种改进，保留了牛顿法中收敛速度快的优点，还克服了再牛顿法中需要计算函数导数值\(f^{'}\)的缺点。
弦截法中，我们用差商
\[\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}} \]
去代替牛顿法中的导数值，于是可以得到以下离散化的迭代
\[x_{k+1} = x_k - \frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1}) \]
这种方法叫做双点弦截法，如图所示：
从图中可以知道，弦截法一直利用两点之间的连线作为迭代的内容，那么，他的合理性在哪里呢？
整个迭代法都离不开中值定理，这里也是这样，事实上，这个差商之所以可以对导数值进行替代，是因为中值定理中说过，连续函数中两函数上的点的连线的斜率必为两点之间某一点的导数值，并且迭代过程中，这两点的中值会越来越逼近函数零点，事实上这已经说明了弦截法是牛顿法的改进方法了。
不过，如果将函数中\(x_{k-1}\)用一个定点\(x_0\)代替，这种方法叫做单点弦截法，几何意义如图所示：
//这里就对双点弦截法进行编程
public static double StringCut(Func<double, double> func, double x1, double x2)
{
var temp = x1 - (func.Invoke(x1) / (func.Invoke(x1) - func.Invoke(x2))) * (x1 - x2);
x2 = x1;
x1 = temp;
if (Math.Abs(func.Invoke(x1)) <= epsilon)
{
return x1;
}
return StringCut(func, x1, x2);
}
收敛迭代的加速
除了普通迭代法意外，我们介绍的牛顿法、弦截法都是针对迭代法的改进，那么对于普通迭代法应该如何去进行加速呢？首先我们应该对收敛速度进行一个研究。
收敛速度的计算
对于收敛的速度，我们很难直接从数值看出速度，因为收敛过程中，误差的变化是越来越小的，因此我们再次进行一个比阶。
\[\exists p\in N \And C>0,\displaystyle \lim_{k \to \infty}\frac{e_{k+1}}{e_k^p} = C \]
其中：
\[e_k = |x_k-x^*| \]
\[p= \begin{cases} p=1,线性收敛\\ p=2,平方收敛\\ p>2,高阶收敛 \end{cases} \]
并且有一个定理：若\(\varphi(x)\)在所求的根\(x^*\)邻域有连续的二阶导数，并且有\(|\varphi^{'}(x)|<1\)，有以下特点：
-
当\(\varphi^{'}(x)\neq0\)，那么迭代过程线性收敛
-
当\(\varphi^{'}(x)=0,\varphi^{''}(x)\neq0\)，那么迭代是平方收敛的。
证明过程留给读者，只需要利用泰勒展开是便可以证明该定理。
普通迭代加速
对于迭代过程，利用中值定理有
\[x^*-x_{k+1} = \varphi(x^*)-\varphi(x_k) = \varphi^{'}(\xi)(x^*-x_k) \]
当\(\Delta x \rightarrow0\)，我们将\(\varphi^{'}(\xi)\)看成定值\(\lambda(\lambda<1)\)
于是有
\[\lambda(x^*-x_k) = x^*-x_{k+1}\longrightarrow x^* = \frac{1}{1-\lambda}\bar{x}_{k+1}-\frac{\lambda}{1-\lambda}x_k \]
故可以推出最终的迭代公式为
\[\begin{cases} \bar{x}_{k+1} = \varphi(x_k)\\ x_{k+1} = \bar{x}_{k+1} + \frac{\lambda}{1-\lambda}(\bar{x}_{k+1}-x_k) \end{cases} \]
利用这种方式的好处就是再求得\(x_k\)时，利用加权的方式，使得迭代法变得有点类似像牛顿迭代法一样变成切线性质的线而不是x=c或y=c的线，你经过画图和简单的代数运算后，你会发现\(\bar{x}_{k+1}<x_{k+1}\)，也就是说达到了我们的加速目的。
Atiken加速法
Atiken加速法其实就是再普通加速迭代公式上在进行一次迭代，这里直接写出公式：
\[\begin{cases} \bar{x}_{k+1} = \varphi(x_k)\\ \bar{\bar{x}}_{k+1} = \varphi(\bar{x}_{k+1})\\ x_{k+1} = \bar{\bar{x}}_{k+1} - \frac{(\bar{\bar{x}}_{k+1}-\bar{x}_{k+1})}{\bar{\bar{x}}_{k+1} - 2\bar{x}_{k+1}+x_k} \end{cases} \]