Численный метод runge kutta

Question

Решите дифференциальное уравнение. dy/dx = xz + 1, dz/dx = -xy при x = 0,3, учитывая, что y = 0 и z = 1, когда x = 0, используя метод рун-кутты четвертого порядка. Я не могу начать эту проблему. Пожалуйста помоги. Мне требуется математическое решение с использованием численного анализа, а не программирования.

Answer 1

От классических Numerical Recipes текста http://www.aip.de/groups/soe/local/numres/bookcpdf/c16-1.pdf:

четвертого порядка методом Рунге-Кутта. На каждом шаге производная оценивается четыре раза: один раз в начальной точке, дважды в пробных средах и один раз на конечной точке испытания. Из этих производных вычисляется значение конечной функции.

k1 = hf(xn, yn) 
k2 = hf(xn + h/2 , yn + k1/2) 
k3 = hf(xn + h/2 , yn + k2/2) 
k4 = hf(xn + h, yn + k3) 
yn+1 = yn + k1/6 + k2/3 + k3/3 + k4/6 + O(h^5) 

Пример в C:

#include “nrutil.h" 
float **y,*xx; 
void rkdumb(float vstart[], int nvar, float x1, float x2, int nstep, void (*derivs)(float, float [], float [])) 
{ 
    void rk4(float y[], float dydx[], int n, float x, float h, float yout[], void (*derivs)(float, float [], float [])); 
    int i,k; 
    float x,h; 
    float *v,*vout,*dv; 
    v=vector(1,nvar); 
    vout=vector(1,nvar); 
    dv=vector(1,nvar); 
    for (i=1;i<=nvar;i++) { 
     v[i]=vstart[i]; 
     y[i][1]=v[i]; 
    } 
    xx[1]=x1; 
    x=x1; 
    h=(x2-x1)/nstep; 
    for (k=1;k<=nstep;k++) { 
     (*derivs)(x,v,dv); 
     rk4(v,dv,nvar,x,h,vout,derivs); 
     if ((float)(x+h) == x) 
      nrerror("Step size too small in routine rkdumb”); 
     x += h; 
     xx[k+1]=x; 
     for (i=1;i<=nvar;i++) { 
      v[i]=vout[i]; y[i][k+1]=v[i]; 
     } 
    } 
    free_vector(dv,1,nvar); 
    free_vector(vout,1,nvar); 
    free_vector(v,1,nvar); 
} 

Это должно вам начать работу. Надеюсь, поможет.