求非线性方程根的 Durand-Kerner 方法答案

【问题标题】：Durand-Kerner method to find roots of nonlinear equation求非线性方程根的 Durand-Kerner 方法
【发布时间】：2014-12-14 12:39:55
【问题描述】：

我被要求找到 f(x) = 5x(e^-mod(x))cos(x) + 1 的根。我之前使用过 Durand-Kerner 方法来查找函数 x^4 -3x^3 + x^2 + x + 1 的根，代码如下所示。我以为我可以简单地重用代码来查找 f(x) 的根，但是每当我用 f(x) 替换 x^4 -3x^3 + x^2 + x + 1 时，程序都会为所有根输出 nan。我的 Durand-Kerner 实现有什么问题，我该如何修改它以适用于 f(x)？如果有任何帮助，我将不胜感激。

#include <iostream>
#include <complex>
#include <math.h>

using namespace std;

typedef complex<double> dcmplx;

dcmplx f(dcmplx x)
{
    // the function we are interested in
    double a4 = 1;
    double a3 = -3;
    double a2 = 1;
    double a1 = 1;
    double a0 = 1;

    return (a4 * pow(x,4) + a3 * pow(x,3) + a2 * pow(x,2) + a1 * x + a0);
}


int main()
{

dcmplx p(.9,2);
dcmplx q(.1, .5);
dcmplx r(.7,1);
dcmplx s(.3, .5);

dcmplx p0, q0, r0, s0;

int max_iterations = 100;
bool done = false;
int i=0;

while (i<max_iterations && done == false)
{
    p0 = p;
    q0 = q;
    r0 = r;
    s0 = s;


p = p0 - f(p0)/((p0-q)*(p0-r)*(p0-s));
q = q0 - f(q0)/((q0-p)*(q0-r)*(q0-s));
r = r0 - f(r0)/((r0-p)*(r0-q)*(r0-s0));
s = s0 - f(s0)/((s0-p)*(s0-q)*(s0-r));

    // if convergence within small epsilon, declare done
    if (abs(p-p0)<1e-5 && abs(q-q0)<1e-5 && abs(r-r0)<1e-5 && abs(s-s0)<1e-5)
        done = true;

    i++;
}

cout<<"roots are :\n";
cout << p << "\n";
cout << q << "\n";
cout << r << "\n";
cout << s << "\n";
cout << "number steps taken: "<< i << endl;

return 0;
}

到目前为止，我唯一更改的是 dcmplx f 函数。我一直把它改成

dcmplx f(dcmplx x)
{
    // the function we are interested in
    double a4 = 5;

    double a0 = 1;

    return (a4 * x * exp(-x) * cos(x) )+ a0;
}

【问题讨论】：

标签： c++ polynomial-math nonlinear-functions newtons-method

【解决方案1】：

您使用的 Durand-Kerner 方法要求函数在您工作的时间间隔内是连续的。

在这里，我们发现数学观点与数值应用的限制之间存在差异。我建议你绘制你的函数（在谷歌中输入公式会给你一个快速的概览，当然对于真实的部分）。你会注意到：

由于余弦的周期性，存在无穷多个根。
由于 x*exp(-x)，绝对值迅速上升，超过了浮点数可以保持的最大精度。

为了了解对代码的影响，我邀请您跟踪不同的迭代。您会注意到 p、r 和 s 收敛得非常快，而 q 正在发散（显然在一个巨大的峰值的轨道上）：

在第二次迭代中 q 已经在 1e74
在第 3 次迭代中，已经超出了 double 可以存储的容量。
由于 q 用于计算 p、r 和 s，误差会传播到其他项
在第 5 次迭代中，所有项都在 NAN
然后它勇敢地继续通过 100 次迭代

也许您可以通过选择不同的起点来使其发挥作用。如果没有，您将不得不使用其他方法并仔细选择您正在使用的中间墙。

【讨论】：

知道我还可以使用什么其他方法吗？
我不是数学家。但我会尝试简单的 Newton-Raphson ( Xi+1 = Xi - f(Xi)/f'(Xi) )。或者布伦特法en.wikipedia.org/wiki/Brent%27s_method
我确认 Brent 方法和维基百科中公开的算法运行良好，只要您可以在根节点上方和下方给出一个点。我尝试过使用双打而不是复杂，它在 5 次迭代中找到了根源

【解决方案2】：

您应该在有关 Durand-Kerner 方法（由 Karl Weierstrass 于 1850 年左右发明）的文档中指出它仅适用于 多项式。您的第二个函数远非多项式。

确实，由于 mod 函数，它必须被声明为数值方法的讨厌函数。它们中的大多数依赖于给定函数的连续性，即，如果值接近于零，则附近很有可能存在根，如果符号在区间内发生变化，则区间中有根。即使是最基本的无导数方法，例如该类复杂端的二分法或 Brents 方法，也预先假定了这些属性。

【讨论】：