【发布时间】：2016-07-13 00:09:12
【问题描述】：
我一直在 C 中使用此代码与不同的方程集成，但今天我对其进行了修改以集成显示的那个，并且 .dat 文件为我提供了所有充满“-nan”的列。这是一个错误的编码问题还是只是这个方程不打算通过这个程序来解决？
这里是集成例程和主要代码。 谢谢
首先是代码（我砍了它解释每个部分的作用，如果它看起来不愉快，我很抱歉）
#include <stdio.h>
#include <math.h>
#define a -1
参数
struct Par{
double mu1, mu2, w1, w2, eps;
} aa;
方程式
void ecuaciones(int n, double v[], double dv[], double t){
double x1,y1,x2,y2;
x1=v[0];
x2=v[1];
y1=v[2];
y2=v[3];
y1=aa.mu1*x1-aa.w1*y1-(x1*x1+y1*y1)*x1+aa.eps*(x1-x2) ;
y2=aa.mu2*x2-aa.w2*y2-(x2*x2+y2*y2)*x2+aa.eps*(x2-x1) ;
dv[0]=y1;
dv[1]=y2;
dv[2]=aa.w1*x1+aa.mu1*y1-y1*(x1*x1+y1*y1)+aa.eps*(y1-y2) ;
dv[3]=aa.w2*x2+aa.mu2*y2-y2*(x2*x2+y2*y2)+aa.eps*(y2-y1) ;
return;
}
主要
int main(){
int i,j;
FILE *ptr;
double v[4],t,dt,t_pre,t_max;
退出
ptr=fopen("NLIAC.dat","w");
dt=0.01;
t_max=20;
初始条件
for (i = 2; i < 6; i++) {
v[0]=i;
v[1]=6;
v[2]=0;
v[3]=1;
参数定义 aa.w1=1; aa.mu1=1； aa.w2=6; aa.mu2=1; aa.eps=0;
t=0.;
集成命令
while(t<t_max){
rk4(ecuaciones,v,4,t,dt);
出口
fprintf(ptr,"%lg\t%lg\t%lg\t%lg\t%lg\n",t,v[0],v[1],v[2],v[3]);
t+=dt;
}}
fprintf(ptr,"\n");
fclose(ptr);
return(0);
}
这是集成例程（很好）
void rk4(void deri(int , double [], double [], double ), \
double h[], int n, double t, double dt)
{
#define naux 26
int i;
double k1[naux],k2[naux],k3[naux],k4[naux],h0[naux];
double dt2, dt6;
dt2=dt/2.;
dt6=dt/6.;
for (i = 0 ; i<n; i++)
h0[i] = h[i];
deri(n,h0,k1,t);
for (i =0 ; i<n ; i++)
h0[i]=h[i]+dt2*k1[i];
deri(n,h0,k2,t+dt2);
for (i =0 ; i<n ; i++)
h0[i]=h[i]+dt2*k2[i];
deri(n,h0,k3,t+dt2);
for (i =0 ; i<n ; i++)
h0[i]=h[i]+dt*k3[i];
deri(n,h0,k4,t+dt);
for (i = 0 ; i<n ; i++)
{h0[i]=h[i]+dt*k4[i];};
for (i =0; i<n ; i++)
h[i]=h[i]+dt6*(2.*(k2[i]+k3[i])+k1[i]+k4[i]);
return;
}
【问题讨论】：
-
我会检查 y1 和 y2 的重新定义。您在第二个公式中使用 y1，但是它的值已经改变。直接定义dv[0]和dv[1]，不用中间赋值。
-
谢谢，我会尝试并评论结果如何
-
我很高兴地宣布，LutzL 的建议确实是解决问题的方法。该线程可以关闭。谢谢你的帮助！！！
-
@GatoMazzei 你不能提前 2 分钟发帖吗？ ;-) 但很严重：如果 LutzL 不想发布答案，您可以自己回答您的问题。
标签： c integration physics runge-kutta