我想计算洛伦茨系统使用龙格库塔方法,但我找不到我的代码有错误的地方。当我运行它时,系统进入静态点,我应该获得一个蝴蝶(洛伦兹吸引子)。我认为这是Runge Kutta方法部分的'for'循环中的一些东西。这里是我的代码龙格库塔在C为洛伦兹方程
#include
#include
#include
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
{
K1[j] = f(j,ti,x);
}
for(j=0;j
{
y[j] = x[j]+(K1[j]/2)*dt;
}
for(j=0;j
{
K2[j] = f(j,ti+dt/2,y);
}
for(j=0;j
{
y[j] = x[j]+(K2[j]/2)*dt;
}
for(j=0;j
{
K3[j] = f(j,ti+dt/2,y);
}
for(j=0;j
{
y[j] = x[j]+(K3[j])*dt;
}
for(j=0;j
{
K4[j] = f(j,ti+dt,y);
}
for(j=0;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];
}
}
在此先感谢
编辑
我在Linux上重新写的代码(以更简化的方式),并且它运行正常,看来我在Windows上编辑了一些东西(也许是编码),当我编译代码时,它抛出了很多infinites值。不知道为什么,我花了很多时间才注意到这一点。 感谢您的帮助
+0
双h作为函数f的输入点是什么?它似乎没有被使用。有必要吗? –
2012-02-20 00:28:40
+0
mm nop,我只是为了做一个更一般化的代码,以防万一diff方程式明确地取决于t –
2012-02-20 03:05:04