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

[DenProx]

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

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

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.

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

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

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]);
}

Буду очень признателен за любую помощь!

 

[DenProx]

Еще забыл написать, что необходимо применить эту функцию для интегрирования дифференциального уравнения 3-го порядка:

y''' + 2y'' + 3y' + y = 5 + x^2

в интервале x [0,2] с шагом delta X = 0,1
и начальными условиями: x=0; y(0) =1; y'(0) = 0,1 ; y''(0) = 0