龙格库塔在C为洛伦兹方程

我想计算洛伦茨系统使用龙格库塔方法，但我找不到我的代码有错误的地方。当我运行它时，系统进入静态点，我应该获得一个蝴蝶（洛伦兹吸引子）。我认为这是Runge Kutta方法部分的'for'循环中的一些东西。这里是我的代码龙格库塔在C为洛伦兹方程

#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上编辑了一些东西（也许是编码），当我编译代码时，它抛出了很多infinites值。不知道为什么，我花了很多时间才注意到这一点。感谢您的帮助

来源

2012-02-18 David Winchester

双h作为函数f的输入点是什么？它似乎没有被使用。有必要吗？ – 2012-02-20 00:28:40

mm nop，我只是为了做一个更一般化的代码，以防万一diff方程式明确地取决于t – 2012-02-20 03:05:04

，看起来像

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)); 
}

龙格 - 库塔是这样工作的。你有一个方程组。

DX/DT = F（T，X，Y） DY/DT = G（T，X，Y）

有关的功能的每个参数，需要一个龙格 - 库塔子步骤（在k1，k2等）。因此，对于每个子步骤“k”，您需要x和y的“子步骤”更新值。以第二个子步骤x + k1/2和y + l1/2为例。

来源

2012-02-18 06:06:52 user1139069

嗯......但我这样做考虑到函数f的j组件。如果你在代码的最后看到，f有3个组件。对x也是如此，考虑到我在方程（x [0]，x [1]，x [2]）上有x，y，z，我有x [j]。如果我像你说的那样编写代码，'for'循环不应该在那里我认为 – 2012-02-18 14:52:01

你写这些代码的方式非常混乱，但我明白你现在正在做什么。问题是你乘以dt两倍。在[y] = x [j] +（K3 [j]）* dt这一行中''你乘以'dt'，所以在这一行中'x [j] + = dt *（（K1 [j]/6）+（K2 [j]/3）+（K3 [j]/3）+（K4 [j]/6））;''你不需要乘以'dt'。 – user1139069 2012-02-19 17:17:12

嗯我想我需要这样做，因为y [j]只是'更新'x []必须在f（） – 2012-02-20 01:53:06

那么忽略任何其他问题，变量'我'初始化？

来源

2012-02-18 05:51:59 Jxj

这甚至不应该编译; 'double'没有'++'运算符 – 2012-02-18 06:10:55

@SethCarnegie我有（至少）两个支持它的编译器。然而，它有点奇怪... – justin 2012-02-18 06:25:03

哦，我没有注意到，但是，无论如何，代码仍然给我错误的输出 – 2012-02-18 14:53:15

我不认为你有什么不妥，只是N = 3是指X，Y的分配，和KS应

X [3]，Y [3]，K1 [3]，.. 。

代替双X [2]，DT，Y [2]，K1 [2]，K2 [2]，K3 [2]，K 4 [2]，TI，T，I;

来源

2013-07-23 17:43:02 GaryX

龙格库塔在C为洛伦兹方程

回答

相关问题