In this notebook we will investigate different methods of numerical integration and see their stability. The Particle class by default uses the velocity Verlet method, a symplectic method (meaning it nearly conserves energy) with error $O(\Delta t^2)$. The algorithm is as follows: $$$$ $$ x_{n+1} = x_n + v_n \Delta t + \frac{1}{2} a_n \Delta t^2 $$ $$$$ $$ v_{n+1} = v_n + \frac{a_n+ a_{n+1}}{2} \Delta t $$ $$$$ We will subclass the Particle class and create an EulerParticle class which implements the Euler method, and see how it compares with the velocity Verlet method. First, we import the necessary libraries:
#NAME: Integration Methods
#DESCRIPTION: Investigation of different numerical integration techniques and their stability in dynamics simulations.
from pycav.mechanics import *
from vpython import *
import numpy as np
Next, we will create the EulerParticle class, which implements the Euler method: $$$$ $$ x_{n+1} = x_n + v_n \Delta t $$ $$$$ $$ v_{n+1} = v_n + a \Delta t $$ $$$$
class EulerParticle(Particle):
def update(self, dt):
self.pos += (self.v * dt)
self.v += (dt * self.total_force) * self.inv_mass
We will also subclass the System class to create a MeasuredSystem class which will be able to tell us interesting properties of a system.
class MeasuredSystem(System):
@property
def kinetic_energy(self):
_ke = 0.
for particle in self.particles:
_ke += (0.5 * np.inner(particle.v, particle.v))/(particle.inv_mass)
return _ke
@property
def gravitational_potential_energy(self):
_gpe = 0.
for particle1 in self.particles:
for particle2 in self.particles:
dist = particle1.pos - particle2.pos
if particle1 is not particle2:
r = np.sqrt(np.inner(dist, dist))
_gpe = -1. / (r * particle1.inv_mass * particle2.inv_mass)
_gpe /= 2.
return _gpe
@property
def total_energy(self):
return (self.gravitational_potential_energy + self.kinetic_energy)
We will then set up the same system using normal Particles and EulerParticles. The EulerPacticle system will be colored white so that we can see the difference.
planet_verlet = Particle(pos = np.array([200.,0.,0.]),
v = np.array([0., 31.62, 0.]),
inv_mass = 1.,
make_trail = True,
radius = 10.)
star_verlet = Particle(pos = np.array([0., 0., 0.]),
v = np.array([0., 0., 0.]),
inv_mass = 1./200000.,
radius = 20.)
planet_euler = EulerParticle(pos = np.array([200.,0.,0.]),
v = np.array([0., 31.62, 0.]),
inv_mass = 1.,
make_trail = True,
color = [1., 1., 1.],
radius = 10.)
star_euler = EulerParticle(pos = np.array([0., 0., 0.]),
v = np.array([0., 0., 0.]),
inv_mass = 1./200000.,
radius = 20.)
verlet = [planet_verlet, star_verlet]
euler = [planet_euler, star_euler]
We will now simulate these two systems side by side, on the same canvas. We will also plot the number of steps taken and the Star-Planet distance for each system. You should see that for any time step, the Verlet system stays stable while the Euler method has the orbits spiralling out. Try looking at other properties of the system (total energy, kinetic energy, etc.) and see how they are preserved by these methods.
scene1 = canvas(title='Comparison of integrating methods')
verlet_sys = MeasuredSystem(collides = False, interacts = True, visualize = True, particles = verlet, canvas = scene1)
euler_sys = MeasuredSystem(collides = False, interacts = True, visualize = True, particles = euler, canvas = scene1)
graph1 = graph(x=0, y=0,
xtitle='Steps', ytitle='Star-Planet distance',
foreground=color.black, background=color.black)
f1 = gcurve(color=color.white)
f2 = gcurve(color=color.red)
dt = 0.1
while True:
rate(150)
verlet_sys.simulate(dt)
euler_sys.simulate(dt)
dis_verlet = planet_verlet.pos - star_verlet.pos
dis_euler = planet_euler.pos - star_euler.pos
f2.plot(verlet_sys.steps, np.sqrt(np.inner(dis_verlet, dis_verlet)))
f1.plot(verlet_sys.steps, np.sqrt(np.inner(dis_euler, dis_euler)))
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) <ipython-input-5-b3f4078dfd5d> in <module>() 11 12 while True: ---> 13 rate(150) 14 verlet_sys.simulate(dt) 15 euler_sys.simulate(dt) /Applications/anaconda/lib/python3.5/site-packages/vpython/vpython.py in __call__(self, maxRate) 132 if (self.rval != maxRate) and (maxRate >= 1.0): 133 self.rval = maxRate --> 134 super(RateKeeper2, self).__call__(maxRate) 135 136 if sys.version > '3': /Applications/anaconda/lib/python3.5/site-packages/vpython/rate_control.py in __call__(self, maxRate) 176 if self.calls > 0: # if there were some calls to rate after the last render 177 dt = self.lastSleep + self.calls*(self.userTime + self.callTime + self.delay) - _clock() --> 178 _sleep(dt) 179 self.buildStrategy(maxRate) 180 nr = self.whenToRender[0] /Applications/anaconda/lib/python3.5/site-packages/vpython/rate_control.py in _sleep(dt) 47 dtsleep = nticks*_tick 48 t = _clock() ---> 49 time.sleep(dtsleep) 50 t = _clock()-t 51 dt -= t KeyboardInterrupt: