洛伦兹方程C中的龙格-库塔

Question

我试图使用Runge-Kutta方法计算Lorenz系统，但我找不到我的代码有错误的地方。当我运行它时，系统进入一个静态点，我应该得到一个蝴蝶（洛伦兹吸引子）。我认为这是龙格-库塔方法部分的"for"部分的内容。这是我的代码

#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];
    }
}

提前感谢

编辑

我在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[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));
}

龙格库塔是这样工作的。

你有一个方程组。

dx/dt=f（t，x，y）dy/dt=g（t，x，y）

对于函数的每个参数，您需要一个龙格-库塔子步骤（k1、k2等）。因此，对于每个子步骤"k"，您需要"子步骤"更新的x和y值。作为第二个子步骤x+k1/2和y+l1/2的示例。

Answer 2

忽略任何其他问题，变量'i'在哪里初始化？

Answer 3

我认为你没有做错什么，除了N=3意味着x、y和Ks的分配应该是

x[3]，y[3]，K1[3]。。。

而不是双x[2]，dt，y[2]，K1[2]，K2[2]，K3[2]，K4[2]，ti，t，i；