非线性方程求根方法总结&附代码（从二分法、试位法到牛顿迭代、二次插值等）

一只特立独行的蚂蚁

已于 2022-07-02 17:16:24 修改

阅读量1.2k

点赞数 2

分类专栏：数值分析算法文章标签：算法数值分析

于 2022-06-24 22:36:33 首次发布

本文链接：https://blog.csdn.net/luolei188/article/details/125453127

版权

数值分析算法专栏收录该内容

10 篇文章 20 订阅

订阅专栏

求方程 $f (x) = 0$ 的根。

1 划界法

函数通常在其根附近改变符号。利用两个初始值，这两个初始值必须位于跟的两侧。利用中值定理，必存在一个根落在区间内。利用这个性质，不断地缩小根的区间范围，直到足够接近结果为止。这类方法称为“划界法”。

1.1 二分法

思路：
若已知 $f (a) < 0$ 且 $f (b) > 0$ ，根据中值定理，必存在一个根落在区间(a, b)内。因此，

取区间(a, b)的中点，计算中点的值 $f (x) = f ((a + b) / 2)$ ；
若 $f ((a + b) / 2)$ 与 $f (b) < 0$ 同号，则将区间终点b替换成(a+b)/2，此时区间变成了 $(a, (a + b) / 2)$ ; 若 $f ((a + b) / 2)$ 与 $f (b) < 0$ 异号，则将区间起点a替换成(a+b)/2，此时区间变成了 $((a + b) / 2, b)$ ;
重复上述步骤1、2直到 $f (x)$ 接近为0，结束。

要求：已知两个点a，b的值 $f (a) < 0$ 且 $f (b) > 0$ 异号

收敛速度：线性收敛

缺点：收敛速度慢

伪代码：
在这里插入图片描述

1.1 试位法

在这里插入图片描述

实现思路与二分法类似，只不过二分法是直接取区间的中点，试位法利用式（5.7）得到一个x值，然后计算在该值下的 $f (x)$ ，其余步骤2、3与二分法相同。

要求：（同二分法一样）已知两个点a，b的值 $f (a) < 0$ 且 $f (b) > 0$ 异号

收敛速度：

缺点：在某些情况下收敛很慢，如下图。
在这里插入图片描述

1.3 改进试位的法思路

缓和试位法的片面性特征的一种方式是让算法检测一个边界是否固定不变。如果出现这种情况，可以将停滞的边界点处的函数值变为原来的一半。
注意，算法使用计数器来确定一个边界在两次迭代过程中是否保持不变，如果是，停滞的边界点的函数值将变为原来的一半。可以考虑当连续三次以上某一边界值保持不变的时候，就将该值减半。

1.4 Ridders方法

Ridders对试位法进行了改进。
在这里插入图片描述

1.5 划界法的问题

划界法（二分法和试位法等）实际应用时存在一个很大的问题就是，必须直到两个符号相反的初始值，有时候很难去找这样两个初始值。

2 开方法

2.1 定点迭代法

根据当前x的预测下一个x的值，从而不断的接近方程的根。
由 $f (x) = 0$ 得： $x = f (x) - x$
一个简单的预测方法就是： $x_{i+1} = f(x_i)-x_i$

收敛条件：令 $g (x) = f (x) - x$ ，当 $g^{'}(x)|<1$ 时收敛。
在这里插入图片描述

缺点：不一定收敛

2.2 牛顿迭代法（牛顿迭代法也是一种定点迭代法）

上面的定点迭代法预测的效果实在太差。
根据泰勒级数展开可得： $f(x_{i+1})=f(x_{i})+f^{'}(x_{i})(x_{i+1}-x_{i})+o(x^2)$
因此： $f(x_{i+1}) \approx f(x_{i})+f^{'}(x_{i})(x_{i+1}-x_{i})$
得：
$x_{i+1}=x_{i}- \frac{f(x_{i})}{f^{'}(x_{i})}$

在这里插入图片描述
最大的优点：二次收敛，收敛速度快！！！！

缺点：也有不收敛的情况。此外，需要函数能够准确计算一阶导数。并且，二重根得位置，收敛也很慢。

在这里插入图片描述

3 弦截法

3.1 弦截法

也称正割法。
牛顿迭代收敛很快，但是也存在一个问题，就是要能方便计算函数的一阶导数。而有些函数计算导数很困难，于是就想，能不能用近似的导数行不行，于是产生了一个想法，就是用有限差商来逼近导数。如下：
在这里插入图片描述
于是就得到了正割法：
$x_{i+1}=x_{i}- \frac{f(x_{i}) * (x_{i} - x_{i-1})}{f(x_{i})-f(x_{i-1})}$

3.2 弦截法与试位法区别

弦截法和试位法有一些相似之处，都是使用两个初始估计值来计算函数斜率的近似值，并将其投射到 x 轴，取得根的一个新估计值。然而，这两个方法的关键不同是新估计值如何取代两个初始值中的一个。在试位法中，根的最新估计值取代的是两个原估计值中函数值符号和 $f (x)$ 相同的一个。因此，这两个估计值总是将根界定在内。所以，对于所有现实情况，这个方法总是收敛的，因为根一直在划界范围内．相反，弦截法按照严格的顺序取代值，用新的 $x_{i+1}$ 的值取代 $x_{i}$ ，用 $x_{i}$ 取代 $x_{i-1}$ 因此，这两个值有时处于根的同一边。对于某些情形，这可能导致发散。
因此，用试位法优点在于一定收敛，优势有时也是劣势，试位法必须要已知两个符号相反的初始点；而弦截法不需要两个点的值必须符号相反，但劣势就是可能不收敛。
在这里插入图片描述

3.3 改进的弦截法

在这里插入图片描述

4 二次插值法

在上面的弦截法中，其实可以看做是用两个点做了一次线性插值，然后用插值出来的直线与x求交得到估计的x的值。因此，弦截法也可以称为线性插值法。
既然有线性插值法，那么就应该有非线性插值法。下面介绍二次插值法。

4.1 二次插值法

假设我们已经有了3个点，那么我们可以通过这3个点确定一个二次函数。
与线性插值法（弦截法）相同，这条抛物线与x轴的交点就表示新的估计值。

但是，由于二次函数不一定与x有交点，这就会导致算法失败。这是二次插值法的劣势。
因此如果与x轴没有交点的话，那就考虑取抛物线的最值点的x值作为替代。

优势：二次插值法有个很大的优势，就是收敛快。是仅次于牛顿迭代法之外的最快的方法之一。并且，对于二重根，收敛也很快。比牛顿迭代快。对于（三重根）多重根，没有测试过。
逆势：有时候会失败。譬如说当三个点的y值相等时，无法拟合成二次曲线。并且，当有抛物线与x轴有两个交点时，取哪一个比较麻烦。还待考虑。

4.2 逆二次插值法

二次插值法有一个问题就是，二次函数不一定与x有交点，并且可能存在三个y值相等。于是就想，要不插值成x关于y的二次函数。
在这里插入图片描述
这样，就不会出现与x轴没有交点的问题了。也不会出现三个x值相同的问题。

但是实际使用中发现，逆二次插值法面对二重根的时候，效果不太好。不如二次插值法好。

5 混合方法实现

由上面的介绍分析可以看出，各种球根的方法会有一些优势，也都会存在一些缺陷。那么，就想着怎么把这些方法的优势利用起来，合成一种综合性的方法。
要把划界法可靠性和开方法的速度优势结合起来。

5.1 分析(Van Wijngaarden和Dekker、Brent等探索结果)

虽然弦截法和试位法一般来讲比二分法收敛速度快，
但是在实践中发现，对于病态函数却是二分法的收敛速度要快些，如锯齿形函数、不连续函数，甚至是二阶导数在根附近有突变的光滑函数等。二分法总是将解区间对分一半，而弦截法和试位法有时要进行很多次迭代才能把两个相距较远的端点向根所在的位置拉近。
Ridders方法虽然很好，但有时反被复杂化。
究竟有没有一种方法，既能超线性收敛、又不失二分法的确定性呢？答案是肯定的。
可以对一个假设为超线性收敛的方法进行跟踪和记录，看它是否真正地按假定的方式收敛。如果不是的话，可以在迭代过程中穿插一些二分法的步骤，以保证这种方法至少具有线性收敛性。这类堪称上策的方法，需要注意记录迭代过程中的细节，并需仔细考虑舍人误差对指导策略的影响，而且还必须清楚知道算法何时收敛。
20世纪60年代，阿姆斯特丹数学中心的Van Wijngaarden和Dekker等学者成功地研究出一种能够实现上述要求的算法，而后该算法又由 Brent 进行了改进。为简单起见，我们称该算法的最后版本为 Brent方法。只要函数在包含根的初始区间内是可计算的， Brent 方法即可保证收敛（这个结论已由 Brent 证明）。Brent 方法将根的划界、二分法以及二次反插值方法相结合，从过零附近的邻域进行迭代收敛。

5.2 布伦特法（Brent）

这就会有很多种组合了。
一种组合就是布伦特方法。思路是尽可能的使用速度较快的开方法（逆二次插值法）；如果生成了不可接受的结果，就使用二分法。
一般的，这个过程中，刚开始二分法会占优势（使用比较多）；当距离根比较接近时，开方法会比较占优势，接近根的速度比较快。
在这里插入图片描述

5.3 “牛顿迭代+”方法

前面Brent方法并没有使用一阶导数的性质，如果能够方便并准确的计算一阶导数，可以使用牛顿迭代和二分法结合起来的方法。
同样的，尽可能地使用牛顿迭代法，如果不行就是用二分法。思路方法同Brent方法差不多。

由于牛顿迭代收敛极快！！！这是它天然的优势，必须用起来。但是也要注意到，牛顿迭代法可能不收敛。
但是不收敛的根本原因也可以归结为：给的初始点的位置不好。如果初始点位置给的好，或者足够接近，牛顿迭代一定会收敛。
但是，遇到二重根、多重根的时候，牛顿迭代收敛会很慢，成线性收敛。

5.4 牛顿迭代法和分形

要考虑收敛到复根区域的问题。在此不做深入介绍，感兴趣的可以去看看相关的介绍。
在这里插入图片描述

Ref 代码实现

Ref1 牛顿迭代法

/**
     * 牛顿迭代法，可能不收敛
     * 备注：迭代收敛较快，但是对于多重根，迭代非常慢，需要改进；并且容易陷入局部极值点位置
     * @param func 目标函数
     * @param dvtFunc 计算目标函数一阶微分的函数
     * @param x0 初值
     * @param tol
     */
    function iterationNewton(
        func: (param: number) => number,
        dtFunc: (param: number) => number,
        x0: number,
        tol: number = 1e-12,
        validEps?: number,
        maxIterNum?: number,
    ) {
        let iter = 0;
        let xi: number = x0;
        let bIsDecrease = true;
        while (iter < NORMAL_ITER_NUM || bIsDecrease) {
            const fx = func(xi);
            const dfx = dtFunc(xi);
            const newx = xi - fx / dfx;
            if (Math.abs(newx - xi) < tol || iter > (maxIterNum || MAX_ITER_NUM)) {
                xi = newx;
                break;
            }

            if (iter >= NORMAL_ITER_NUM) {
                const newFx = func(newx);
                bIsDecrease = Math.abs(newFx) < Math.abs(fx) - 1e-12; // 如果迭代趋势收敛(距离0更近了)，继续迭代
            }

            xi = newx;
            iter++;
        }

        if (Math.abs(func(xi)) < (validEps || tol)) {
            return xi;
        }
        return undefined;
    }

Ref2 改进的截弦法

/**
     * 改进的正割法（近似牛顿迭代法，可能不收敛）
     * 备注：迭代收敛较快，但是对于多重根，迭代非常慢，需要改进；并且容易陷入局部极值点位置
     * @param x0 初值
     * @param tol
     */
    function linearInterpolation(
        func: (param: number) => number,
        x0: number,
        tol: number = 1e-12,
    ) {
        let iter = 0;
        let xi: number = x0;
        let bIsDecrease = true;
        const delta = 0.0001;
        while (iter < NORMAL_ITER_NUM || bIsDecrease) {
            const fx0 = func(xi);
            const fx1 = func(xi + delta); // 多重根的时候，delta需要更小，否则导数计算不准确，会导致不收敛
            const newx = xi - (fx0 * delta) / (fx1 - fx0);
            if (Math.abs(newx - xi) < tol || iter > MAX_ITER_NUM) {
                xi = newx;
                break;
            }

            if (iter >= NORMAL_ITER_NUM) {
                bIsDecrease = Math.abs(func(newx)) < Math.abs(func(xi)) - 1e-12; // 如果迭代趋势收敛(距离0更近了)，继续迭代
            }

            xi = newx;
            iter++;
        }

        if (Math.abs(func(xi)) < tol) {
            return xi;
        }
        return undefined;
    }

Ref3 二次插值法（米勒法）

/**
     * 二次插值法（米勒法），求方程的根
     * 备注：稳定性不错，但是给的初值不好，容易陷入局部极值点附近无意义循环
     * @param x0 初值
     * @param h 初始步长
     * @param tol
     */
    function quadraticInterpolation(
        func: (param: number) => number,
        xMin: number,
        xMax: number,
        tol: number = 1e-12,
    ): number | undefined {
        let iter = 0;
        let x3 = xMax;
        const xs: number[] = [xMin, xMax, (xMin + xMax) / 2];
        const fxs: number[] = [func(xs[0]), func(xs[1]), func(xs[2])];
        let bIsDecrease = true;
        while (iter < NORMAL_ITER_NUM || bIsDecrease) {
            // 三个点(param0, fx0)，(param1, fx1)，(param2, fx2)二次方程求根g(x) = 0
            // const x0x1Diff = xs[0] - xs[1];
            // const x1x2Diff = xs[1] - xs[2];
            // const x2x0Diff = xs[2] - xs[0];
            // const gx = (fx0 * (x - x1) * (x - x2)) / (x0x1Diff * -x2x0Diff)
            //          + (fx1 * (x - x0) * (x - x2)) / (x1x2Diff * -x0x1Diff)
            //          + (fx2 * (x - x0) * (x - x1)) / (x2x0Diff * -x1x2Diff);
            // 得到g(x) = ax^2 + bx + c
            const h0 = xs[1] - xs[0];
            const h1 = xs[2] - xs[1];
            const d0 = (fxs[1] - fxs[0]) / h0;
            const d1 = (fxs[2] - fxs[1]) / h1;
            const a = (d1 - d0) / (h0 + h1);
            const b = a * h1 + d1;
            const c = fxs[2];

            const bsqr = b * b;
            const ac4 = 4 * a * c;
            if (bsqr - ac4 < 0) {
                // return undefined; // 多重根求不到：需要处理复数根的情况
                x3 = -b / (2 * a);
                if (Math.abs(x3 - xs[2]) < tol) {
                    break;
                }
            } else {
                let delta = Math.sqrt(bsqr - ac4);
                delta = b > 0 ? delta : -delta;

                const den = -(2 * c) / (b + delta);
                x3 = xs[2] + den;

                if (Math.abs(den) < tol) {
                    break;
                }
            }

            if (Math.abs(x3 - xs[1]) < Math.abs(x3 - xs[0])) {
                [xs[0], xs[1], xs[2]] = [xs[1], xs[2], x3];
                [fxs[0], fxs[1], fxs[2]] = [fxs[1], fxs[2], func(x3)];
            } else {
                [xs[1], xs[2]] = [xs[2], x3];
                [fxs[1], fxs[2]] = [fxs[2], func(x3)];
            }

            if (iter >= NORMAL_ITER_NUM) {
                bIsDecrease = Math.abs(fxs[2]) < Math.abs(fxs[1]) - 1e-12; // 如果迭代趋势收敛(距离0更近了)，继续迭代
            }
            if (iter > MAX_ITER_NUM) {
                break;
            }

            iter++;
        }

        if (Math.abs(func(x3)) < tol) {
            return x3;
        }
        return undefined;
    }

Ref4 逆二次插值法

/**
     * 逆二次插值法，求方程的根
     * 备注：能保证能计算一个根，即使初值在局部极值点附近也能迭代最后到一个根（如果存在根），但是对于多重根的情况，（y值相等）会使用割线法，计算效率较低，并且达到的精度不够
     * @param x0 初值
     * @param h 初始步长
     * @param tol
     */
    function inverseQuadraticInterpolation(
        func: (param: number) => number,
        x0: number = 0.5,
        h: number,
        tol: number = 1e-12,
    ): number | undefined {
        let iter = 0;
        let x3 = x0;
        const xs: number[] = [x0 - h, x0 + h, x0];
        const fxs: number[] = [func(xs[0]), func(xs[1]), func(xs[2])];
        let bIsDecrease = true;
        while (iter < NORMAL_ITER_NUM || bIsDecrease) {
            // 如果fx0，fx1，fx2其中两个相等，采用正割法计算第三个点
            if (Math.abs(fxs[0] - fxs[1]) < tol || Math.abs(fxs[0] - fxs[2]) < tol) {
                // 使用(x1, fx1)，(x2, fx2)正割法
                x3 = xs[2] - (fxs[2] * (xs[1] - xs[2])) / (fxs[1] - fxs[2]);
                if (Math.abs(x3 - xs[2]) < tol) {
                    break;
                }

                [xs[0], xs[1], xs[2]] = [xs[1], xs[2], x3];
                [fxs[0], fxs[1], fxs[2]] = [fxs[1], fxs[2], func(x3)];
                iter++;
                continue;
            } else if (Math.abs(fxs[1] - fxs[2]) < tol) {
                // 选取(x0, fx0)，(x1, fx1)使用正割法？？？还是选取(x0, fx0)，(x2, fx2)使用正割法？？？怎么选取问题
                x3 = xs[2] - (fxs[2] * (xs[0] - xs[2])) / (fxs[0] - fxs[2]);
                if (Math.abs(x3 - xs[2]) < tol) {
                    break;
                }

                [xs[1], xs[2]] = [xs[2], x3];
                [fxs[1], fxs[2]] = [fxs[2], func(x3)];
                iter++;
                continue;
                // if (Math.abs(xs[0] - xs[2]) < Math.abs(xs[0] - xs[1])) {
                //     // 选取(x0, fx0)，(x2, fx2)使用正割法
                //     x3 = xs[2] - (fx2 * (xs[0] - xs[2])) / (fx0 - fx2);
                //     if (Math.abs(x3 - xs[2]) < tol) {
                //         break;
                //     }

                //     [xs[1], xs[2]] = [xs[2], x3];
                //     iter++;
                //     continue;
                // } else {
                //     // 选取(x0, fx0)，(x1, fx1)使用正割法
                //     x3 = xs[1] - (fx1 * (xs[0] - xs[1])) / (fx0 - fx1);
                //     if (Math.abs(x3 - xs[1]) < tol) {
                //         break;
                //     }

                //     xs[2] = x3;
                //     iter++;
                //     continue;
                // }
            }

            // 三个点(param0, fx0)，(param1, fx1)，(param2, fx2)逆二次方程x = g(y)求根g(y) = 0：保证曲线总与x轴有交
            const fx0fx1Diff = fxs[0] - fxs[1];
            const fx1fx2Diff = fxs[1] - fxs[2];
            const fx2fx0Diff = fxs[2] - fxs[0];
            // // 逆二次插值函数
            // const gy = (x0 * (y - fx1) * (y - fx2)) / (fx0fx1Diff * -fx2fx0Diff)
            //          + (x1 * (y - fx0) * (y - fx2)) / (fx1fx2Diff * -fx0fx1Diff)
            //          + (x2 * (y - fx0) * (y - fx1)) / (fx2fx0Diff * -fx1fx2Diff);
            // 令y = 0得到x3:
            x3 =
                (xs[0] * fxs[1] * fxs[2]) / (fx0fx1Diff * -fx2fx0Diff) +
                (xs[1] * fxs[0] * fxs[2]) / (fx1fx2Diff * -fx0fx1Diff) +
                (xs[2] * fxs[0] * fxs[1]) / (fx2fx0Diff * -fx1fx2Diff);

            // 容差多大还待考虑
            if (Math.abs(x3 - xs[2]) < tol) {
                break;
            }

            if (Math.abs(x3 - xs[1]) < Math.abs(x3 - xs[0])) {
                [xs[0], xs[1], xs[2]] = [xs[1], xs[2], x3];
                [fxs[0], fxs[1], fxs[2]] = [fxs[1], fxs[2], func(x3)];
            } else {
                [xs[1], xs[2]] = [xs[2], x3];
                [fxs[1], fxs[2]] = [fxs[2], func(x3)];
            }

            if (iter >= NORMAL_ITER_NUM) {
                bIsDecrease = Math.abs(fxs[2]) < Math.abs(fxs[1]) - 1e-12; // 如果迭代趋势收敛(距离0更近了)，继续迭代
            }
            if (iter > MAX_ITER_NUM) {
                break;
            }

            iter++;
        }

        if (Math.abs(func(x3)) < tol) {
            return x3;
        }
        return undefined;
    }