Что случилось с моей реализацией durand-kerner?

Question

Реализация этого простого алгоритма поиска корней. 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; 
} 

Answer 1

Половина года: решение загадки заключается в том, что знаменатель должен быть приближением производной от полинома и, следовательно, должен содержать ведущий коэффициент a4 как фактор.

В качестве альтернативы, можно разделить полиномиальное значение на а4 в ответном заявлении, так что многочлен эффективно нормированное, то есть, имеет старший коэффициент 1.

Обратите внимание, что пример кода в википедии Бо Джейкоби является Вариант метода сейделского типа, классическая формулировка - метод, подобный Иордану, где все новые аппроксимации одновременно вычисляются из старого приближения. Зейдел может иметь более быструю сходимость, чем порядок 2, что формулировка как многомерный метод Ньютона обеспечивает Якоби.

Однако для больших степеней Якоби можно ускорить с использованием алгоритмов быстрого полиномиального умножения для требуемых многоточечных оценок полиномиальных значений и продуктов в знаменателях.

Answer 2

Вы неправильно применили формулы. Например

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

должен быть

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

Посмотрите еще раз на вики-статьи

Answer 3

Ах, проблема заключалась в том, что коэффициенты в N-степени полинома должны быть определены как

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

где 1 - коэффициент наибольшего члена степени. В противном случае первый корень никогда не сходится.