Численное Интегрирование Системы Дифференциальных Уравнений Методом Ру

Тема в разделе "C/C++/C#", создана пользователем DenProx, 2 окт 2012.

  1. DenProx

    DenProx Гость

    Доброго времени суток. У меня есть вот такая задачка:

    Разработать функцию для численного интегрирования системы дифференциальных уравнений методом Рунге - Кутта. Прототип функции:

    Код (C++):
    void runge_k(void f(double *y, double *ys, double t), double *y, int n, duble tn, double tk, int m, double delt);
    где:

    f - функция вычисления правых частей системы дифференциальных уравнений;
    y - массив размера n значений зависимых переменных;
    ys - массив размера n значений производных;
    n - порядок системы дифференциальных уравнений;
    t - независимая переменная;
    tn - начальное значение интервала интегрирования;
    tk - конечное значение интервала интегрирования;
    m - начальное число разбиений отрезка интегрирования ;
    delt - шаг интегрирования.

    Шаг интегрирования для метода использовать 0,0001.

    Помогите решить, пожалуйста !

    У меня есть примерное решение на С++ (решено не до конца)

    Код (C++):
    void runge_k(
    void f(double *y, double *ys, double t),
    double *y, int n,
    duble tn, double tk,
    int m, double delt)
    {
    double k1[n],k2[n],k3[n],k4[n],Yh[n+1];
    double h=delt;
    double* ys=y+n+1;
    double X;
    for(X=tn;X<tk;X+=h){
    f(y,ys,X);

    for( int i=0;i<n;i++) k1[i]=h*y[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+0.5*k1[i];      
    f(Yh,Yh+n+1,X+h*0.5);

    for( int i=0;i<n;i++) k2[i]=h*Yh[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+0.5*k2[i];
    f(Yh,Yh+n+1,X+h*0.5);

    for( int i=0;i<n;i++) k3[i]=h*Yh[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+k3[i];
    f(Yh,Yh+n+1,X+h);

    for( int i=0;i<n;i++) k4[i]=h*Yh[i+1];

    for( int i=0;i<n;i++) y[i]+=1./6.*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);
    }
    h=tk-X;
    f(y,ys,X);

    for( int i=0;i<n;i++) k1[i]=h*y[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+0.5*k1[i];      
    f(Yh,Yh+n+1,X+h*0.5);

    for( int i=0;i<n;i++) k2[i]=h*Yh[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+0.5*k2[i];
    f(Yh,Yh+n+1,X+h*0.5);

    for( int i=0;i<n;i++) k3[i]=h*Yh[i+1];
    for( int i=0;i<n;i++) Yh[i]=y[i]+k3[i];
    f(Yh,Yh+n+1,X+h);

    for( int i=0;i<n;i++) k4[i]=h*Yh[i+1];

    for( int i=0;i<n;i++) y[i]+=1./6.*(k1[i]+2*k2[i]+2*k3[i]+k4[i]);
    }
    Буду очень признателен за любую помощь!
     

Поделиться этой страницей