为什么这个简单的 C 代码集成代码会失败？

Question

我一直在 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;
 }

Answer 1

你很快就会在 y2 中发生大量溢出，它会更快地雪崩到 NaNs

永远不要低估 printf 在数字代码中的用处，所以让我们在这里做：

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];
  printf("aa.mu1: %g, aa.w1: %g, aa.eps: %g, aa.mu2: %g, aa.w2: %g \n", aa.mu1,
     aa.w1, aa.eps, aa.mu2, aa.w2);

  printf("x1: %g, x2: %g, y1: %g, y2: %g\n", x1, x2, y1, y2);
  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;


  printf("y1: %g, y2: %g\n", y1, y2);
  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);
  printf("%g %g %g %g %g\n\n", dv[0], dv[1], dv[2], dv[3], y1);
  return;
}

显示：

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: 2, x2: 6, y1: 0, y2: 1
y1: -6, y2: -222
-6 -222 236 1.09489e+07 -6

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: 1.97, x2: 4.89, y1: 1.18, y2: 54745.3
y1: -9.5984, y2: -1.46559e+10
-9.5984 -1.46559e+10 913.916 3.148e+30 -9.5984

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: 1.95201, x2: -7.32794e+07, y1: 4.56958, y2: 1.574e+28
y1: -50.8154, y2: 1.81548e+64
-50.8154 1.81548e+64 131360 -5.98381e+192 -50.8154

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: 1.49185, x2: 1.81548e+62, y1: 1313.6, y2: -5.98381e+190
y1: -2.57558e+06, y2: -inf
-2.57558e+06 -inf -nan -nan -2.57558e+06

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: -4290.84, x2: -inf, y1: -nan, y2: -nan
y1: -nan, y2: -nan
-nan -nan -nan -nan -nan

aa.mu1: 1, aa.w1: 1, aa.eps: 0, aa.mu2: 1, aa.w2: 6 
x1: -nan, x2: -nan, y1: -nan, y2: -nan
y1: -nan, y2: -nan
-nan -nan -nan -nan -nan

为什么这个简单的 C 代码集成代码会失败？

Why is this simple C-code integration code failing?

c

integration

physics

runge-kutta