Runge Kutta in C for Lorenz equation
-
30-04-2021 - |
Question
i'm trying to compute the Lorenz system using Runge Kutta method, but i can't find where my code have an error. When I run it, the system goes to a static point and i should obtain a butterfly (the Lorenz attractor). I think it's something on the 'for'loops in the Runge Kutta method section. Here's my code
#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];
}
}
Thanks in advance
EDIT
I wrote the code again (in a more simplified way) on linux and it runs ok, it seems my editor on Windows had something (maybe the encoding) that when i compiled the code, it throwed a lot of infinites values. Don't know why and it took me a lot of time to notice that. Thanks for your help
Solution
The lines that look like
for(j=0;j<N;j++)
{
y[j] = x[j]+(K1[j]/2)*dt;
}
are wrong.
It should look like,
for(j=0;j<N;j++)
{
K1[j] = f(j,ti,x[j],y[j]);
L1[j] = 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));
}
Runge-Kutta works like this.
You have a system of equations.
dx/dt = f(t,x,y) dy/dt = g(t,x,y)
For each argument of the functions, you need a Runge-Kutta sub step (the k1, k2, etc...). So for each sub-step "k", you need the "sub-step" updated values of x and y. As an example for the 2nd sub-step x+k1/2 and y+l1/2.
OTHER TIPS
Well ignoring any other issue, where is the variable 'i' initialised?
I do not think you got anything wrong except that N=3 means the allocation of x, y, and Ks should be
x[3], y[3], K1[3], ...
instead of double x[2],dt,y[2],K1[2],K2[2],K3[2],K4[2],ti,t,i;