2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


В этом разделе нельзя создавать новые темы.



Начать новую тему Ответить на тему На страницу 1, 2  След.
 
 Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение13.01.2011, 05:30 


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

Мне требуется посчитать значения этой системы уравнений:
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 ];
         
      }
   }
   
}

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение13.01.2011, 09:03 
Заслуженный участник


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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение13.01.2011, 11:33 
Заслуженный участник


11/05/08
32166
JustAMan в сообщении #399143 писал(а):
Ошибка именно в математике. Сомневаюсь, что в алгоритме :)

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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение13.01.2011, 18:03 


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

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

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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение13.01.2011, 18:26 
Заслуженный участник


26/07/09
1559
Алматы
Хотя-бы разберитесь с вашими индексами и размерами массивов...

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение13.01.2011, 19:50 


06/10/10
106
Circiter в сообщении #399436 писал(а):
Хотя-бы разберитесь с вашими индексами и размерами массивов...

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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение13.01.2011, 23:27 


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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение13.01.2011, 23:40 
Заслуженный участник
Аватара пользователя


30/10/10
1481
Ереван(3-й участок)
JustAMan
Если со скалярной функцией все работает, введите структуру
Код:
struct vectorvalue
{
    double x1;
    double x2;
};

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

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

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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение14.01.2011, 00:04 
Заслуженный участник


11/05/08
32166
JustAMan в сообщении #399584 писал(а):
в случае с одним уравнением - с этими же индексами, оно работает верно,

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

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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение14.01.2011, 00:10 
Заслуженный участник
Аватара пользователя


30/10/10
1481
Ереван(3-й участок)
ewert в сообщении #399606 писал(а):
Только вот систем с такой правой частью обычно не бывает

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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение14.01.2011, 00:12 
Заслуженный участник


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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение14.01.2011, 00:49 
Заслуженный участник
Аватара пользователя


30/10/10
1481
Ереван(3-й участок)
Код:
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;
};

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение14.01.2011, 01:12 


06/10/10
106
Ну тут, кстати, 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
Спасибо большое за код! Мне бы сперва порядок вычислений "на пальцах" понять.. :-)

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение14.01.2011, 01:21 
Заслуженный участник


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

 Профиль  
                  
 
 Re: Метод Рунге-Кутта - не получается реализовать для системы.
Сообщение14.01.2011, 01:27 
Заслуженный участник
Аватара пользователя


30/10/10
1481
Ереван(3-й участок)
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);
}


(Оффтоп)

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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 30 ]  На страницу 1, 2  След.

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group