Runge-Kutta法模拟Lorenz吸引子异常：变量演化后趋稳停滞求助

我尝试用Runge-Kutta4方法模拟Lorenz吸引子，但未得到预期结果：系统先演化一段时间，随后所有变量变为常数，演化停滞。我模拟了多个粒子，初始条件通过伪随机生成器在0到1之间生成，以下是编写的代码，恳请帮忙排查问题。

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

#define N 1000000
#define Nt 3
#define fran rand()/((double)RAND_MAX+1)

double h,s,r,b;

void rk4(double *x, double *y, double *z);
double fx(double x,double y,double z);
double fy(double x,double y,double z);
double fz(double x,double y,double z);


int main()
{
    double x[Nt], y[Nt], z[Nt];
    int i,j;
    FILE *f;
    char name[20]="Lorenz.txt";
    srand(time(NULL));

    f = fopen(name,"w");
    h = 0.01;
    s = 10.0;
    r = 18.0;
    b = 8.0/3.0;

    for(i=0;i<Nt;i++)
    {
        x[i] = fran;
        y[i] = fran;
        z[i] = fran;
    }

    if (f==NULL) return 1;

    for(i=0;i<N;i++)
    {
        rk4(x,y,z);

        for(j=0;j<Nt;j++)
        {
            fprintf(f,"%lf\t%lf\t%lf\t%lf\t",h*i,x[j],y[j],z[j]);
        }

        fprintf(f,"\n");
    }

    fclose(f);

}

double fx(double x,double y,double z)
{
    return s * (y - x);
}

double fy(double x,double y,double z)
{
    return x * (r - z) - y;
}

double fz(double x,double y,double z)
{
    return x * y - b * z;
}


void rk4(double *x, double *y, double *z)
{
    int i;
    double xn,yn,zn;
    double k1[3], k2[3], k3[3], k4[3];

    for(i=0;i<Nt;i++)
    {

        xn = x[i];
        yn = y[i];
        zn = z[i];

        k1[0] = h * fx(xn,yn,zn);
        k1[1] = h * fy(xn,yn,zn);
        k1[2] = h * fz(xn,yn,zn);
        k2[0] = h * fx(xn+0.5*k1[0],yn+0.5*k1[1],zn+0.5*k1[2]);
        k2[1] = h * fy(xn+0.5*k1[0],yn+0.5*k1[1],zn+k1[2]*0.5);
        k2[2] = h * fz(xn+0.5*k1[0],yn+0.5*k1[1],zn+k1[2]*0.5);
        k3[0] = h * fx(xn+0.5*k2[0],yn+0.5*k2[1],zn+k2[2]*0.5);
        k3[1] = h * fy(xn+0.5*k2[0],yn+0.5*k2[1],zn+k2[2]*0.5);
        k3[2] = h * fz(xn+0.5*k2[0],yn+0.5*k2[1],zn+k2[2]*0.5);
        k4[0] = h * fx(xn+k3[0],yn+k3[1],zn+k3[2]);
        k4[1] = h * fy(xn+k3[0],yn+k3[1],zn+k3[2]);
        k4[2] = h * fz(xn+k3[0],yn+k3[1],zn+k3[2]);

        x[i] = xn + (k1[0]+2*k2[0]+2*k3[0]+k4[0])/6.0;
        y[i] = yn + (k1[1]+2*k2[1]+2*k3[1]+k4[1])/6.0;
        z[i] = zn + (k1[2]+2*k2[2]+2*k3[2]+k4[2])/6.0;

    }
}

问题排查与解决 #

核心原因：参数选择导致系统收敛到不动点 #

Lorenz系统的行为由参数s、r、b决定：

• 当r < 24.74时，系统会收敛到稳定的不动点，不会出现混沌的吸引子行为。你当前设置的r=18.0正处于这个区间，所以粒子演化一段时间后会停在不动点，变量不再变化。

修复方案 #

将参数r修改为经典的混沌阈值以上的值，比如r=28.0（这是Lorenz吸引子的标准参数之一）：

r = 28.0;

修改后重新运行代码，就能观察到粒子在混沌吸引子上的演化轨迹，不会再出现停滞现象。

代码其他细节验证 #

• 随机数生成：fran的定义rand()/((double)RAND_MAX+1)正确，能生成0到1之间的均匀分布随机数。
• RK4实现：四阶龙格-库塔的k1~k4计算步骤、状态更新公式均符合标准，代码逻辑无错误。
• 文件输出：每次迭代输出所有粒子的时间步和状态，格式正常。

内容的提问来源于stack exchange，提问作者Ander Gabarrus