Exploring the 3-body Problem (Part I, Jacobi coordinates and Symplectic Integration)¶
Joseph E. McEwen (c) 2017¶
Please report all errors, typos, etc.
What is this notebook?¶
- some fun with the 3-body problem, symplectic integrators, and learning Julia
- disclaimer: I am still learning Julia - please point out areas of improvement
- the contents are:
(code for video is below)
The 3-Body Hamiltonian¶
The Hamiltonian for the 3-body problem is \begin{align} H(p,q) = \frac{p_1^2}{2 m_1} + \frac{p_2^2}{2 m_2} + \frac{p_3^2}{2 m_3} - G \frac{m_1 m_2}{q_{12}} - G \frac{m_1 m_3}{q_{13}} - G \frac{m_2 m_3}{q_{23}}~, \end{align} where $q_{ij} = \big| \mathbf{q}_i - \mathbf{q}_j \big|$ (like wise for the vector $\mathbf{q}_{ij} = \mathbf{q}_i - \mathbf{q}_j$). The Hamiltonian can be split into kinetic energy and potential energy parts \begin{align} H(p,q) = K(p) + \Phi(q)~, \end{align} where $\Phi$ is the potential energy function. The force on the ith particle from the jth particle is given by the gradient of the potential function \begin{align} \begin{split} \mathbf{F}_{ij} &= - \boldsymbol{\nabla_{q_i}} \Phi(q_{ij}) \\ & =- \boldsymbol{\nabla_{q_i}} \left[ -G \frac{m_i m_j}{q_{ij}} \right] \\ & = -G \frac{m_i m_j}{q_{ij}^2} \boldsymbol{\nabla_{q_i}} \big| \mathbf{q}_i - \mathbf{q}_j \big| \\ & = -G \frac{m_i m_j}{q_{ij}^2} \frac{ \mathbf{q}_i - \mathbf{q}_j }{q_{ij} } \\ & = -G \frac{m_i m_j}{q_{ij}^2} \hat{q}_{ij} ~. \end{split} \end{align}
Jacobi coordinates¶
Jacobi coordinates are a coordinate transformation particulary suited for describing the dynamics of many body system.
- For two particles Jacobi coordinates are the vector from the 2nd particle to the first along with the center of mass vector.
- For three particles Jacobi coordinates can be: the vector from the 2nd particle to the first, the vector from the 3rd particle to the center of mass of the 2nd and 1st pair, and the center of mass vector.
- This procedure of clustering particles can be carried up to any number of particles.
- As long as the center of mass remains constant, it can be ignored in calculations.
For this notebook, the coordinate transformation between configuration coordinates $\mathbf{q}_i$ to Jacobi coordinates $\mathbf{s}_i$ is
\begin{align} \left( \begin{array}{c} \mathbf{s}_1 \\ \mathbf{s}_2 \\ \mathbf{R}_\text{cm} \end{array} \right) = \left( \begin{array}{ccc} 1 & 0 & -1 \\ - \frac{m_1}{m_1 + m_3} & 1 & -\frac{m_3}{m_1 + m_3} \\ \frac{m_1}{M} & \frac{m_2}{M} & \frac{m_2}{M} \end{array} \right) \left( \begin{array}{c} \mathbf{q}_1 \\ \mathbf{q}_2 \\ \mathbf{q}_3 \end{array} \right) ~, \end{align} where $M= m_1 + m_2 + m_3 $ and $R_\text{cm}$ is the center of mass coordinate.
The inverse transformation is given by \begin{align} \left( \begin{array}{c} \mathbf{q}_1 \\ \mathbf{q}_2 \\ \mathbf{q}_3 \end{array} \right)= \left( \begin{array}{ccc} \frac{m_3}{m_1+m_3} & -\frac{m_2}{M} & 1 \\ 0 & \frac{m_1 +m_3}{M} & 1 \\ \frac{-m_1}{m_1 +m_3} & \frac{-m_2}{M} & 1 \end{array} \right) \left( \begin{array}{c} \mathbf{s}_1 \\ \mathbf{s}_2 \\ \mathbf{R}_\text{cm} \end{array} \right)~. \end{align}
Ignoring the center of mass the coordinates and their time derivatves are \begin{align} \begin{split} & \mathbf{s}_1 = \mathbf{q}_1 - \mathbf{q}_3 \\ & \dot{\mathbf{s}}_1 = \dot{\mathbf{q}}_1 - \dot{\mathbf{q}}_3~, \end{split} \end{align}
\begin{align} \begin{split} & \mathbf{s}_2 = \frac{m_1\mathbf{q}_1 - m_3\mathbf{q}_3}{m_1 + m_3} + \mathbf{q}_2 \\ & \dot{\mathbf{s}}_2 = \frac{m_1 \dot{\mathbf{q}}_1 - m_3\dot{\mathbf{q}}_3}{m_1 + m_3} + \dot{\mathbf{q}}_2 ~. \end{split} \end{align}
It is useful to define the reduced masses \begin{align} \mu_1= \frac{ m_1 m_3}{m_1 + m_3} ~, \quad \mu_2= \frac{m_2(m_1 + m_3)}{M} ~. \end{align}
The kinetic energy in Jacobi coordinates is: \begin{align} \begin{split} \frac{1}{2} m_1 \dot{q}_1^2 & = \frac{1}{2} m_1 \left( \frac{m_3^2}{(m_1 + m_3)^2} \dot{s}_1^2 + \frac{m_2^2}{M^2} \dot{s}_2^2 + \dot{R}^2_\text{cm} - \frac{2 m_3 m_2}{M(m_1 + m_3)} \dot{\mathbf{s}}_1 \cdot \dot{\mathbf{s}}_2 + \frac{2 m_3}{m_1 + m_3} \dot{\mathbf{s}}_1 \cdot \dot{\mathbf{R}}_\text{cm} - \frac{2 m_2}{M} \dot{\mathbf{s}}_2 \cdot \dot{\mathbf{R}}_\text{cm} \right) ~, \end{split} \end{align}
\begin{align} \begin{split} \frac{1}{2} m_2 \dot{q}_2^2 & = \frac{1}{2} m_2 \left( \frac{ (m_1^2 + 2m_1 m_3 + m_3^2)}{M^2} \dot{s}_2^2 + \dot{R}_\text{cm}^2 + \frac{2(m_1 + m_3)}{M} \dot{\mathbf{s}}_2 \cdot \dot{\mathbf{R}}_\text{cm}\right) ~, \end{split} \end{align}
\begin{align} \begin{split} \frac{1}{2} m_3 \dot{q}_3^2 & = \frac{1}{2} m_3 \left(\frac{m_1^2}{(m_1 +m_3)^2}\dot{s}_1^2 + \frac{m_2^2}{M^2} \dot{s}_2^2 + \dot{R}^2_\text{cm} + \frac{2 m_1 m_2}{M(m_1 + m_3)}\dot{\mathbf{s}}_1 \cdot \dot{\mathbf{s}}_2 - \frac{2 m_1}{m_1 + m_3}\dot{\mathbf{s}}_1 \cdot \dot{\mathbf{R}}_\text{cm} - \frac{2 m_2}{M} \dot{\mathbf{s}}_2 \cdot \dot{\mathbf{R}}_\text{cm}\right) ~, \end{split} \end{align}
\begin{align} \begin{split} \frac{1}{2} \left( m_1 \dot{q}_1^2 + m_2 \dot{q}_2^2 + m_3 \dot{q}_3^2 \right) & = \frac{1}{2} \left( \frac{m_1m_3}{m_1 + m_3} \dot{s}_1^2 + \frac{m_2(m_1 + m_3)}{M} \dot{s}_2^2 + M \dot{R}^2_\text{cm} \right)~. \end{split} \end{align}
The differences in configuration vectors is \begin{align} \mathbf{q}_{21}= \mathbf{s}_2 - \frac{m_3}{m_1 + m_3} \mathbf{s}_1~, \end{align}
\begin{align} \mathbf{q}_{31} = - \mathbf{s}_1~, \end{align}
\begin{align} \mathbf{q}_{23} = \mathbf{s}_2 + \frac{m_1}{m_1 + m_3} \mathbf{s}_1 ~. \end{align}
The Hamiltonian in Jacobi coordinates is \begin{align} \begin{split} H = \frac{p_1^2}{2 \mu_1} + \frac{p^2_2}{2 \mu_2} - G m_1 m_2 \frac{1}{\Big| \mathbf{s}_2 - \frac{m_3}{m_1 + m_3} \mathbf{s}_1\Big|} - G m_1 m_3 \frac{1}{\Big| \mathbf{s}_1 \Big| } - G m_2 m_3 \frac{1}{\Big| \mathbf{s}_2 - \frac{m_1}{m_1 + m_3} \mathbf{s}_1 \Big| } ~. \end{split} \end{align}
Hamilton's Equations¶
In Jacobi coordinates Hamilton's equations are \begin{align} \dot{\mathbf{s}}_1 = \frac{ \partial H}{\partial \mathbf{p}_1} = \frac{ \mathbf{p}_1}{\mu_1}~, \end{align}
\begin{align} \dot{\mathbf{s}}_2 = \frac{ \partial H}{\partial \mathbf{p}_2} = \frac{ \mathbf{p}_2}{\mu_2}~, \end{align}
\begin{align} \begin{split} \dot{\mathbf{p}}_1 & = -\frac{\partial H}{\partial \mathbf{s}_1} \\ & = \frac{m_3}{m_1 + m_3}\frac{G m_1 m_2}{ \Big| \mathbf{s}_2 - \frac{m_3}{m_1 + m_3} \mathbf{s}_1\Big|^3} \left( \mathbf{s}_2 - \frac{m_3}{m_1 + m_3} \mathbf{s}_1\right) \\ & \;\;\;\;\; - \frac{G m_1 m_3}{\Big| \mathbf{s}_1 \Big|^3 } \mathbf{s}_1 - \frac{m_1}{m_1 + m_3}\frac{G m_2 m_3}{ \Big| \mathbf{s}_2 + \frac{m_1}{m_1 + m_3} \mathbf{s}_1 \Big|^3} \left( \mathbf{s}_2 + \frac{m_1}{m_1 + m_3} \mathbf{s}_1 \right) ~, \end{split} \end{align}
\begin{align} \begin{split} \dot{\mathbf{p}}_2 & = -\frac{\partial H}{\partial \mathbf{s}_2} \\ & = -\frac{G m_1 m_2}{ \Big| \mathbf{s}_2 - \frac{m_3}{m_1 + m_3} \mathbf{s}_1\Big|^3} \left( \mathbf{s}_2 - \frac{m_3}{m_1 + m_3} \mathbf{s}_1\right) - \frac{G m_2 m_3}{ \Big| \mathbf{s}_2 + \frac{m_1}{m_1 + m_3} \mathbf{s}_1 \Big|^3} \left( \mathbf{s}_2 + \frac{m_1}{m_1 + m_3} \mathbf{s}_1 \right) ~. \end{split} \end{align}
A good way for me to get the signs right was to calculate $\frac{\partial H}{\partial \mathbf{s}_1}= \frac{\partial H}{\partial \mathbf{q}_1}\frac{\partial \mathbf{q}_1}{\partial \mathbf{s}_1 } + \frac{\partial H}{\partial \mathbf{q}_2}\frac{\partial \mathbf{q}_2}{\partial \mathbf{s}_1 } + \frac{\partial H}{\partial \mathbf{q}_3}\frac{\partial \mathbf{q}_3}{\partial \mathbf{s}_1 }$ (similar for other terms).
Symplectic Integration¶
A good reference on symplectic integration is the book by Haier, Wanner, and Lubich.
What does symplectic mean (the briefest of explanations)?
Symplectic referes to a geometric property of Hamiltonian systems. In one dimension the 2-form $\omega= dp_1 \wedge dq^1$ describes an oriented area in phase space. More specifically, for two vectors \begin{align} X=\left( \begin{array}{c} X^p \\ X^q \end{array} \right)~\quad \&\quad ~Y=\left( \begin{array}{c} Y^p \\ Y^q \end{array} \right)~, \end{align} the symplectic 2-form on $X$ and $Y$ gives \begin{align} \begin{split} \omega(X,Y) & = \left( \begin{array}{c} X^p \\ X^q \end{array} \right)^T \left( \begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array} \right) \left( \begin{array}{c} Y^p \\ Y^q \end{array} \right) = X^pY^q - X^q Y^p \\ & =\text{the oriented area of the parallelogram spaned by $X~ \&~ Y$.} \end{split} \end{align}
In dimensions greater than one the symplectic 2-form is $\omega = dp_1 \wedge dq^1 + dp_2 \wedge dq^2 + ... + dp_n \wedge dq^n$. Evaluated on vectors $X$ and $Y$ the symplectic two form gives the sum of projected areas.
A sympelcitc transformation is one that preserves the symplectic structure, i.e., the quantitiy $\omega(X,Y)$ is invariant under a symplectic transformation.
The evolution of phase space is described by Hamilton's equations.
Hamiltonian evolution is symplectic, it is a very special kind of a symplectic (or canonical) transformation.
In one dimension we get the blob story (see page 89 in Isereles). A blob of phase space points desrcibes an area in phase space. Under Hamiltonian flow, this area is conserved.
What about conservation of energy?
- Symplectic integrators do not necessarily conserve energy. But, the error in energy from its original value tends to be small or bounded. Other integrators can accumulate errors.
- Symplectic integration is advantages for long time integration or when the step size is large - the generally do better at conservation of energy under these types of integration conditions.
The symplectic integrator presented in this notebook follows the paper of Yoshida 1993 section 3.2, which takes a Lie algebra perspective.
Hamilton's equations are \begin{align} \dot{q}^i= \frac{\partial H}{\partial p_i} ~, \quad \dot{p}_i = - \frac{\partial H}{\partial q^i} ~. \end{align} Letting $\eta=(q,p)$ Hamilton's equations can be written as
\begin{align} \begin{split} \frac{d\eta}{dt} &= \frac{ \partial \eta}{\partial q^i} \dot{q}^i + \frac{ \partial \eta}{\partial p_i} \dot{p}_i\\ & = \frac{ \partial \eta}{\partial q^i} \frac{\partial H}{\partial p_i} - \frac{ \partial \eta}{\partial p_i} \frac{\partial H}{\partial q^i} \\ & = \{ \eta, H \} ~, \end{split} \end{align} where the last line makes use of the Poisson bracket \begin{align} \{f, g\} = \frac{ \partial f}{\partial q^i} \frac{\partial g}{\partial p_i} - \frac{ \partial f}{\partial p_i} \frac{\partial g}{\partial q^i}~. \end{align}
Define the operator \begin{align} D_H = \frac{ \partial H}{\partial p_i} \frac{\partial }{\partial q^i} - \frac{\partial H}{\partial q^i} \frac{\partial }{\partial p_i} ~, \end{align} so that the time evolution of $\eta$ can be written as \begin{align} \frac{d \eta}{dt} = D_H ~\eta~. \end{align} For a time step of size $\tau$ the solution is an exponential \begin{align} \eta(\tau) = \exp\left[ \tau D_H \right] ~\eta(0) ~. \end{align}
Many Hamiltonians can be split into parts, for instance $H=H_0 + H_\text{interaction}$. For separable Hamiltonians the Hamiltonain operator can be split into a kinetic term and a potential term \begin{align} D_H = K_p + V_q~. \end{align} In this case the solution takes the form of \begin{align} \eta(\tau) = \exp \left[ \tau ( A + B) \right] ~\eta(0)~, \end{align} where $A= T_p$ (the kinetic term) and $B=V_p$ (the potential term). Up to order $\tau^n$ the solution can be determined by \begin{align} \exp(\tau[ A + B]) = \displaystyle \prod_{i}^n \exp(\tau~ c_i A)~ \exp(\tau~ d_i B) + o(\tau^{n+1}) ~, \end{align} for a set of real numbers $c_i$ and $d_i$.
A second order solution is $c_1 =1/2,~ c_2=1/2,~ d_1=1,~ d_2=0$, with expansion given by
\begin{align} \begin{split} \displaystyle \prod_{i}^2 \exp(\tau~ c_i A)~ \exp(\tau~ d_i B) + o(\tau^2) & =\exp(\tau/2~ A) ~\exp(\tau B)~ \exp(\tau/2~A) + o(\tau^3) \\ & =\big( 1+ \frac{\tau}{2} A + \frac{\tau^2}{4 2!} A^2 + ...\big)~\big( 1+ \tau B + \frac{\tau^2}{2!} B^2 + ...\big) ~\big( 1+ \frac{\tau}{2} A + \frac{\tau^2}{4 2!} A^2 + ...\big) + o(\tau^3) \\ & = 1 + \tau(A + B) + \frac{\tau^2}{2!} \left[ A^2 + [A,B] + B^2\right] + ... + o(\tau^3) \\ & = \exp(\tau~\left[ A+B\right]) + o(\tau^3)~. \end{split} \end{align}
For $A= \frac{\partial H}{\partial p}$ and $B= \frac{\partial H}{\partial q}$[where is the minus sign?], the above corresponds to a half step update in $q$ followed by a full step update in $p$ followed by a half step update in $q$, or \begin{align} \begin{split} & q_\star= q_0 + \frac{\tau}{2} \frac{\partial H(p_0,q_0)}{\partial p} ~, \\ & p=p_0 - \tau \frac{\partial H(p_0,q_\star)}{\partial q}~, \\ & q= q_\star + \frac{\tau}{2} \frac{\partial H(p,q_\star)}{\partial p}~, \end{split} \end{align} which is the Stormer-Verlet or leapfrog methods.
This notebook employs the third order method by Ruth \begin{align} \begin{split} & c_1 = \frac{7}{24}~, \quad ~ c_2 = \frac{3}{4} \\ & c_3=\frac{-1}{24} \\ & d_1 = \frac{2}{3} ~, \quad ~ d_2 = -d1 \\ & d_3=1~. \end{split} \end{align}
And, also the fourth order method by Ruth \begin{align} \begin{split} & c_1 = \frac{1}{2(2 - 2^{1/3})} ~, \quad ~ c_2 = \frac{1-2^{1/3}}{2(2 - 2^{1/3})} \\ & c_3=c_2 ~, \quad ~ c_4=c_1 \\ & d_1 = \frac{1}{2(2 - 2^{1/3})} ~, \quad ~ d_2 = \frac{-2^{1/3}}{2(2 - 2^{1/3})} \\ & d_3=d_1 ~, \quad ~ d_4=0 ~. \end{split} \end{align}
Julia code¶
function s1_dot(p1::Vector{Float64},p2::Vector{Float64}, params::Array{Float64})
m1, m2, m3, mu1, mu2=params
return p1/mu1
end;
function s2_dot(p1::Vector{Float64},p2::Vector{Float64}, params::Array{Float64})
m1, m2, m3, mu1, mu2=params
return p2/mu2
end;
function p1_dot(s1::Vector{Float64}, s2::Vector{Float64}, params::Array{Float64})
m1, m2, m3, mu1, mu2=params
q23=s2+ m1/(m1+m3)*s1
q21=s2- m3/(m1+m3)*s1
result= m1*m2*mu1/m1*q21/norm(q21)^3. - m2*m3*mu1/m3*q23/norm(q23)^3. - m1*m3/norm(s1)^3*s1
return result
end;
function p2_dot(s1::Vector{Float64}, s2::Vector{Float64},params::Array{Float64})
m1, m2, m3, mu1, mu2=params
q23=s2+ m1/(m1+m3)*s1
q21=s2- m3/(m1+m3)*s1
result= -m1*m2/norm(q21)^3.*q21 - m2*m3/norm(q23)^3.*q23
return result
end;
function map_to_Jacobi(q1::Vector{Float64},q2::Vector{Float64},
q3::Vector{Float64},params::Vector{Float64})
m1, m2, m3, mu1, mu2=params
s1=q1-q3
s2=(-m1*q1-m3*q3)/(m1+m3) + q2
R=(m1*q1+m2*q2+m3*q3)/(m1+m2+m3)
return s1, s2, R
end;
function map_to_Cartesian(s1::Vector{Float64},s2::Vector{Float64},
R::Vector{Float64},params::Vector{Float64})
m1, m2, m3, mu1, mu2=params
M=m1+m2+m3
q1=m3/(m1+m3)*s1 - m2/M*s2 + R
q2=(m1+m3)/M*s2 + R
q3=-m1/(m1+m3)*s1 -m2/M*s2 + R
return q1, q2, q3
end;
function get_mu1(m1,m2,m3)
return m1*m3/(m1+m3)
end;
function get_mu2(m1,m2,m3)
return m2*(m1+m3)/(m1+m2+m3)
end;
function dsdt(s,p)
p1=p[1:3]
p2=p[4:6]
return vcat(s1_dot(p1,p2,params), s2_dot(p1,p2,params))
end;
function dpdt(s,p)
s1=s[1:3]
s2=s[4:6]
return vcat(p1_dot(s1,s2,params), p2_dot(s1,s2,params))
end;
function H(s,p,params)
m1, m2, m3, mu1, mu2=params
p1=p[1:3]
p2=p[4:6]
s1=s[1:3]
s2=s[4:6]
q23=s2+ m1/(m1+m3)*s1
q21=s2- m3/(m1+m3)*s1
K=1/2.*norm(p[1:3])^2/mu1 + 1/2.*norm(p[4:6])^2/mu2
Phi=- m1*m2/norm(q21)-m1*m3/norm(s1)-m2*m3/norm(q23)
return K + Phi
end;
function s_to_q(s,R,params)
n=size(s,1)
println(size(s))
q=zeros(n,9)
for i in 1:size(s,1)
q[i,1:3],q[i,4:6],q[i,7:9]=map_to_Cartesian(s[i,1:3],s[i,4:6], R, params)
end
return q
end;
function Ruth3rd(dqdt::Function,dpdt::Function,q_0::Vector{Float64},p_0::Vector{Float64},tau::Float64)
c1=7./24.
c2=3./4.
c3=-1./24.
d1=2./3.
d2=-d1
d3=1.
# the sign in front of d_i is positive
# since I use dpdt instead of V_q
q_star=q_0 + c1*dqdt(q_0,p_0)*tau
p_star=p_0 + d1*dpdt(q_star,p_0)*tau
q_star=q_star + c2*dqdt(q_star,p_star)*tau
p_star=p_star + d2*dpdt(q_star,p_star)*tau
q_star=q_star + c3*dqdt(q_star,p_star)*tau
p_star=p_star + d3*dpdt(q_star,p_star)*tau
return q_star, p_star
end;
function Ruth4th(dqdt::Function,dpdt::Function,q_0::Vector{Float64},p_0::Vector{Float64},tau::Float64)
c1=1/2./(2-2^1/3.)
c2=(1-2^1/3.)/2./(2-2^1/3.)
c3=c2
c4=c1
d1=1/(2-2^1/3.)
d2=-2^1/3./(2-2^1/3.)
d3=d1
#d4=0
# the sign in front of d_i is positive
# since I use dpdt instead of V_q
q_star=q_0 + c1*dqdt(q_0,p_0)*tau
p_star=p_0 + d1*dpdt(q_star,p_0)*tau
q_star=q_star + c2*dqdt(q_star,p_star)*tau
p_star=p_star + d2*dpdt(q_star,p_star)*tau
q_star=q_star + c3*dqdt(q_star,p_star)*tau
p_star=p_star + d3*dpdt(q_star,p_star)*tau
q_star=q_star + c4*dqdt(q_star,p_star)*tau
return q_star, p_star
end;
function Hamiltonian_int(q_0,p_0,dqdt,dpdt,time,Integrator)
n=size(q_0,1)*2 # number of phase space varaibles
n_t=size(time,1) # number of time steps
q=zeros(n_t,n//2)
p=zeros(n_t,n//2)
q[1,:]=q_0
p[1,:]=p_0
for t in 2:n_t
dt=time[t]-time[t-1]
q[t,:],p[t,:]=Integrator(dqdt,dpdt,q[t-1,:],p[t-1,:],dt)
end
return q, p
end;
Figure Eight¶
Some interesting solutions of the three body problem can be found on the associated scholarpedia website. Amoung them is the figure eight solution.
# initial parameters for figure eight
m1=1.0
m2=1.0
m3=1.0
mu1=get_mu1(m1,m2,m3)
mu2=get_mu2(m1,m2,m3)
params=[m1,m2,m3,mu1,mu2]
q1_0=[-0.97000436; 0.24308753; 0.0]
q2_0=[0.97000436; -0.24308753; 0.0]
q3_0=[0.0 ;0.0; 0.0]
q3_dot_0=[0.93240737; 0.8647314; 0.0]
q2_dot_0=[-0.46620368; -0.4323657; 0.0]
q1_dot_0=[-0.46620368 ; -0.4323657 ; 0.0]
s1,s2, R=map_to_Jacobi(q1_0,q2_0,q3_0, params)
v1,v2,v_R=map_to_Jacobi(q1_dot_0,q2_dot_0, q3_dot_0, params)
p1=v1*mu1;
p2=v2*mu2;
s_0=vcat(s1,s2);
p_0=vcat(p1,p2);
time=linspace(0,50,5000)
@time s, p=Hamiltonian_int(s_0,p_0,dsdt,dpdt,time,Ruth4th);
println("integrated")
q=s_to_q(s,R,params);
n_t=1000
save_name="figure8.mp4"
using PyPlot
using PyCall
@pyimport matplotlib.animation as anim
function showanim(filename)
base64_video = base64encode(open(read, filename))
display("text/html", """<video controls src="data:video/x-m4v;base64,$base64_video">""")
end
fig = figure(figsize=(6,6))
ax = axes(xlim = (-1.5,1.5),ylim=(-1,1))
ax[:axis]("off")
global line1 = ax[:plot]([],[],"black", alpha=.5)[1]
global line2 = ax[:plot]([],[],"black", alpha=.5)[1]
global line3 = ax[:plot]([],[],"black", alpha=.5)[1]
global pt1 = ax[:plot]([],[],"o", color="black")[1]
global pt2 = ax[:plot]([],[],"o",color="black")[1]
global pt3 = ax[:plot]([],[],"o",color="black")[1]
function init()
global line1
global line2
global line3
global pt1
line1[:set_data]([],[])
line2[:set_data]([],[])
line3[:set_data]([],[])
pt1[:set_data]([],[])
pt2[:set_data]([],[])
pt3[:set_data]([],[])
return (line1,line2,line3,pt1,pt2,pt3,Union{}) # Union{} is the new word for None
end
function animate(i)
# because FuncAnimation starts indexing at zero
# and Julia starts indexing at 1
i=i+1
global line1
global line2
global line3
global pt1
line1[:set_data](q[1:i,1],q[1:i,2])
pt1[:set_data](q[i,1],q[i,2])
line2[:set_data](q[1:i,4],q[1:i,5])
pt2[:set_data](q[i,4],q[i,5])
line3[:set_data](q[1:i,7],q[1:i,8])
pt3[:set_data](q[i,7],q[i,8])
return (line1,line2,line3,pt1,pt2,pt3,Union{})
end
withfig(fig, clear=false) do
# Create the animation object by calling the Python function FuncAnimaton
myanim = anim.FuncAnimation(fig, animate, init_func=init, frames=n_t, interval=20)
# Convert it to an MP4 movie file and saved on disk in this format.
myanim[:save](save_name, bitrate=-1, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
end;
# Function for creating an embedded video given a filename
function html_video(filename)
#empty!(IJulia.displayqueue)
open(filename) do f
base64_video = base64encode(f)
"""<video controls src="data:video/x-m4v;base64,$base64_video">"""
end
end
# Display the movie in a Julia cell as follows. Note it has animation controls for the user.
#display("text/html", html_video(save_name))
showanim(save_name)
function anim_plot(N,q,lim,save_name)
n=3
colors=["blue", "red", "green"]
# find way to set axis
fig = figure()
ax =axes(xlim = (-lim,lim),ylim=(-lim,lim),zlim=(-lim,lim), projection="3d")
ax[:axis]("off")
ax[:set_facecolor]("black")
global lines = Array{Any}(n)
global pts = Array{Any}(n)
for i = 1:n
lines[i] = ax[:plot]([],[],[],color = colors[i], alpha=.8)[1]
pts[i]= ax[:plot]([],[],[],"o", color = colors[i])[1]
end
function init()
global lines
global pts
for i = 1:n
lines[i][:set_data]([],[])
lines[i][:set_3d_properties]([])
pts[i][:set_data]([],[])
pts[i][:set_3d_properties]([])
end
result = tuple(tuple([lines[i] for i = 1:n]),tuple([pts[i] for i = 1:n]))
return result
end
function animate(i)
i=i+1
global lines
global pts
for k = 1:n
j=1+(k-1)*3
x = q[1:i,j]
y = q[1:i,j+1]
z = q[1:i,j+2]
lines[k][:set_data](x,y)
lines[k][:set_3d_properties](z)
pts[k][:set_data](last(x),last(y))
pts[k][:set_3d_properties](last(z))
#Rotate axes
ax[:view_init](30, 0.3 * i)
end
result = tuple(tuple([lines[i] for i = 1:n]),tuple([pts[i] for i = 1:n]))
return result
end
# set clear equal to false or you
# don't get a plot
withfig(fig, clear=false) do
myanim = anim.FuncAnimation(fig, animate, init_func=init, frames=N, interval=15)
# Convert it to an MP4 movie file and saved on disk in this format.
myanim[:save](save_name, bitrate=-1, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
end
display("text/html", html_video(save_name))
end
Fun with initial conditions¶
# random selection of initial parameters
m1=.8
m2=1.5
m3=4.5
mu1=get_mu1(m1,m2,m3)
mu2=get_mu2(m1,m2,m3)
params=[m1,m2,m3,mu1,mu2]
q1_0=[-1.0; 0.5; -2.0]
q2_0=[1.0; -1.0; 3.0]
q3_0=[0.0 ;0.0; 0.0]
q3_dot_0=[1.0; 0.8; 0.5]
q2_dot_0=[-0.5; -0.2; 0.0]
q1_dot_0=[-0.5 ; -0.1 ; 0.0]
s1,s2, R=map_to_Jacobi(q1_0,q2_0,q3_0, params)
v1,v2,v_R=map_to_Jacobi(q1_dot_0,q2_dot_0, q3_dot_0, params)
p1=v1*mu1;
p2=v2*mu2;
s_0=vcat(s1,s2);
p_0=vcat(p1,p2);
println("Initial energy ", H(s_0,p_0,params))
#time=linspace(0,350,25600)
time=linspace(0,350,24001)
@time s, p=Hamiltonian_int(s_0,p_0,dsdt,dpdt,time,Ruth4th);
q=s_to_q(s,R,params);
q=q[1:8:end,:];
println("size of q", size(q))
#anim_plot(3200,q, 3.5,"3body_run1.mp4")
anim_plot(3001,q, 3.5,"3body_run1.mp4")
Symplectic versus RK4¶
As mentioned in the notes, symplectic integrators are good for long time integrations. As an example, I integrated the Keppler problem for many orbits with a large time step using the Ruth 3rd order scheme and the Runge-Kutta 4th order (RK4) scheme. The code is below and the results are plotted. The plots show that the errors in the RK4 method grows. The error in the Ruth 3rd order method is much more stable (even though it is 3rd order.)
function H_Keppler(q,p)
return .5*(p[1]^2 + p[2]^2)- 1/norm(q)
end;
function dpdt_Kep(q,p)
return -1/norm(q)^2*q
end;
function dqdt_Kep(q,p)
return p
end;
function RHS_Kep(Y)
q=Y[1:2]
p=Y[3:4]
return vcat(dqdt_Kep(q,p),dpdt_Kep(q,p))
end;
function RK4(F,Y,dt)
k1=F(Y)
k2=F(Y + 1/2.*dt*k1)
k3=F(Y + 1/2.*dt*k2)
k4=F(Y + dt*k3)
return Y + 1/6.*dt*(k1 + 2*k2 + 2*k3 + k4)
end;
q_0=[1,0]
p_0=[0,1]
n_t=650
time=linspace(0,300,n_t)
H_rk4=zeros(time)
H_symp=zeros(time)
H_rk4[1]=H_Keppler(q_0,p_0)
H_symp[1]=H_rk4[1]
Y=zeros(n_t,4)
q=zeros(n_t,2)
p=zeros(n_t,2)
q[1,:]=q_0
p[1,:]=p_0
Y[1,1:2]=q_0
Y[1,3:4]=p_0
Y_0=Y[1,:]
for i in 2:n_t
dt=time[i]-time[i-1]
Y[i,:]=RK4(RHS_Kep,Y_0,dt)
q[i,:],p[i,:]=Ruth3rd(dqdt_Kep, dpdt_Kep, q[i-1,:],p[i-1,:],dt)
Y_0=Y[i,:]
H_rk4[i]=H_Keppler(Y_0[1:2],Y_0[3:4])
H_symp[i]=H_Keppler(q[i,:],p[i,:])
end;
fig = figure(figsize=(14,6))
subplot(121)
ax=gca()
xlabel("X", size=20)
ylabel("Y", size=20)
ax[:grid]()
ax[:plot](Y[:,1],Y[:,2],"--", color="black", alpha=.5,label="RK4")
ax[:plot](q[:,1],q[:,2], color="blue",label="Ruth 3rd")
ax[:legend](loc="upper left", fontsize="x-large")
subplot(122)
ax=gca()
xlabel("time", size=20)
ylabel(L"$E/E_0$", size=20)
ax[:grid]()
ax[:plot](time,H_rk4/H_rk4[1],"--",color="black", label="RK4")
ax[:plot](time,H_symp/H_symp[1],color="blue", label="Ruth 3rd")
ax[:legend](loc="upper left", fontsize="x-large")
savefig("Kepler_error.pdf")
