Метод Рунге-Кутта - не получается реализовать для системы.

[JustAMan]

Реализовал и опробовал для одного уравнения. А для системы уравнений он у меня не те числа выдаёт, что приводятся в примере :(

Мне требуется посчитать значения этой системы уравнений:
x_1'=x_2
x_2'=-0.2*x_2 - 10*sin(x_1)
При начальных условиях:
x_1(0) = 0.3; x_1'(0)=0
Скорее всего, я именно начальные условия неправильно задаю. Посмотрите пожалуйста.. как-то иначе задавать их тут нужно?

Ниже я привожу свой код, реализующий это. Думаю, тут ошибку заметят и люди не понимающие в си-шном коде, поскольку он интуитивно понятен.. Ошибка именно в математике. Сомневаюсь, что в алгоритме :)

double myF( unsigned int n, double x1, double x2 ){
   double result = 0;
   if ( n == 0 ) {
      result = x2;
   } else if ( n == 1 ) {
      result = -0.2 * x2 - 10 * sin( x1 );
   }   
   return result;   
}

void RungeKutt() {
   
   unsigned int i;
   double x0, x1, x2[ 2 ], x2_t[ 2 ], k1[ 2 ], k2[ 2 ], k3[ 2 ], k4[ 2 ];
   
   double h = 0.1;
   x0 = 0.0;
   x1 = 0.0;
   x2[ 0 ] = 0.3;

   for( i = 0; i < 4; i++ ){

      x1 = x0 + i*h;            
      for( unsigned int n = 0; n < 2; n++ ){
      
         k1[ n ] = h * myF( n, x1, x2[ n ] );
         k2[ n ] = h * myF( n, x1 + ( h / 2 ), x2[ n ] + ( k1[ n ] / 2 ) );
         k3[ n ] = h * myF( n, x1 + ( h / 2 ), x2[ n ] + ( k2[ n ] / 2 ) );
         k4[ n ] = h * myF( n, x1 + ( h / 2 ), x2[ n ] + k3[ n ] );
         
         x2_t[ n ] = ( k1[ n ] + 2 * k2[ n ] + 2 * k3[ n ] + k4[ n ] ) / 6;

         printf( "%i: x1=%f, x2=%f, k1=%f, k2=%f, k3=%f, k4=%f, delta_x2=%f\n", n, x1, x2[ n ], k1[ n ], k2[ n ], k3[ n ], k4[ n ], x2_t[ n ] );
         
         x2[ n ] = x2[ n ] + x2_t[ n ];
         
      }
   }
   
}

[Circiter]

Код изобилует ошибками. Предлагаю вам для начала переписать его аккуратно и без странностей. :)

[ewert]

JustAMan в сообщении #399143 писал(а):

Ошибка именно в математике. Сомневаюсь, что в алгоритме :)

Тут математика и алгоритм неразделимы. Принципиально, что функция myF должна быть тоже векторной, а не скалярной, как у Вас. Соответственно, вспомогательные переменные k1, k2, k3 и k4 также должны быть векторными.

[JustAMan]

А можно в терминах "попроще" объяснить в чём тут ошибки? :)
У нас ведь в системе двух уравнений должно быть два набора переменных k1, k2, k3, k4, правильно? Соответственно для каждой рассчитывается своё приращение, которое у меня именуется: x2_t.

Видимо с реализацией самой функции что-то совсем не так? (ну и с её вызовом - это тогда уж понятно). У нас ведь для каждого уравнения подсчитывается свой набор k1, k2, k3, k4? Поэтому я и разделил так в функции этот участок..

Circiter
А почему он неаккуратно написан? :) Я бы мог, конечно, через двумерные массивы это всё описать - но так сделал только потому, что это у меня работало для случая одного уравнения, ну и чтобы быстро как-то переделать в случай с системой - я пошёл вот по такому методу написания кода) Да и это напорядок усложнило бы сейчас понимание алгоритма..

[Circiter]

Хотя-бы разберитесь с вашими индексами и размерами массивов...

[JustAMan]

Circiter в сообщении #399436 писал(а):

Хотя-бы разберитесь с вашими индексами и размерами массивов...

А что там разбираться? Не понимаю.. Размерность массивов - два элемента с индексами: 0, 1.
Итерации в цикле идут по n от 0 до 1 (два раза - по разу на каждое уравнение). В чём здесь ошибка в алгоритме? Честно говоря, пока не вижу ошибок, укажите пожалуйста :)

[JustAMan]

Тут ведь понимаете в чём проблема.. в случае с одним уравнением - с этими же индексами, оно работает верно, а вот в случае системы чего-то уже не срабатывает. (естественно, с разными функциями myF). Может в ней как раз что-то не так?

[Bulinator]

JustAMan
Если со скалярной функцией все работает, введите структуру

struct vectorvalue
{
    double x1;
    double x2;
};

перегрузите для нее операторы +/-/* и просто подставьте везде вместо double свою структуру.

-- Пт янв 14, 2011 01:45:02 --

Вообще, отличная задача:
написать программу, которая будет считывать систему уравнений(решаемых методом Рунге-Кутты), и решать ее. Притом можно написать ее для любого(вменяемого) числа уравнений.

[ewert]

JustAMan в сообщении #399584 писал(а):

в случае с одним уравнением - с этими же индексами, оно работает верно,

Она не может верно работать -- переменная k4 неверно вычисляется. Т.е. она должна давать правдоподобные результаты, но точность будет сильно занижена.

В остальном код похож на правду (я в тот раз невнимательно вчитался -- уж очень нелепо выбраны обозначения x1, x2, да и x2_t тоже, и это сильно отвлекает). Но вот что помимо этого -- сильные подозрения вызывает код
функции myF. Нет,синтаксически он верен, и что-то на выходе выдаст. Только вот систем с такой правой частью обычно не бывает (т.е. их не приходится рассматривать). Вы явно закодировали совсем не то, что требовалось. Напишите здесь Вышу систему в обычной математической форме.

[Bulinator]

ewert в сообщении #399606 писал(а):

Только вот систем с такой правой частью обычно не бывает

Частица под действием силы пропорциональной $cos{x_1}$ и силы сопротивления, скажем, воздуха(пропорциональна скорости $x_2=\dot{x}_1$).

[ewert]

А, пардон увидел систему. Естественно, Ваше программное описание правой части системы -- и впрямь не имеет ничего общего с заданием. И основная причина тому -- действительно, именно безобразный выбор идентификаторов. Ваши переменные x1 и x2 не имеют ничего общего с тем, что под этим понималось в условии.

[Bulinator]

struct vectorvalue
{
    vectorvalue()
    {
        x1=x2=0;
    }
    vectorvalue operator +(vectorvalue a)
    {
        a.x1+=this->x1;
        a.x2+=this->x2;
        return a;
    }
    vectorvalue operator -(vectorvalue a)
    {
        a.x1=this->x1-a.x1;
        a.x2=this->x2-a.x2;
        return a;
    }
    void operator =(vectorvalue a)
    {
        this->x1=a.x1;
        this->x2=a.x2;
    }
    vectorvalue operator *(double a)
    {
        vectorvalue res;
        res.x1=a*this->x1;
        res.x2=a*this->x2;
        return res;
    }
    
    double x1;
    double x2;
};

[JustAMan]

Ну тут, кстати, x2_t[ n ] - совершенно не нужно было описывать массивом, просто промежуточное вычисление.

По поводу k4 - ну возможно у меня не совсем точные формулы были в той книге, из которой я брал. Ну да ладно, уточним их) Но тут что-то совершенно иное у меня считается.. тут одно уточнение явно не поможет исправлению ошибки))

Bulinator
Аа.. у нас же там где myF теперь должен быть вектор значений, а не такой возврат одного числа, как в моём случае сейчас.. Ага, тогда я примерно понял в чём проблема.

Но погодите, я тогда совершенно не понимаю порядка вычислений. Например в случае с k1:
$k1 = h * f(x,y)$
где f(x,y) у нас теперь вектор значений (к примеру): [x1, x2]. И как тут обсчитываются две этих k1? У нас ведь k1 тоже вектором чтоль получится? Первое значение вектора это первый k1, а второй - это второй k1? Вряд ли... у меня, кажется, так и сделано уже А как же иначе они обсчитываются?

-- Пт янв 14, 2011 02:13:25 --

Bulinator
Спасибо большое за код! Мне бы сперва порядок вычислений "на пальцах" понять..

[Joker_vD]

Тьфу, правильно нас на Pascal ABC заставляют это писать. И правильно, что дают методичку, в которой русскими буквами сказано, как должны именоваться процедуры и параметры :)

[Bulinator]

JustAMan
ловите. В коде сами разберетесь?
Формулы взяты из вики.

#include <stdlib.h>
#include <iostream>
#include <math.h>
#include <stdio.h>
using namespace std;
/*
 * 
 */
struct vectorvalue
{
    vectorvalue()
    {
        x1=x2=0;
    }
    vectorvalue(double xx1, double xx2)
    {
        x1=xx1;
        x2=xx2;
    }
    vectorvalue operator +(vectorvalue a)
    {
        a.x1+=this->x1;
        a.x2+=this->x2;
        return a;
    }
    vectorvalue operator -(vectorvalue a)
    {
        a.x1=this->x1-a.x1;
        a.x2=this->x2-a.x2;
        return a;
    }
    void operator =(vectorvalue a)
    {
        this->x1=a.x1;
        this->x2=a.x2;
    }
    vectorvalue operator *(double a)
    {
        vectorvalue res;
        res.x1=a*this->x1;
        res.x2=a*this->x2;
        return res;
    }
    
    double x1;
    double x2;
};
//Сама функция, вообще-то, должна еще зависеть и от параметра t. В нашем частном случае его нет.
vectorvalue function(vectorvalue Y)
{
     return vectorvalue(Y.x2,-0.2*Y.x2-10*sin(Y.x1));
}
void RungeKutta()
{
    double endpoint=2; //конец отрезка
    vectorvalue func; //параметры x_1,x_2

    vectorvalue k1;
    vectorvalue k2;
    vectorvalue k3;
    vectorvalue k4;

    double h=0.01;
    double t=0;

    func.x1=0.3;
    func.x2=0;
    printf( "0: x1=%f, x2=%f\n", func.x1, func.x2);

    for(int i=0;t<endpoint;t+=h)
    {
        i++;
        k1=function(func);
        k2=function(func+k1*(h/2));
        k3=function(func+k2*(h/2));
        k4=function(func+k3*h);
        func=func+(k1+k2*2+k3*2+k4)*(h/6);
         printf( "%i: x1=%f, x2=%f\n", i, func.x1, func.x2);
    }


}


int main(int argc, char** argv)
{    
    RungeKutta();
    return (EXIT_SUCCESS);
}

(Оффтоп)

Эх... и черт меня дернул в науку податься. Получал бы сейчас тысячи американских и жил бы себе, не заботясь о великом объединении