Runge Kutta in C для уравнения Лоренца

Question

Я пытаюсь вычислить систему Лоренца с помощью метода Рунге Кутты, но я не могу найти, где мой код имеет ошибку. Когда я запускаю его, система переходит в статическую точку, и я должен получить бабочку (аттрактор Лоренца). Я думаю, что это что-то вроде «for'loops» в разделе метода Runge Kutta. Вот мой код

#include<stdio.h> 
#include<stdlib.h> 
#include<math.h> 
double f(int,double,double []); 
double s,r,b; 
FILE *output; 
int main() 
{ 
    output=fopen("lorenzdata.dat","w"); 
    int j,N=3; 
    double x[2],dt,y[2],K1[2],K2[2],K3[2],K4[2],ti,t,i; 
    printf("\n Introduce parameters s, r and b sepparated by a space:"); 
    scanf("%lf %lf %lf",&s,&r,&b); 
    printf("\n Introduce initial conditions t0,x0,y0 and z0:"); 
    scanf("%lf %lf %lf %lf",&ti,&x[0],&x[1],&x[2]); 
    printf("\n Introduce time step:"); 
    scanf("%lf",&dt); 
    printf("\n Specify time to compute in the equations:"); 
    scanf("%lf",&t); 
/* Loop para Runge Kutta 4 */ 
    do 
    { 
     /* printf("%9.4f %9.4f %9.4f %9.4f \n",ti,x[0],x[1],x[2]); */ 
     i++; 
     for(j=0;j<N;j++) 
     { 
      K1[j] = f(j,ti,x); 
     } 
     for(j=0;j<N;j++) 
     { 
      y[j] = x[j]+(K1[j]/2)*dt; 
     } 
     for(j=0;j<N;j++) 
     { 
      K2[j] = f(j,ti+dt/2,y); 
     } 
     for(j=0;j<N;j++) 
     { 
      y[j] = x[j]+(K2[j]/2)*dt; 
     } 
     for(j=0;j<N;j++) 
     { 
      K3[j] = f(j,ti+dt/2,y); 
     } 
     for(j=0;j<N;j++) 
     { 
      y[j] = x[j]+(K3[j])*dt; 
     } 
     for(j=0;j<N;j++) 
     { 
      K4[j] = f(j,ti+dt,y); 
     } 
     for(j=0;j<N;j++) 
     { 
      x[j] += dt*((K1[j]/6)+(K2[j]/3)+(K3[j]/3)+(K4[j]/6)); 
     } 
     ti +=dt; 
     fprintf(output,"%9.4f %9.4f %9.4f \n",x[0],x[1],x[2]); 
    }while(i*dt <= t); 
    fclose(output); 
    return 0; 
} 
/* Definimos la funcion */ 
double f(int m,double h,double x[]) 
{ 
    if(m==0) 
    { 
     return s*(x[1]-x[0]); 
    } 
    else if(m==1) 
    { 
     return x[0]*(r-x[2])-x[1]; 
    } 
    else if(m==2) 
    { 
     return x[0]*x[1]-b*x[2]; 
    } 
} 

Заранее спасибо

EDIT

Я написал код еще раз (в более упрощенном виде) на Linux, и она работает нормально, это, кажется, мой редактор в Windows, было что-то (возможно, кодирование), что когда я скомпилировал код, он выбрасывал много бесконечных значений. Не знаю, почему, и мне потребовалось много времени, чтобы заметить это. Спасибо за вашу помощь

Answer 1

линии, которые выглядят как

for(j=0;j<N;j++) 
{ 
    y[j] = x[j]+(K1[j]/2)*dt; 
} 

неверны.

Это должно выглядеть,

for(j=0;j<N;j++) 
{ 
    K1[j] = f(j,ti,x[j],y[j]); 
    L1[y] = g(j,ti,x[j],y[j]); 
} 
for(j=0;j<N;j++) 
{ 
    K2[j] = f(j,ti+dt/2,x+k1[j]/2,y+L1[j]/2); 
    L2[j] = g(j,ti+dt/2,x+k1[j]/2,y+L1[j]/2); 
} 
for(j=0;j<N;j++) 
{ 
    K3[j] = f(j,ti+dt/2,x+K2[j]/2,y+L2[j]/2); 
    L3[j] = g(j,ti+dt/2,x+K2[j]/2,y+L2[j]/2); 
} 
for(j=0;j<N;j++) 
{ 
    K4[j] = f(j,ti+dt,x+K3[j],y+L3[j]); 
    L4[j] = g(j,ti+dt,x+K3[j],y+L3[j]); 
} 
for(j=0;j<N;j++) 
{ 
    x[j] += dt*((K1[j]/6)+(K2[j]/3)+(K3[j]/3)+(K4[j]/6)); 
    y[j] += dt*((L1[j]/6)+(L2[j]/3)+(L3[j]/3)+(L4[j]/6)); 
} 

Рунге-Кутта как это работает.

У вас есть система уравнений.

дх/дт = е (т, х, у) ду/дт = г (т, х, у)

Для каждого аргумента функции, необходимо вспомогательный этап Рунге-Кутта (k1, k2 и т. д.). Итак, для каждого подшага «k» вам нужны «подшаговые» обновленные значения x и y. В качестве примера для второго подшагания x + k1/2 и y + l1/2.

Answer 2

Хорошо игнорируя любую другую проблему, где переменная «i» инициализирована?

Answer 3

Я не думаю, что вы получили ничего плохого кроме того, что N = 3 означает, что распределение х, у, и Ks должны быть

х [3], у [3], K1 [3], .. .

вместо double x [2], dt, y [2], K1 [2], K2 [2], K3 [2], K4 [2], ti, t, i;