spenam spenam - 11 days ago 8
C++ Question

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

My first post on Stack Overflow, be gentle. I wrote a code to follow the position on the x,y plane of a particle of mass M on a potential V(r) described by a four-dimensional system of equations of motion

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


Which are solved by using the Runge Kutta 4th Order method, where r=(x,y) and v=(vx,vy); now the state of the particle is defined by x, y and the angle theta between the vector v and the positive x-axis where the magnitude of the velocity is given by

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


where E is the energy in the plane and the potential V(r) is given by

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


now here is the code I made for the initial values

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


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

}


The problem here is that when I run the code with the initial conditions given, the values do not match with what it is supposed to according to the case given by the problem

what the graphic should look like with the initial conditions given
what the graphic should look like with the initial conditions given

now, I think i got right the implementation of the method but i do not know why the graphs do not match because when i run the code the particle goes away from the potential.

Any help will be appreciated.

Answer

The paths look chaotic with sharp turns. This requires an adaptive step size, you will need to implement some step size control. Either by comparing each step with two steps of half the step length or by using a method with embedded methods of higher order like Fehlberg or Dormand-Price.


More immediate errors:

  • define Em as V(1,1) to avoid unnecessary magic numbers
  • your initial position is, if you read the chart right,

    y0 = 3.0;
    y1 = -0.3905+k*0.0010; 
    

    with k=-1,0,1, note the minus sign.

  • your initial velocity is horizontal, and the kinetic energy is computed to complement the potential energy at that position. Thus

    y2 = v0 = -(std::sqrt((2.0*(E-V(y0, y1)))/M));
    y3 = v1 = 0.0;
    

With these changes and an adaptive solver I get (in python) the plot

reproduction of plot from cited article

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()