Метод Рунге-Кутты

[smbd2011]

Здравствуйте!

Мне нужно решить систему обыкновенных дифференциальных уравнений с использованием метода рунге-кутты. Фото системы прикреплено в картинке. № 2.2

Код написан согласно всем инструкциям учебника, однако практически сразу в результат выдается "1.#QNAN0". Что это значит? На 0 нигде не делю... Все честно вроде бы.

#include "stdio.h"
#include "stdlib.h"
double func1 (double x,double y, double t);
double func2 (double x,double y, double t);
double func3 (double x, double v, double t);
double v_ (double v, double x, double y, double t);
double w_ (double w, double x, double y, double t);
double x_ (double f, double x, double y, double t, double v);
double y_ (double f, double x, double y, double t, double w);
int main (void)
{
FILE *g;
int i = 0;
double *x, *y, *w, *v;
double alpha = 0.1, t = 0;
double k1, k2, k3, k4, k5, k6, tmp;
g = fopen("res.dat","w+");
x = (double*)malloc(sizeof(double)*1000 );
y = (double*)malloc(sizeof(double)*1000 );
v = (double*)malloc(sizeof(double)*1000 );
w = (double*)malloc(sizeof(double)*1000 );
x[0] = alpha;
y[0] = alpha*alpha-alpha*alpha*alpha*alpha;
v[0] = w[0] = 0;
while (t<=50)
{
t = t + 0.1;
v[i+1] = v_(v[i], x[i], y[i], t);
fprintf (g,"%lf\n", v[i+1]);
w[i+1] = w_(w[i], x[i], y[i], t);
x[i+1] = x_(x[i], x[i-1], y[i-1], t, v[i-1]);
y[i+1] = y_(y[i], x[i-1], y[i-1], t, w[i-1]);
//x[i+1] = v[i]*t + x[i];
//y[i+1] = w[i]*t + y[i];
fprintf(g," %lf %lf %lf\n", x[i], y[i], t);
i++;

}
free(x);
free(y);
free(v);
free(w);
}

double func1 (double x,double y,double t)
{
double res;
res = 2*x*(2*x*x - 1)*(1+x*x-x*x*x*x - y);
return res;
}

double func2 (double x,double y, double t)
{
double res;
res = (1 + x*x - x*x*x*x -y);
return res;
}

double func3 (double x, double v, double t)
{
double res;
res = v*t + x;
return res;
}

double v_ (double v, double x, double y, double t)
{
double tmp, res, k1, k2, k3, k4, k5, k6;
tmp = func1(x,y, t);
k1 = t*tmp;
tmp = func1(x + t*1/5,y + k1*1/5, t);
k2 = t*tmp;
tmp = func1(x + t*3/10,y + k1*3/40 + k2*9/40, t);
k3 = t*tmp;
tmp = func1(x + t*4/5,y + k1*44/45 - k2*56/15 + k3*32/9, t);
k4 = t*tmp;
tmp = func1(x + t*8/9,y + k1*19372/6561 - k2*25360/2187 + k3*64448/6561 - k4*212/729, t);
k5 = t*tmp;
tmp = func1(x + t,y + k1*9017/3168 - k2*355/33 + k3*46732/5247 + k4*49/176 - k5*5103/18656, t);
k6 = t*tmp;
res = v + k1*35/384 + k3*500/1113 + k4*125/192 - k5*2187/6784 + k6*11/84;
return res;
}

double w_ (double w, double x, double y, double t)
{
double tmp, res, k1, k2, k3, k4, k5, k6;
tmp = func2(x,y,t);
k1 = t*tmp;
tmp = func2(x + t*1/5,y + k1*1/5, t);
k2 = t*tmp;
tmp = func2(x + t*3/10,y + k1*3/40 + k2*9/40, t);
k3 = t*tmp;
tmp = func2(x + t*4/5,y + k1*44/45 - k2*56/15 + k3*32/9, t);
k4 = t*tmp;
tmp = func2(x + t*8/9,y + k1*19372/6561 - k2*25360/2187 + k3*64448/6561 - k4*212/729, t);
k5 = t*tmp;
tmp = func2(x + t,y + k1*9017/3168 - k2*355/33 + k3*46732/5247 + k4*49/176 - k5*5103/18656, t);
k6 = t*tmp;
res = w + k1*35/384 + k3*500/1113 + k4*125/192 - k5*2187/6784 + k6*11/84;
}

double x_ (double f, double x, double y, double t, double v)
{
double tmp, res, k1, k2, k3, k4, k5, k6;
tmp = v_(v, x,y,t);
k1 = t*tmp;
tmp = v_(v, x + t*1/5,y + k1*1/5, t);
k2 = t*tmp;
tmp = v_(v,x + t*3/10,y + k1*3/40 + k2*9/40, t);
k3 = t*tmp;
tmp = v_(v,x + t*4/5,y + k1*44/45 - k2*56/15 + k3*32/9, t);
k4 = t*tmp;
tmp = v_(v,x + t*8/9,y + k1*19372/6561 - k2*25360/2187 + k3*64448/6561 - k4*212/729, t);
k5 = t*tmp;
tmp = v_(v,x + t,y + k1*9017/3168 - k2*355/33 + k3*46732/5247 + k4*49/176 - k5*5103/18656, t);
k6 = t*tmp;
res = f + k1*35/384 + k3*500/1113 + k4*125/192 - k5*2187/6784 + k6*11/84;
return res;
}

double y_ (double f, double x, double y, double t, double w)
{
double tmp, res, k1, k2, k3, k4, k5, k6;
tmp = w_(w, x,y,t);
k1 = t*tmp;

tmp = w_(w, x + t*1/5,y + k1*1/5, t);
k2 = t*tmp;
tmp = w_(w,x + t*3/10,y + k1*3/40 + k2*9/40, t);
k3 = t*tmp;
tmp = w_(w,x + t*4/5,y + k1*44/45 - k2*56/15 + k3*32/9, t);
k4 = t*tmp;
tmp = w_(w,x + t*8/9,y + k1*19372/6561 - k2*25360/2187 + k3*64448/6561 - k4*212/729, t);
k5 = t*tmp;
tmp = w_(w,x + t,y + k1*9017/3168 - k2*355/33 + k3*46732/5247 + k4*49/176 - k5*5103/18656, t);
k6 = t*tmp;
res = f + k1*35/384 + k3*500/1113 + k4*125/192 - k5*2187/6784 + k6*11/84;
return res;
}