0

この単純なルート検索アルゴリズムを実装します。 http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method 実装の何が問題なのか、一生わからない。根は膨らみ続け、収束の兆しはありません。助言がありますか？

ありがとう。

#include <iostream>
#include <complex>

using namespace std;

typedef complex<double> dcmplx;

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

    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 = 20;
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-q0)*(p0-r0)*(p0-s0));
q = q0 - f(q0)/((q0-p)*(q0-r0)*(q0-s0));
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;
}
4

3 に答える 3

1

ああ、問題は、N 次多項式の係数を次のように指定する必要があることでした。

1*x^N + a*x^(N-1) + b*x^(N-2) ... etc + z;

ここで、1 は最大次数項の係数です。そうしないと、最初のルートは決して収束しません。

于 2013-11-10T17:08:11.630 に答える
1

半年遅れ: 謎の解決策は、分母が多項式の導関数の近似値である必要があるため、先頭の係数 a4 を係数として含める必要があることです。

別の方法として、return ステートメントで多項式の値を a4 で除算して、多項式が効果的に正規化されるようにすることもできます。つまり、先頭の係数が 1 になります。

Bo Jacoby によるウィキペディアのサンプル コードは、メソッドの Seidel 型の変形であることに注意してください。古典的な定式化は、すべての新しい近似が古い近似から同時に計算される Jordan のような方法です。ザイデルは、多次元ニュートン法としての定式化がヤコビに提供する次数 2 よりも高速に収束する可能性があります。

ただし、次数が大きい場合は、高速多項式乗算アルゴリズムを使用して多項式値​​と分母の積の必要な多点評価を行うことで、ヤコビを高速化できます。

于 2014-04-07T18:07:25.057 に答える