Почему этот простой код интеграции кода C не работает?

Question

Я использовал этот код в C для интеграции с различными уравнениями, но сегодня я изменил его, чтобы интегрировать показанный, а файл .dat дает мне все столбцы, полные «-nan». Это проблема неправильного кодирования или ее просто то, что эти уравнения не должны решаться этой процедурой?

Вот как процедура интеграции, так и основной код. Благодаря

ПЕРВОГО КОДА (я нарезанное это объяснить, что делает каждую часть делать, я прошу прощения, если его неприятно в поле зрения)

#include <stdio.h>  
#include <math.h> 

#define a -1 

ПАРАМЕТРЫ

struct Par{ 
double mu1, mu2, w1, w2, eps; 
} aa; 

УРАВНЕНИЕ

void ecuaciones(int n, double v[], double dv[], double t){ 
double x1,y1,x2,y2; 

    x1=v[0]; 
    x2=v[1]; 
    y1=v[2]; 
    y2=v[3]; 
    y1=aa.mu1*x1-aa.w1*y1-(x1*x1+y1*y1)*x1+aa.eps*(x1-x2) ; 
    y2=aa.mu2*x2-aa.w2*y2-(x2*x2+y2*y2)*x2+aa.eps*(x2-x1) ; 
    dv[0]=y1; 
    dv[1]=y2; 
    dv[2]=aa.w1*x1+aa.mu1*y1-y1*(x1*x1+y1*y1)+aa.eps*(y1-y2) ; 
    dv[3]=aa.w2*x2+aa.mu2*y2-y2*(x2*x2+y2*y2)+aa.eps*(y2-y1) ; 

    return; 
} 

ГЛАВНАЯ

int main(){ 

    int i,j; 
    FILE *ptr; 
    double v[4],t,dt,t_pre,t_max; 

EXIT

ptr=fopen("NLIAC.dat","w"); 

    dt=0.01; 
    t_max=20; 

начальных условий

for (i = 2; i < 6; i++)  { 


    v[0]=i; 
    v[1]=6; 
    v[2]=0; 
    v[3]=1; 

ПАРАМЕТРЫ ОПРЕДЕЛЕНИЕ aa.w1 = 1; aa.mu1 = 1; aa.w2 = 6; aa.mu2 = 1; aa.eps = 0;

t=0.; 

ИНТЕГРАЦИЯ COMMAND

while(t<t_max){ 



rk4(ecuaciones,v,4,t,dt); 

EXPORT

fprintf(ptr,"%lg\t%lg\t%lg\t%lg\t%lg\n",t,v[0],v[1],v[2],v[3]); 

t+=dt; 
    }} 
     fprintf(ptr,"\n"); 

    fclose(ptr); 
    return(0); 
} 

И ВОТ ИНТЕГРАЦИЯ РЕГЛАМЕНТНОЕ (его штраф)

void rk4(void deri(int , double [], double [], double), \ 
double h[], int n, double t, double dt) 
{ 
#define naux 26 

int i; 
double k1[naux],k2[naux],k3[naux],k4[naux],h0[naux]; 
double dt2, dt6; 

dt2=dt/2.; 
dt6=dt/6.; 

for (i = 0 ; i<n; i++) 
h0[i] = h[i]; 

deri(n,h0,k1,t); 
for (i =0 ; i<n ; i++) 
    h0[i]=h[i]+dt2*k1[i]; 

deri(n,h0,k2,t+dt2); 
for (i =0 ; i<n ; i++) 
h0[i]=h[i]+dt2*k2[i]; 

deri(n,h0,k3,t+dt2); 
for (i =0 ; i<n ; i++) 
    h0[i]=h[i]+dt*k3[i]; 

deri(n,h0,k4,t+dt); 
for (i = 0 ; i<n ; i++) 
{h0[i]=h[i]+dt*k4[i];}; 

for (i =0; i<n ; i++) 
h[i]=h[i]+dt6*(2.*(k2[i]+k3[i])+k1[i]+k4[i]); 

return; 
} 

Answer 1

Вы получаете массивное переполнение в y2 довольно быстро, что лавины еще быстрее в NaN s

Никогда не undestimate полезности printf с в цифровом коде, так что давайте сделаем это здесь:

void ecuaciones(int n, double v[], double dv[], double t) 
{ 
    double x1, y1, x2, y2; 

    x1 = v[0]; 
    x2 = v[1]; 
    y1 = v[2]; 
    y2 = v[3]; 
    printf("aa.mu1: %g, aa.w1: %g, aa.eps: %g, aa.mu2: %g, aa.w2: %g \n", aa.mu1, 
    aa.w1, aa.eps, aa.mu2, aa.w2); 

    printf("x1: %g, x2: %g, y1: %g, y2: %g\n", x1, x2, y1, y2); 
    y1 = aa.mu1 * x1 - aa.w1 * y1 - (x1 * x1 + y1 * y1) * x1 + aa.eps * (x1 - x2); 
    y2 = aa.mu2 * x2 - aa.w2 * y2 - (x2 * x2 + y2 * y2) * x2 + aa.eps * (x2 - x1); 
    dv[0] = y1; 


    printf("y1: %g, y2: %g\n", y1, y2); 
    dv[1] = y2; 
    dv[2] = 
     aa.w1 * x1 + aa.mu1 * y1 - y1 * (x1 * x1 + y1 * y1) + aa.eps * (y1 - y2); 
    dv[3] = 
     aa.w2 * x2 + aa.mu2 * y2 - y2 * (x2 * x2 + y2 * y2) + aa.eps * (y2 - y1); 
    printf("%g %g %g %g %g\n\n", dv[0], dv[1], dv[2], dv[3], y1); 
    return; 
} 

шоу:

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: 2, x2: 6, y1: 0, y2: 1 
y1: -6, y2: -222 
-6 -222 236 1.09489e+07 -6 

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: 1.97, x2: 4.89, y1: 1.18, y2: 54745.3 
y1: -9.5984, y2: -1.46559e+10 
-9.5984 -1.46559e+10 913.916 3.148e+30 -9.5984 

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: 1.95201, x2: -7.32794e+07, y1: 4.56958, y2: 1.574e+28 
y1: -50.8154, y2: 1.81548e+64 
-50.8154 1.81548e+64 131360 -5.98381e+192 -50.8154 

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: 1.49185, x2: 1.81548e+62, y1: 1313.6, y2: -5.98381e+190 
y1: -2.57558e+06, y2: -inf 
-2.57558e+06 -inf -nan -nan -2.57558e+06 

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: -4290.84, x2: -inf, y1: -nan, y2: -nan 
y1: -nan, y2: -nan 
-nan -nan -nan -nan -nan 

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: -nan, x2: -nan, y1: -nan, y2: -nan 
y1: -nan, y2: -nan 
-nan -nan -nan -nan -nan