Метод сходимости Ньютона не работает

Я пытаюсь аппроксимировать корень многочлена, используя метод Ньютона-Рафсона. Код, который я написал для этого, выглядит так:

#include <stdio.h>
#include <math.h>

int main (){

double c, nq, nnq, nnnq, v, h, q, o;

o=2;
c=-0.55;
nq=-0.04625;
nnq=-0.55;
nnnq=1;

   while(fabs(c-v) > 0.000001)
   {
      nq=((2*(o-1)+1)*(c)*(nnq)-(o-1)*nnnq)/((o-1)+1);                                                     
      q=(-o*c*nq+o*nnq)/(1-(c*c));
      h=(c-(nq/q));
      printf("The better root is %.15lf\n",h);
      v=c;
      c=h;
   }

}

Я знаю, что нет необходимости записывать переменные o, c, nq и т. д., потому что я мог бы просто использовать их точные значения. Это часть более крупной проблемы, и мне нужны эти переменные, так что не обращайте на это внимания.

Эта программа выводит это:

The better root is -0.578030303030303
The better root is -0.591696792857493
The better root is -0.598609887802599
The better root is -0.602171714355970
The better root is -0.604024260228500
The better root is -0.604992519745332
The better root is -0.605499890229896
The better root is -0.605766110042157
The better root is -0.605905895095070
The better root is -0.605979319651017
The better root is -0.606017894664121
The better root is -0.606038162857992
The better root is -0.606048812800124
The better root is -0.606054408979837
The better root is -0.606057349623975
The better root is -0.606058894866533
The better root is -0.606059706860161

Когда вместо этого он должен сходиться к точке -0,57735026918963. Я точно знаю, что Ньютон-Рафсон сходится, поэтому ошибка должна быть в коде. Я также пытался локализовать проблему с помощью printf, и я думаю, что проблема возникает во второй итерации. Я думаю, что программа не может правильно вычислить nq, но я не знаю, почему.


c while-loop newtons-method
person John Keeper    schedule 14.05.2018    source источник
comment
Одна серьезная проблема (хотя и не главная) заключается в том, что v не инициализирован. Это заставляет меня поверить, что вы компилируете без включенных предупреждений (или, возможно, просто игнорируете предупреждения)?   -  person Paul R    schedule 14.05.2018
comment
Вы можете добавить полином, который вы хотите вычислить, к вопросу?   -  person Afshin    schedule 14.05.2018
comment
@PaulR Правда, в другой программе у меня есть v = -20, чтобы убедиться, что цикл while запускается. В любом случае, выход одинаков.   -  person John Keeper    schedule 14.05.2018
comment
@Afshin Это многочлен Лежандра для n=2: (1/2)*(3 *х^2-1)   -  person John Keeper    schedule 14.05.2018
comment
Вы уверены, что первая производная верна? Мне кажется забавным.   -  person Bathsheba    schedule 14.05.2018
comment
@Bathsheba Да, для вычисления производных полиномов Лежандра вы можете использовать этот рекурсивный метод, зная что p_0=1 и p_1=x. В моем случае n=2.   -  person John Keeper    schedule 14.05.2018
comment
@ Вирсавия, да, мне тоже это показалось забавным. Вот почему я попросил уравнение.   -  person Afshin    schedule 14.05.2018
comment
@JohnKeeper, тогда вы используете 2 разных рекурсивных уравнения одновременно. Я думаю, именно поэтому вы не можете сходиться. Попробуйте использовать метод Ньютона только в своем уравнении (вручную добавьте дифференциал в код). Я думаю, вы должны сойтись в это время.   -  person Afshin    schedule 14.05.2018
comment
Почему вы написали ((o-1)+1) вместо o? Возможно, скобки неуместны.   -  person rici    schedule 14.05.2018
comment
@rici: Точно! Говорил же тебе, что первая производная была забавной.   -  person Bathsheba    schedule 14.05.2018

Ответы (2)


arrow_upward
1
arrow_downward

Вместо того, чтобы просто вычислять x = sqrt(1.0/3), вы хотите включить рекурсивную оценку полиномов Лежандра до порядка o, вероятно, для расширения метода позже до значений o больше 2. Итерация

P(0,c) = 1; P(1,c) = c;
(n+1)*P(n+1,c) = (2*n+1)*c*P(n,c) - n*P(n-1,c),   n=1, ... , o-1

и производная может быть вычислена как

(1-c^2)*P'(o,c) = -n*c*P(o,x) + n*P(o-1,c).

Эту итерацию вам нужно будет включить в цикл Ньютона, в идеальном случае используя объект для полинома Лежандра с методами для значения и производной. Я изменил вашу структуру для работы в JavaScript:

var my_div = document.getElementById("my_div");
var c = -0.55;
var v = -20;
var o = 2;
while( Math.abs(v-c) > 1e-12 ) {
    p0 = 1; 
    p1 = c;
    for(n=1; n<o; n++) {
        // p0, p1, p2 stand for P(n-1,c), P(n,c), P(n+1,c)
        p2 = ((2*n+1)*c*p1 - n*p0) / (n+1)
        p0 = p1; p1 = p2;
    }
    // p0, p1 now contain p(o-1,x), p(o,x)
    p1prime = ( -o*c*p1 + o*p0) / (1-c*c);
    h = c - p1/p1prime;
    my_div.innerHTML += "<p>The better root is "+h+"</p>";
    v = c; c = h;
 }
<div id="my_div"></div>

person Lutz Lehmann    schedule 22.05.2018

arrow_upward
1
arrow_downward

Это метод Ньютона для вашего уравнения (это быстрый код, не проверяйте имя переменной):

#include <stdio.h>
#include <math.h>

int main ()
{    
    double s = 2.0, fx = 0, dfx = 0, p = 0;

    while(fabs(s - p) > 0.000001)
    {
        fx = 0.5 * (3 * s * s - 1);
        dfx = 3 * s;
        p = s;
        s = s - (fx / dfx);
        printf("The better root is %.15lf\n", s);
    }    
    return 0;
}

и он сходился к 0.577350269189626. Ваша проблема в том, что вы пытаетесь вычислить 2 рекурсии одновременно. Кстати, в своем вопросе вы сказали, что хотите вычислить «корень многочлена». Я не совсем понял, что вы имели в виду. Если из корня вы имеете в виду квадратный корень вашего уравнения, вам необходимо обновить этот код и изменить fx и dfx соответственно.

person Afshin    schedule 14.05.2018
comment
> 0.000001 на самом деле нужно больше думать. - person Bathsheba; 14.05.2018
comment
@ Вирсавия, да, я просто копирую/вставляю его код и редактирую его. :) кроме того, важна и начальная точка (я начал с 2.0). - person Afshin; 14.05.2018