混沌散射模拟的值与基本情况不匹配

Question

我在 Stack Overflow 上的第一个 post，要温柔。我编写了一个代码来跟踪质量为 M 的粒子在由四维运动方程组描述的势 V(r) 上的 x,y 平面上的位置

M(dv/dt)=-grad V(r),      dr/dt=v,

使用 Runge Kutta 四阶方法求解，其中 r=(x,y) 和 v=(vx,vy)；现在粒子的状态由 x、y 和矢量 v 与正 x 轴之间的角度 theta 定义，其中速度的大小由

给出

|v|=sqrt(2(E-V(r))/M)

其中 E 是平面中的能量，势能 V(r) 由

给出

V(r)=x^2y^2exp[-(x^2+y^2)],

现在这是我为初始值制作的代码

x(0)=3,
y(0)=0.3905,
vx(0)=0,
vy(0)=-sqrt(2*(E-V(x(0), y(0)))),

其中 E=0.260*(1/exp(2))

    // RK4
    #include <iostream> 
    #include <cmath> 

    // constant global variables
    const double M = 1.0;
    const double DeltaT = 1.0;

    // function declaration
    double f0(double t, double y0, double y1, double y2, double y3); // derivative of y0
    double f1(double t, double y0, double y1, double y2, double y3); // derivative of y1
    double f2(double t, double y0, double y1, double y2, double y3); // derivative of y2
    double f3(double t, double y0, double y1, double y2, double y3); // derivative of y3
    void rk4(double t, double h, double &y0, double  &y1, double &y2, double &y3); // method of runge kutta 4th order
    double f(double y0, double y1); //function to use

    int main(void) 
    {
      double y0, y1, y2, y3, time, E, Em;
      Em = (1.0/(exp(2.0)));
      E = 0.260*Em;
      y0 = 3.0; //x
      y1 = 0.3905; //y
      y2 = 0.0; //vx
      y3 = -(std::sqrt((2.0*(E-f(3.0, 0.0)))/M)); //vy
      for(time = 0.0; time <= 400.0; time += DeltaT) 
        {
          std::cout << time << "\t\t" << y0 << "\t\t" << y1 << "\t\t" << y2 << "\t\t" << y3 << std::endl;
          rk4(time, DeltaT, y0, y1, y2, y3);
        }


      return 0;
    }

    double f(double y0, double y1)
    {
      return y0*y0*y1*y1*(exp(-(y0*y0)-(y1*y1)));
    }

    double f0(double t, double y0, double y1, double y2, double y3)
    {
      return y2;
    }

    double f1(double t, double y0, double y1, double y2, double y3)
    {
      return y3;
    }


    double f2(double t, double y0, double y1, double y2, double y3)
    {
      return 2*y0*((y0*y0)-1)*(y1*y1)*(exp(-(y0*y0)-(y1*y1)))/M;
    }

    double f3(double t, double y0, double y1, double y2, double y3)
    {
      return 2*(y0*y0)*y1*((y1*y1)-1)*(exp(-(y0*y0)-(y1*y1)))/M;
    }


    void rk4(double t, double h, double &y0, double &y1, double &y2, double &y3) // method of runge kutta 4th order
    {
      double k10, k11, k12, k13, k20, k21, k22, k23, k30, k31, k32, k33, k40, k41, k42, k43;
      k10 = h*f0(t, y0, y1, y2, y3);
      k11 = h*f1(t, y0, y1, y2, y3);
      k12 = h*f2(t, y0, y1, y2, y3);
      k13 = h*f3(t, y0, y1, y2, y3);
      k20 = h*f0(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2);
      k21 = h*f1(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2);
      k22 = h*f2(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2);
      k23 = h*f3(t+h/2, y0 + k10/2, y1 + k11/2, y2 + k12/2, y3 + k13/2);
      k30 = h*f0(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2);
      k31 = h*f1(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2);
      k32 = h*f2(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2);
      k33 = h*f3(t+h/2, y0 + k20/2, y1 + k21/2, y2 + k22/2, y3 + k23/2);
      k40 = h*f0(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33);
      k41 = h*f1(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33);
      k42 = h*f2(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33);
      k43 = h*f3(t + h, y0 + k30, y1 + k31, y2 + k32, y3 + k33);


      y0 = y0 + (1.0/6.0)*(k10 + 2*k20 + 2*k30 + k40);
      y1 = y1 + (1.0/6.0)*(k11 + 2*k21 + 2*k31 + k41);
      y2 = y2 + (1.0/6.0)*(k12 + 2*k22 + 2*k32 + k42);
      y3 = y3 + (1.0/6.0)*(k13 + 2*k23 + 2*k33 + k43);

    }

这里的问题是，当我运行给出初始条件的代码时，根据问题给出的情况，值与应该的值不匹配

在给定的初始条件下图形应该是什么样子

现在，我认为我正确地实施了该方法，但我不知道为什么图形不匹配，因为当我运行代码时，粒子远离势能。

任何帮助将不胜感激。

Answer 1

路径看起来很混乱，急转弯。这需要自适应步长，您将需要实施一些步长控制。通过将每个步骤与步长一半的两个步骤进行比较，或者使用具有更高阶嵌入方法的方法，如 Fehlberg 或 Dormand-Price。

更多即时错误：

将Em定义为V(1,1)以避免不必要的幻数
你的初始位置是，如果你没看错图表，
```
y0 = 3.0;
y1 = -0.3905+k*0.0010; 
```
和k=-1,0,1，注意减号。
你的初始速度是水平的，计算动能是为了补充那个位置的势能。于是
```
y2 = v0 = -(std::sqrt((2.0*(E-V(y0, y1)))/M));
y3 = v1 = 0.0;
```

通过这些更改和自适应求解器，我得到了（在 python 中）图

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# capture the structure of the potential
f  = lambda     r :     r*np.exp(-r);
df = lambda     r : (1-r)*np.exp(-r);
V  = lambda y1,y2 : f(y1*y1)*f(y2*y2);

M= 1.0
Em = V(1.0,1.0);
E = 0.260*Em;

def prime(t,y):
    x1,x2,v1,v2 = y        
    dV_dx1 = 2*x1*df(x1*x1)*f(x2*x2);
    dV_dx2 = 2*x2*df(x2*x2)*f(x1*x1);
    return [ v1, v2, -dV_dx1/M, -dV_dx2/M ];

# prepare and draw the contour plot
X1,X0=np.ogrid[-4:3:100j,-4:3:100j]
plt.contour(X0.ravel(), X1.ravel(), V(X0,X1), Em*np.arange(0,1,0.1), colors='k', linewidths=0.3)
# display grid and fix the coordinate ranges
plt.grid();plt.autoscale(False)

for k in range(-1,1+1):

    x01 = 3.0; 
    x02 = b = -0.3905 + 0.0010*k; 
    v01 = -( ( E-V(x01,x02) )*2.0/M )**0.5; 
    v02 = 0.0; 
    print "initial position (%.4f, %.4f), initial velocity (%.4f, %.4f)" % ( x01, x02, v01, v02 )

    t0 = 0.0
    tf = 50.0
    tol = 1e-10
    y0 = [ x01, x02, v01, v02 ]
    t = np.linspace(t0,tf,501); y = odeint(lambda y,t: prime(t,y) , y0, t)

    plt.plot(y[:,0], y[:,1], label="b=%.4f" % b, linewidth=2)

plt.legend(loc='best')
plt.show()

混沌散射模拟的值与基本情况不匹配

Values for chaotic scattering simulation do not match with the base case

c++

physics

runge-kutta