【数值计算】数值解析--非线性方程的根

线性方程与非线性方程

当我们求关于的方程的解时，如果，是像

这样的线性形方程(1次方程)的话，其解为，

这里的。但是，是非线性方程的时候解法要复杂的多。比如，像下面这样的次代数方程(algebraic equation)的情况，

的2次方程我们很容易求解，3次或4次方程可以通过Cardano公式或者Ferrari公式求解，然而5次以上却无法直接求解。

更不用说，还有像超越方程式(transcendental equation)这样的无穷次的代数方程()。例如，用无穷幂级数表示的形式便是超越方程。超越方程一般无法直接求解，只能求近似解。

本文，我们对如下的非线性方程的数值求解方法进行逐一介绍。

• 2分法
• 牛顿-拉弗森方法
• 多次元牛顿-拉弗森方法
• 霍纳方法
• DKA法

2分法 

假设是在区间上的函数。和符号相反的时候()，区间内至少存在一个解。我们先求区间的中点的函数值，如果这个值的符号与同号， 则在内存在解。 这样的话，解的存在区间便变成了原来的。同样，再使用的中点的函数值，解的存在区间就变成了。把上述步骤反复处理求得方程式的解的方法就叫做2分法(bisection method)。

图1 2分法

图1为2分法求近似解的说明图。图中解的存在区间为。用对进行初始赋值。 这时的中点为。 由于图1，解在内存在， 令，这时的中点为， 由图1，解在内存在。上述步骤反复迭代，求。

2分法的处理顺序如下：

1. 令, 。
2. 求中点的函数值。
3. 的时候，时，令。
4. 返回步骤2。

迭代终止的收敛条件如下：

下面为2分法的代码示例。

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

 

/*! * 2分法(bisection method) * @param[in] func  * @param[in] xl,xr 检索范围 * @param[out] x 解 * @param[inout] max_iter 最大迭代次数(循环终了后,返回实际迭代次数) * @param[inout] eps 允许误差(循环终了后,返回实际的误差值)  * @return 1:成功,0:失败 */int bisection(double func(const double), double xl, double xr, double &x, int &max_iter, double &eps){    double f = func(xl);    double fmid = func(xr);     if(f*fmid >= 0.0) return 0.0;     int k;    double dx = fabs(xr-xl), xmid;    for(k = 0; k < max_iter; ++k){        xmid = 0.5*(xl+xr);    // 中点        dx *= 0.5;         // 中点的函数值        fmid = func(xmid);         // 收敛判定        if(dx < eps || fmid == 0.0){            x = xmid;            max_iter = k; eps = dx;            return 1;        }         // 新的区间        if(f*fmid < 0){            xr = xmid;        }        else{            xl = xmid;            f = fmid;        }    }     max_iter = k; eps = dx;    return 0;}

牛顿-拉弗森方法

用2分法虽然可以求出指定区间内的解，但是收敛太慢。如果是区间内连续可微的函数的话，用函数的梯度值可以更快收敛。

解的初期值设为。在 的函数的切线为，

这条切线与相交()的点为，

把作为下一次的近似值，用同样方法逐渐向解趋近。上述处理顺序如图2所示。这种方法便称作牛顿-拉弗森法(Newton-Raphson method)，或者牛顿法(Newton method)。牛顿法的公式如下：

迭代次数足够多的话会无限趋近的解。 迭代终止的收敛条件为，

牛顿法只指定一个初始值，并且收敛速度比2分法快很多。然而，这种方法需要求及其微分值，而且根据设定的初始值，也存在解无法收敛的情况。

图2 牛顿法

牛顿法代码示例如下：

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

 

/*! * 牛顿法(Newton-Raphson method) * @param[in] func  * @param[inout] x 设定初始值，之后反复迭代 * @param[inout] max_iter 最大迭代数(迭代终了,返回实际迭代数) * @param[inout] eps 允许误差(迭代终了，返回实际误差)  * @return 1:成功,0:失败 */int newton(double func(const double), double dfunc(const double),            double &x, int &max_iter, double &eps){    double f, df, dx;    int k;    for(k = 0; k < max_iter; ++k){        f = func(x);        df = dfunc(x);        x = xn-f/df;         // 收敛判定        dx = fabs(f/df);        if(dx < eps || fabs(f) < eps){            max_iter = k; eps = dx;            return 1;        }    }     max_iter = k; eps = dx;    return 0;}




多次元牛顿-拉弗森方法

把牛顿法表示成联立非线性方程的一般化形式：

向量表示为，

这里的，

第步的近似值表示成，把上式用 的附近点进行泰勒展开：

这里的是的雅可比行列。 忽略掉2次以上的项，联立非线性方程变为：

设，可得

这个式子是以为未知数的线性联立方程，通过LU分解等解法可以得到。 然后再按照下式计算。

例)2元联立非线性方程

然后，

由此，与相关的公式如下：

求解，得：

然后再更新。

霍纳法

首先我们来计算一下用2分法或牛顿法解下述代数方程的运算次数，

计算的值时，乘法运算次，加法运算次，乘法运算次数呈的2次增加，比如，时55次，时5050次。

为了减少乘法运算次数，对原式做如下变形：

此时，从最里面的括号内开始计算的话，可知

最后。 此时的乘法运算次数是次，加法运算也是次。当很大时，计算次数可以大幅减小。例如当乘法运算次数时为10次，时为100次。这种计算方法便被称作霍纳(Horner)。

用霍纳法求代数方程的值及导函数的代码示例如下。

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

 

/*! * 霍纳法求代数方程的值 * @param[in] x 变量 * @param[in] b 系数 * @param[in] n 方程的次数 * @return 代数方程式的值 */template<class T>inline T func_h(double x, const vector<T> &b, int n){    T f = b[0];    for(int i = 1; i <= n; ++i){        f = b[i]+f*x;    }    return f;}/*! * 导函数计算 * @param[in] x 变量 * @param[in] b 系数 * @param[in] n 方程的次数 * @return 代数方程式的导函数值 */template<class T>inline T dfunc_h(double x, const vector<T> &b, int n){    T df = n*b[0];    for(int i = 1; i <= n-1; ++i){        df = (n-i)*b[i]+df*x;    }    return df;}

示例函数使用了复数计算，可以对应各种类型的方程。

DKA法 

代数方程包括重根在内有复数解(实数也是复数的一部分)。以上的方法，通过赋予合适的初始值和范围求得了一个解。当求代数方程的复数解时不得不重新设定初始值和检索范围。以下对求个复数解的方法进行讲解。

首先，有复数解的代数方程用未知数的表示如下：

用牛顿法分析代数方程。上述方程的解表示成，方程发生如下变形：

将其微分，

当等于时，由于，上式只剩下第一项。同理时左边仅剩第项。用公式表示，即：

这里的。由牛顿法公式，用的各近似值逼近，

这里的。这个式子就是Weierstrass逼近法，或者，称作Durand-Kerner法[1]。以下把上式简称DK公式。

虽然我们可以通过DK公式计算个复数解，由于使用了牛顿法，有必要选择一个合适的的初始值。这里，我们使用下面的(Aberth)的初始值。

这里的，是复数平面上包含这个解的圆的半径。这里的为。

我们把使用Aberth的初始值，通过DK公式的迭代计算，获得个根的方法叫做DKA(Durand-Kerner-Aberth)法。

DKA法的处理顺序

1. 计算中心，半径及初始值
2. 用DK公式更新，在收敛前反复执行2，直到获得近似解。

关于Aberth初始值半径的求法有很多种，这里介绍的是最初Aberth提议的方法[2]。 令得到以下方程。

这里的是通过霍纳法获得的系数。  且所有系数非零时，的解在，求

的方程时以解为半径的圆内存在。由此，在内有解。使用计算的初始值的半径。 由于有单一的实数解，所以使用牛顿法求解。

使用DKA法计算多个根的代码如下。首先是计算Aberth初始值的函数。

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

 

/*! * 使用Aberth法计算初始值 * @param[out] z 初始值 * @param[in] c 多项式系数 * @param[in] n 方程的次数 * @param[in] max_iter 半径计算用最大迭代次数 * @param[in] eps 半径计算用的允许误差 * @return 1:成功,0:失败 */int aberth(vector<complexf> &z, const vector<complexf> &c, int n, int max_iter, double eps){    // 系数    vector<complexf> a;    complexf zc = -c[1]/(c[0]*(double)n);    horner(c, a, n, zc);     vector<double> b(a.size());    b[0] = abs(a[0]);    for(int i = 1; i <= n; ++i){        b[i] = -abs(a[i]);    }     // 通过牛顿法计算Aberth的初始值的半径    double r = 100.0;    newton(b, n, r, max_iter, eps);     // Aberth的初始值    for(int j = 0; j < n; ++j){        double theta = (2*RX_PI/(double)n)*j+RX_PI/(2.0*n);        z[j] = zc+r*complexf(cos(theta), sin(theta));    }     return 1;}

这里由于涉及到了复数计算，所以使用C++中的complex库。

#include <complex>// 复数typedef complex<double> complexf;

计算半径时，的系数计算用到的霍纳函数如下。

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18

 

/*! * 求霍纳法的系数值 * @param[in] a 代数方程的系数 * @param[out] b x=x0时系数 * @param[in] x0  */template<class T>inline void horner(const vector<T> &a, vector<T> &b, int n, T x0){    if(n <= 2) return;     b = a;    for(int i = 0; i <= n; ++i){        for(int j = 1; j <= n-i; ++j){            b[j] += x0*b[j-1];        }    }}

得到初始值后，利用以下的Weierstrass函数求解。

  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

 

/*! * Weierstrass法(DK公式) * @param[inout] z 获得初始值，返回解 * @param[in] c 多项式的系数 * @param[in] n 方程的次数 * @param[inout] max_iter  * @param[inout] eps  * @return 1:成功,0:失败 */int weierstrass(vector<complexf> &z, const vector<complexf> &c, int n, int &max_iter, double &eps){    double e = 0.0, ej;     vector<complexf> zp;    complexf f, df;    int k;    for(k = 0; k < max_iter; ++k){        zp = z;         // DK式的计算        for(int j = 0; j < n; ++j){            f = func(z[j], c, n);            df = c[0];            for(int i = 0; i < n; ++i){                if(i != j){                    df *= zp[j]-zp[i];                }            }             z[j] = zp[j]-f/df;        }         // 误差计算        e = 0.0;        for(int j = 0; j < n; ++j){            if((ej = abs(func(z[j], c, n))) > e){                e = ej;            }        }         // 收敛判定        if(e < eps){            max_iter = k; eps = e;            return 1;        }    }     eps = e;    return 0;}




关于Aberth的初始值

下面推导一下代数方程的解与系数之间的关系。首先，以2次方程为例。

使用解把方程改写成如下形式。

与之前的系数做对比，得到下面的关系式。

同理，n次方程表示如下，

于是，n次方程的解与系数的关系为，

这里的。如，当方程为3次方程时()，

由时的关系式，

Aberth的初始值中，把解所在复数平面的重心：

当作中心的半径为的圆周上等间距的配置。也就是说，

由欧拉公式()，可得：

这里的为。另外，角度值加上是为了防止，初始值在实轴上对称分布(防止成为共轭复数)。

[1] http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method
[2] O. Aberth, ``Iteration methods for finding all zeros of a polynomial simultaneously'', Math. Comput. 27, pp.339-344, 1973.