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

[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;
}

[Сыроежка]

Почему вы считаете, что кто-то должен разбираться в вашем коде? Если у вас возникает ошибка. то локализуйте ее по крайней мере! И укажите то предложение программы, где возникает эта ошибка.

[smbd2011]

Я не считаю, что кто-то должен разбираться в моем коде. В теме о халяве написано, что никто просто так тут не поможет и нужно прикреплять код.

Вопрос, если Вы прочитали, заключен в том, почему возникает такой результат:
1.#QNAN0?

Если бы я знала, где возникает ошибка, я бы ее локализовала.

[Сыроежка]

Сообщение от smbd2011

Я не считаю, что кто-то должен разбираться в моем коде. В теме о халяве написано, что никто просто так тут не поможет и нужно прикреплять код.

Вопрос, если Вы прочитали, заключен в том, почему возникает такой результат:
1.#QNAN0?

Если бы я знала, где возникает ошибка, я бы ее локализовала.

Это, очевидно, ошибка. Невозможно представить число с плавающей запятой.

Надо собственные программы писать, а не воровать. Тогда будете знать, где возникает ошибка!

[smbd2011]

С чем может быть связана проблема с представлением числа? Все функции и вычисления крайне простые. Числа не уходят в бесконечность и не совсем маленькие получаются сначала. Числа, как числа...
Это МОЯ программа. Помогли бы лучше, чем хамить.

[smbd2011]

В общем, учебник подкачал с коэффициентами.
От этого начинаются скачки.Мб, кому это поможет.