Notebook

Ordinary Differential Equations (ODE)¶

Ordinary differential equations (ODEs) arise in many fields of physics, but also in social and other natural sciences. Often, quantities are defined as the rate of change of other quantities (for example, derivatives of displacement with respect to time $\dot{x}(t)$, $\dot{y}(t)$ are related to $x(t)$, $y(t)$ and $t$). This results in ODE.

Definition¶

If differential equation contains only one independent variable (such as time $t$) and hence only derivatives with respect to this single independent variable are present (like $\dot{x}$, $\ddot{x}$,...), we classify it as an ordinary differential equation (ODE). If derivatives with respect to more than one independent variables appear, we call it a partial differential equation. The latter is more challenging to solve (for example using Finite element method) and we will not discuss it here.

In ODE only total derivatives (and no partial derivatives) appear. In general, it might be written in a form $$F(t, y(t), \frac{d y(t)}{dt}, \frac{d^2 y(t)}{dt^2}, ... \frac{d^n y(t)}{dt^n})=0$$ where $F$ is an arbitrary function, and it could be non-linear in its arguments.

One of the simplest examples is the Newton's law:

$$m \frac{d^2 x_i(t)}{dt^2} = F_i(x_1,x_2,x_3,t)$$

Here $(x_1,x_2,x_3)=(x,y,z)$ and $F_i$ is the force.

Classification¶

ODE's are usually classified into:

Linear ODE : in which F is a linear function of its arguments, for example $$ \sum_{n=0}^N a_n \frac{d^n y(t)}{dt^n} -r (t)= 0 $$ Here $r(t)$ is an arbitrary function, which we usually call source term.
Homogeneous ODE : In this case source term $r(t)=0$ and hence $y(t)=0$ is also a trivial solution. But note that homogeneous ODE can be non-linear.
Autonomous ODE : If it does not explicitely depent on $t$.
Non-linear ODE : at least one derivative appears at higher power than 1, i.e., $\left(\frac{d^n y(t)}{dt^n}\right)^2$.
inhomogeneous : r(t) is nonzero.

What can we say about the Newton's equation and this classification?

Numerical solvers¶

Many general purpose numerical solvers have been developed, which are easy to use, and can solve an arbitrary ODE.

We will discuss a few below: Runge Kutta method, Numerov, Verlet, etc.

While very generic solvers exist in scipy.integrate.ode, the precision of the solution can be a problem. For example, it is hard to follow Newton's eq. of motion for a very long time as the precision might deteriorate (molecular dynamics is challenging). If ODE is non-linear, it tends to be even more challenging to solve numerically. For example, the numerical solution might exhibit instability in some regimes, which is called stiffness. This is often caused by the presence of different time scales in the underlying problem. For example, a non-elastic collision in a mechanical system typically occurs at much smaller time scale than the time for the motion of objects, and this discrepancy causes a very "sharp turns" in the curves of the state parameters. These are difficult to capture numerically.

Reduction to a first-order system¶

Most ODE integration algorithms (integrators) solve first order ODE's only, but with any number of components. It is simple to see that any high order ODE can be reduced to a larger system of first order ODE.

For a generic form $$F(t, y(t), \frac{d y(t)}{dt}, \frac{d^2 y(t)}{dt^2}, ... \frac{d^n y(t)}{dt^n})=0$$ we could, always rewrite into $$\frac{d^n y(t)}{dt^n} = \widetilde{F}(t, y(t), \frac{d y(t)}{dt}, \frac{d^2 y(t)}{dt^2}, ...\frac{d^{n-1} y(t)}{dt^{n-1}})$$

Than we can choose \begin{eqnarray} &&y(t) \equiv y_0(t) \\ &&\frac{d y(t)}{dt}=\frac{d y_0(t)}{dt}\equiv y_1(t)\\ &&\frac{d^2 y(t)}{dt^2}=\frac{d y_1(t)}{dt}\equiv y_2(t)\\ &&\frac{d^3 y(t)}{dt^3}=\frac{d y_2(t)}{dt}\equiv y_3(t)\\ &&...\\ &&\frac{d^n y(t)}{dt^n}=\frac{d y_{n-1}(t)}{dt}=\widetilde{F}(t, y_0(t), y_1(t), y_2(t), ...y_{n-1}(t)) \end{eqnarray} This proves that $n$-th order ODE is equivalent to $n$ component first order ODE. Here $y(t)$ could itself be multi-component function, for example $(x,y,z)$, and we would than have $3 n$ component first order ODE.

A good example of such transformation in physics is transformation of Lagrangian $L$ equations, which involve the second order derivatives $\ddot{x}(t)$, to Hamiltonian $H$ equations, which involve only $\dot{p}(t)$ and $\dot{x}(t)$.

Since it is always possible to transform n-other ODE into n-component first order ODE, we are going to discuss numeric treatment of the first order ODE of the form:

$$\frac{d y(t)}{dt}=\overline{F}(t, y(t))$$

where $y=[y_0,y_1,...y_{n-1}]$ and $\overline{F}=[\dot{y}_0,\dot{y}_1,...,\widetilde{F}(t,y_0,....y_{n-1})]$

Boundary conditions¶

To solve ODE, the boundary conditions need to be specified. For $n$-th component first order ODE (or n-th order ODE) we need exactly $n$ boundary condistions.

Boundary conditions can be classified into:

Initial value problems : all necessary conditions are specified at the starting point ($t_0$):

\begin{eqnarray} &&y_0(t_0)=x_0\\ &&y_1(t_0)=x_1\\ &&...\\ &&y_{n-1}(t_0)=x_{n-1} \end{eqnarray}

For example, in projectile motion, we give the initial position $x_0=y(t_0)$ and the initial velocity $v_0=dy(t_0)/dt$.

Two point boundary problems : part of the conditions are specified at the starting point ($t_0$), and the rest at the end point ($t_f$)

For example, projectile motion, in which we are given the initial position $x_0=y(t_0)$ and the end position $x_f=y(t_f)$. The initial velocity needs to be guessed by so-called shooting method.
More complicated boundary conditions : an arbitrary set of nonlinear equations relating values of $y_i(t)$ and their derivatives $d y_i(t)/dt$ at any points $t$.

Standard numeric solvers can only solve Initial value problem. All other boundary conditions have to be solved by iterations or shooting, in which we attempt to solve ODE many times, and try to satisfy the boundary conditions better and better with iterations.

Euler's method¶

From any point on a curve, you can find an approximation of a nearby point on the curve by moving a short distance along a line tangent to the curve.

We replace the derivative $\frac{dy(t)}{dt}$ by the finite difference approximation: $$\frac{dy(t)}{dt}\approx \frac{y(t+h)-y(t)}{h}$$ When inserting this approximation into ODE, we get $$y(t+h) \approx y(t) + h \overline{F}(t, y(t))$$

If we choose a constant step $h$, so that $t_1 = t_0+h$, $t_2=t_0+2 h$,..., $t_n=t_0+n h=t_f$, and we define $y_i \equiv y(t_i)=y(t_0+i h)$ we have simple recurrence relation:

$$y_{i+1}=y_i + h \overline{F}(t, y_i) + O(h^2)$$

which shows that each step give an error of the order $h^2$. The error grows very quickly, and Euler's method is not being used in practice.

Runge-Kutta methods¶

The second order Runge-Kutta is only slightly better than Euler's, but demonstrates the idea of this method.

In Euler's method the derivative is taken at the beginning of the interval. The precision would increase if one could estimate derivative in the middle of the interval $$y_{i+1}=y_i + h \overline{F}(t_i+h/2,y(t_i + h/2))$$.

The second order Runge-Kutta (RK2) method implements precisely this idea, namely, \begin{eqnarray} && k_1 = h \overline{F}(t_i,y_i)\\ && y_{i+1} = y_i + h \overline{F}(t_i+\frac{1}{2}h,y_i + \frac{1}{2}k_1) + O(h^3) \end{eqnarray} It is called the second order, because error is of the order of $h^3$. In general, the error of the $n$-th order routine is $O(h^{n+1})$.

Most popular is the forth-order Runge Kutta (RK4) method, which builds on the above idea of evaluating derivative somewhere inside the interval: \begin{eqnarray} && k_1 = h \overline{F}(t_i,y_i)\\ && k_2 = h \overline{F}(t_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1)\\ && k_3 = h \overline{F}(t_i+\frac{1}{2}h,y_i+\frac{1}{2}k_2)\\ && k_4 = h \overline{F}(t_i+h, y_i+k_3)\\ && y_{i+1} = y_i + \frac{1}{6}k_1+\frac{1}{3}k_2+\frac{1}{3}k_3+\frac{1}{6}k_4+O(h^5) \end{eqnarray}

How do we understand the method? Looking at the above figure, we see:

$k_1$ is the slope at the beginning of the interval;
$k_2$ is the slope at the midpoint of the interval, using slope $k_1$ to determine the value of $y$ at the point $t_i + h/2$ using Euler's method;
$k_3$ is again the slope at the midpoint, but now using the slope $k_2$ to determine the $y$-value;
$k_4$ is the slope at the end of the interval, with its $y$-value determined using $k_3$;
in averaging the four slopes, greater weight is given to the slopes at the midpoint.

The RK4 method is a fourth-order method, meaning that the error per step is on the order of $h^5$, while the total accumulated error has order $h^4$. With only four function evaluations, for fourth order accuracy is extremely good.

Example: Bouncing Ball¶

The bouncing ball has Newton's Eq: $$ \frac{d^2 y(t)}{dt^2}=-g$$ which can be reduce to two first order ODE's: $$\frac{dy}{dt} = v$$ $$\frac{dv}{dt} = -g$$

The initial condition is $v(0)=0$ and $y(0)=y_0$.

We will first solve this with direct Euler's method. Next we will use RK4 and other more advanced method to get more precise solution.

Euler's method requires: \begin{eqnarray} &&y(t+h) = y(t) + h\; v(t)\\ &&v(t+h) = v(t) - h\; g \end{eqnarray}

Free fall first¶

In [1]:

import numpy as np

g = 9.81

# intital conditions
y = 10.
v = 0.0

ti = 0
tf = 2.
Nt=201

# fixed time point
t = np.linspace(ti,tf,Nt)
dt=t[1]-t[0]

data = np.zeros((len(t),2))
data[0,:]=[y,v]
for i in range(Nt-1):
    y = y + v*dt
    v = v - g*dt
    data[i+1,:] = [y, v] 

In [2]:

import matplotlib.pyplot as plt

plt.plot(t, data[:,0],'.-')

Out[2]:

[<matplotlib.lines.Line2D at 0x7fbad8e1af10>]

Bouncing¶

Add a floor at $y=0$.

What happens at the floor? – The ball bounces back elastically, hence velocity is reversed on collision $v\rightarrow -v$.

In [3]:

import numpy as np

g = 9.81

# intital conditions
y = 10.
v = 0.0
y_floor = 0


ti = 0
tf = 20.
Nt=2001

# fixed time point
t = np.linspace(ti,tf,Nt)
dt=t[1]-t[0]

data = np.zeros((len(t),2))
data[0,:]=[y,v]
for i in range(Nt-1):
    y = y + v*dt
    if y > y_floor:
        v = v - g*dt
    else:
        v = -v   # bounce off
    data[i+1,:] = [y, v] 
 

In [4]:

plt.plot(t, data[:,0])
plt.xlabel("time (s)")
plt.ylabel("position (m)");

Example: Oscilators¶

Harmonic oscilator : $V=\frac{1}{2} k x^2$ and $F=-\frac{dV(x)}{dx}=-k x$

Anharmonic oscilator : $V=\frac{1}{2} k x^2 (1-\frac{2}{3}\alpha x)$ and $F=-kx (1-\alpha x)$

$$m \ddot{x} = F$$

Harmonic oscilator : $$\ddot{x}+ \omega^2 x=0$$ Anharmonic oscilator: $$\ddot{x}+\omega^2 x(1-\alpha x)=0$$ where we used $\omega=\sqrt{\frac{k}{m}}$ and we can just set $\alpha=0$ to recover harmonic oscilator from anharmonic.

Corresponding First order Eq: \begin{eqnarray} &&\frac{d x(t)}{dt} = v\\ &&\frac{d v(t)}{dt}= -\omega^2 x (1-\alpha x) \end{eqnarray}

Total energy should be conserved. The expression for total energy is $$E=T+V = \frac{1}{2} m \dot{x}^2+k(\frac{1}{2} x^2-\frac{1}{3}\alpha x^3),$$ which can also be written as $$\frac{E}{m} = \frac{1}{2}\dot{x}^2+\omega^2(\frac{1}{2} x^2-\frac{1}{3}\alpha x^3)$$

Recall the integration methods:¶

Euler:

\begin{eqnarray} y_{i+1}=y_i + h \overline{F}(t, y_i) + O(h^2) \end{eqnarray}

RK2:

\begin{eqnarray} && k_1 = h \overline{F}(t_i,y_i)\\ && k_2 = h \overline{F}(t_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1)\\ && y_{i+1} = y_i + h k_2+O(h^3) \end{eqnarray}

RK4:

\begin{eqnarray} && k_1 = h \overline{F}(t_i,y_i)\\ && k_2 = h \overline{F}(t_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1)\\ && k_3 = h \overline{F}(t_i+\frac{1}{2}h,y_i+\frac{1}{2}k_2)\\ && k_4 = h \overline{F}(t_i+h, y_i+k_3)\\ && y_{i+1} = y_i + \frac{1}{6}k_1+\frac{1}{3}k_2+\frac{1}{3}k_3+\frac{1}{6}k_4+O(h^5) \end{eqnarray}

In [5]:

## Three simple integrators implementations

def euler(y, f, t, h):
    """Euler integrator.
    Returns new y(t+h), where y=y(t).
    """
    return y + h * f(t, y)

def rk2(y, f, t, h):
    """Runge-Kutta RK2 midpoint"""
    k1 = f(t, y)
    k2 = f(t + 0.5*h, y + 0.5*h*k1)
    return y + h*k2

def rk4(y, f, t, h):
    """Runge-Kutta RK4"""
    k1 = f(t, y)
    k2 = f(t + 0.5*h, y + 0.5*h*k1)
    k3 = f(t + 0.5*h, y + 0.5*h*k2)
    k4 = f(t + h, y + h*k3)
    return y + h/6 * (k1 + 2*k2 + 2*k3 + k4)

Recall First order Eq: \begin{eqnarray} \begin{bmatrix}\frac{d x(t)}{dt}\\ \frac{d v(t)}{dt}\end{bmatrix}= \begin{bmatrix}v \\ -\omega^2 x (1-\alpha x)\end{bmatrix} \end{eqnarray}

Total energy : $$\frac{E}{m} = \frac{1}{2}v^2+\frac{1}{2} x^2 \omega^2(1-\frac{2}{3}\alpha\; x)$$

In [6]:

def Harmonic(t, y, w2=1):
    """Harmonic EOM"""
    return np.array([y[1],-w2*y[0]])

def anHarmonic(t, y, alpha=0.5, w2=1):
    """anharmonic EOM"""
    return np.array([y[1],-w2*y[0]*(1-alpha*y[0])])

def E_harmonic(y, w2=1):
    """Harmonic total energy  E(v,x)/m = 1/2*(v^2+ w^2 x^2) """
    T = 0.5*y[1]**2
    V = 0.5*w2*y[0]**2
    return np.array([T+V,T,V])

def E_anharmonic(y, alpha=0.5, w2=1):
    """Anharmonic total energy E(v,x)/m = 1/2*(v^2 + w^2 x^2 (1 - 2/3 alpha x)"""
    T = 0.5*y[1]**2
    V = 0.5*w2*y[0]**2*(1 - 2./3.*alpha*y[0])
    return np.array([T+V,T,V])

Next we code a generic Solve for fixed step integrators, which we just implemented. It will take arguments solver, which can be euler or rk4 and derivs, which can be Harmonic or anHarmonic.

In [7]:

def Solve(t, y0, solver=euler, derivs=Harmonic):
    """ t is independent variable, which needs to be a predefined mesh of points (could be non-equidistant).
        y0 is the initial value for [x,v]
        solver can be euler/rk2/rk4/verlet, etc
        derivs needs to provide derivatives, for example Harmonic or anHarmonic
    """
    data = np.zeros((len(t),2))  # storage for solution
    y = np.array(y0)
    data[0,:]=y
    for i in range(len(t)-1):
        y = solver(y,derivs,t[i],t[i+1]-t[i])
        data[i+1,:] = y 
    return data

# Next we evaluate for Harmonic oscilator using Euler and RK4
# intital conditions
x = 0    # initial position
v = 1.0  # initial velocity
# time mesh
ti = 0    # initial time
tf = 100. # end time
Nt=2001   # number of points, gives step of 0.05
t = np.linspace(ti,tf,Nt)

dataE = Solve(t, [x,v], euler, Harmonic)
dataR = Solve(t, [x,v], rk4, Harmonic)

Ene = E_harmonic(dataE.T)
Enr = E_harmonic(dataR.T)

In [8]:

plt.plot(t, dataE[:,0], label='Euler')
plt.plot(t, dataR[:,0], label='RK4')
plt.plot(t, np.sin(t), label='exact')
plt.legend(loc='best');

We notice that $RK4$ is orders of magnitude better than Euler in solving the harmonic oscillator (HO). Euler's method increases the amplitude of oscillations, which means that energy is erroneously increasing.

Next we plot the energy change with time

In [9]:

E0=Ene[0,0]
plt.plot(t,Ene[0]-E0)
plt.plot(t,Enr[0]-E0);

Indeed, the energy shows huge increase in Euler's method, but mostly constant in RK4. Using Euler's method we are clearly not solving the equations of motions in a physically meaningful way.

Next we plot RK4 energy alone to see the error better:

In [10]:

plt.plot(t,Enr[0]-E0)
print(Enr[0,-1]-E0)

-2.1694602519994888e-07

The error of RK4 method is small, and is excellent for small times. Of course the error depends on timestep, but for 2000 function evaluations, this is small error.

But we do see a problem: The energy is monotonically decreasing. If we simulate Newton's motion long enough we will loose all energy, and system will halt. This is very unphysical, therefore RK4 is not acceptable for molecular dynamics simulation.

Next we check the anharmonic oscilator.

In [11]:

data2E = Solve(t, [x,v], euler, anHarmonic)
data2R = Solve(t, [x,v], rk4, anHarmonic)

/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/ipykernel_56755/1125370763.py:7: RuntimeWarning: overflow encountered in scalar multiply
  return np.array([y[1],-w2*y[0]*(1-alpha*y[0])])

In [12]:

plt.plot(t, data2E[:,0], label='x:Euler')
plt.plot(t, data2R[:,0], label='x:RK4')
plt.plot(t, data2R[:,1], label='v:RK4')
plt.legend(loc='best')
plt.ylim([-1.2,1.7])

Out[12]:

(-1.2, 1.7)

The Euler solution becomes unstable after just one period, while RK4 remains reasonable many periods.

Next we plot Phase-space portrait $(x,p)$ for HO and anharmonic oscilator.

In [13]:

plt.plot(dataR[:,0], dataR[:,1], label='Harmonic')
plt.plot(data2R[:,0],data2R[:,1], label='anHarmonic')
plt.legend(loc='best')

Out[13]:

<matplotlib.legend.Legend at 0x7fbb1a020070>

How good is the energy of the ahharmonic oscilator?

In [14]:

# Energy for the anharmonic oscilator
En2 = E_anharmonic(data2R.T)
plt.plot(t, En2[1], 'r-', label='$E_{kin}$')
plt.plot(t, En2[2], 'b-', label='$E_{pot}$')
plt.plot(t, En2[0], 'k-', label='$E_{tot}$')
plt.legend(loc='best')

Out[14]:

<matplotlib.legend.Legend at 0x7fbb1a020ee0>

But zoom-in shows that the algorithm is loosing energy with time. The error is small with small time-step, but long time unstable.

In [15]:

E0=En2[0,0]
plt.plot(t,En2[0]-E0)
print(En2[0,-1]-E0)

-3.023807507718246e-07

Let's try scipy recomended solver scipy.solve_ivp, which is also RK4, but with variable step and more powerful algorithm RK45, oe even eith order RK (DOP853).

In [16]:

from scipy.integrate import solve_ivp
help(solve_ivp)

Help on function solve_ivp in module scipy.integrate._ivp.ivp:

solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)
    Solve an initial value problem for a system of ODEs.
    
    This function numerically integrates a system of ordinary differential
    equations given an initial value::
    
        dy / dt = f(t, y)
        y(t0) = y0
    
    Here t is a 1-D independent variable (time), y(t) is an
    N-D vector-valued function (state), and an N-D
    vector-valued function f(t, y) determines the differential equations.
    The goal is to find y(t) approximately satisfying the differential
    equations, given an initial value y(t0)=y0.
    
    Some of the solvers support integration in the complex domain, but note
    that for stiff ODE solvers, the right-hand side must be
    complex-differentiable (satisfy Cauchy-Riemann equations [11]_).
    To solve a problem in the complex domain, pass y0 with a complex data type.
    Another option always available is to rewrite your problem for real and
    imaginary parts separately.
    
    Parameters
    ----------
    fun : callable
        Right-hand side of the system. The calling signature is ``fun(t, y)``.
        Here `t` is a scalar, and there are two options for the ndarray `y`:
        It can either have shape (n,); then `fun` must return array_like with
        shape (n,). Alternatively, it can have shape (n, k); then `fun`
        must return an array_like with shape (n, k), i.e., each column
        corresponds to a single column in `y`. The choice between the two
        options is determined by `vectorized` argument (see below). The
        vectorized implementation allows a faster approximation of the Jacobian
        by finite differences (required for stiff solvers).
    t_span : 2-tuple of floats
        Interval of integration (t0, tf). The solver starts with t=t0 and
        integrates until it reaches t=tf.
    y0 : array_like, shape (n,)
        Initial state. For problems in the complex domain, pass `y0` with a
        complex data type (even if the initial value is purely real).
    method : string or `OdeSolver`, optional
        Integration method to use:
    
            * 'RK45' (default): Explicit Runge-Kutta method of order 5(4) [1]_.
              The error is controlled assuming accuracy of the fourth-order
              method, but steps are taken using the fifth-order accurate
              formula (local extrapolation is done). A quartic interpolation
              polynomial is used for the dense output [2]_. Can be applied in
              the complex domain.
            * 'RK23': Explicit Runge-Kutta method of order 3(2) [3]_. The error
              is controlled assuming accuracy of the second-order method, but
              steps are taken using the third-order accurate formula (local
              extrapolation is done). A cubic Hermite polynomial is used for the
              dense output. Can be applied in the complex domain.
            * 'DOP853': Explicit Runge-Kutta method of order 8 [13]_.
              Python implementation of the "DOP853" algorithm originally
              written in Fortran [14]_. A 7-th order interpolation polynomial
              accurate to 7-th order is used for the dense output.
              Can be applied in the complex domain.
            * 'Radau': Implicit Runge-Kutta method of the Radau IIA family of
              order 5 [4]_. The error is controlled with a third-order accurate
              embedded formula. A cubic polynomial which satisfies the
              collocation conditions is used for the dense output.
            * 'BDF': Implicit multi-step variable-order (1 to 5) method based
              on a backward differentiation formula for the derivative
              approximation [5]_. The implementation follows the one described
              in [6]_. A quasi-constant step scheme is used and accuracy is
              enhanced using the NDF modification. Can be applied in the
              complex domain.
            * 'LSODA': Adams/BDF method with automatic stiffness detection and
              switching [7]_, [8]_. This is a wrapper of the Fortran solver
              from ODEPACK.
    
        Explicit Runge-Kutta methods ('RK23', 'RK45', 'DOP853') should be used
        for non-stiff problems and implicit methods ('Radau', 'BDF') for
        stiff problems [9]_. Among Runge-Kutta methods, 'DOP853' is recommended
        for solving with high precision (low values of `rtol` and `atol`).
    
        If not sure, first try to run 'RK45'. If it makes unusually many
        iterations, diverges, or fails, your problem is likely to be stiff and
        you should use 'Radau' or 'BDF'. 'LSODA' can also be a good universal
        choice, but it might be somewhat less convenient to work with as it
        wraps old Fortran code.
    
        You can also pass an arbitrary class derived from `OdeSolver` which
        implements the solver.
    t_eval : array_like or None, optional
        Times at which to store the computed solution, must be sorted and lie
        within `t_span`. If None (default), use points selected by the solver.
    dense_output : bool, optional
        Whether to compute a continuous solution. Default is False.
    events : callable, or list of callables, optional
        Events to track. If None (default), no events will be tracked.
        Each event occurs at the zeros of a continuous function of time and
        state. Each function must have the signature ``event(t, y)`` and return
        a float. The solver will find an accurate value of `t` at which
        ``event(t, y(t)) = 0`` using a root-finding algorithm. By default, all
        zeros will be found. The solver looks for a sign change over each step,
        so if multiple zero crossings occur within one step, events may be
        missed. Additionally each `event` function might have the following
        attributes:
    
            terminal: bool, optional
                Whether to terminate integration if this event occurs.
                Implicitly False if not assigned.
            direction: float, optional
                Direction of a zero crossing. If `direction` is positive,
                `event` will only trigger when going from negative to positive,
                and vice versa if `direction` is negative. If 0, then either
                direction will trigger event. Implicitly 0 if not assigned.
    
        You can assign attributes like ``event.terminal = True`` to any
        function in Python.
    vectorized : bool, optional
        Whether `fun` is implemented in a vectorized fashion. Default is False.
    args : tuple, optional
        Additional arguments to pass to the user-defined functions.  If given,
        the additional arguments are passed to all user-defined functions.
        So if, for example, `fun` has the signature ``fun(t, y, a, b, c)``,
        then `jac` (if given) and any event functions must have the same
        signature, and `args` must be a tuple of length 3.
    **options
        Options passed to a chosen solver. All options available for already
        implemented solvers are listed below.
    first_step : float or None, optional
        Initial step size. Default is `None` which means that the algorithm
        should choose.
    max_step : float, optional
        Maximum allowed step size. Default is np.inf, i.e., the step size is not
        bounded and determined solely by the solver.
    rtol, atol : float or array_like, optional
        Relative and absolute tolerances. The solver keeps the local error
        estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
        relative accuracy (number of correct digits), while `atol` controls
        absolute accuracy (number of correct decimal places). To achieve the
        desired `rtol`, set `atol` to be smaller than the smallest value that
        can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
        allowable error. If `atol` is larger than ``rtol * abs(y)`` the
        number of correct digits is not guaranteed. Conversely, to achieve the
        desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
        than `atol`. If components of y have different scales, it might be
        beneficial to set different `atol` values for different components by
        passing array_like with shape (n,) for `atol`. Default values are
        1e-3 for `rtol` and 1e-6 for `atol`.
    jac : array_like, sparse_matrix, callable or None, optional
        Jacobian matrix of the right-hand side of the system with respect
        to y, required by the 'Radau', 'BDF' and 'LSODA' method. The
        Jacobian matrix has shape (n, n) and its element (i, j) is equal to
        ``d f_i / d y_j``.  There are three ways to define the Jacobian:
    
            * If array_like or sparse_matrix, the Jacobian is assumed to
              be constant. Not supported by 'LSODA'.
            * If callable, the Jacobian is assumed to depend on both
              t and y; it will be called as ``jac(t, y)``, as necessary.
              For 'Radau' and 'BDF' methods, the return value might be a
              sparse matrix.
            * If None (default), the Jacobian will be approximated by
              finite differences.
    
        It is generally recommended to provide the Jacobian rather than
        relying on a finite-difference approximation.
    jac_sparsity : array_like, sparse matrix or None, optional
        Defines a sparsity structure of the Jacobian matrix for a finite-
        difference approximation. Its shape must be (n, n). This argument
        is ignored if `jac` is not `None`. If the Jacobian has only few
        non-zero elements in *each* row, providing the sparsity structure
        will greatly speed up the computations [10]_. A zero entry means that
        a corresponding element in the Jacobian is always zero. If None
        (default), the Jacobian is assumed to be dense.
        Not supported by 'LSODA', see `lband` and `uband` instead.
    lband, uband : int or None, optional
        Parameters defining the bandwidth of the Jacobian for the 'LSODA'
        method, i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``.
        Default is None. Setting these requires your jac routine to return the
        Jacobian in the packed format: the returned array must have ``n``
        columns and ``uband + lband + 1`` rows in which Jacobian diagonals are
        written. Specifically ``jac_packed[uband + i - j , j] = jac[i, j]``.
        The same format is used in `scipy.linalg.solve_banded` (check for an
        illustration).  These parameters can be also used with ``jac=None`` to
        reduce the number of Jacobian elements estimated by finite differences.
    min_step : float, optional
        The minimum allowed step size for 'LSODA' method.
        By default `min_step` is zero.
    
    Returns
    -------
    Bunch object with the following fields defined:
    t : ndarray, shape (n_points,)
        Time points.
    y : ndarray, shape (n, n_points)
        Values of the solution at `t`.
    sol : `OdeSolution` or None
        Found solution as `OdeSolution` instance; None if `dense_output` was
        set to False.
    t_events : list of ndarray or None
        Contains for each event type a list of arrays at which an event of
        that type event was detected. None if `events` was None.
    y_events : list of ndarray or None
        For each value of `t_events`, the corresponding value of the solution.
        None if `events` was None.
    nfev : int
        Number of evaluations of the right-hand side.
    njev : int
        Number of evaluations of the Jacobian.
    nlu : int
        Number of LU decompositions.
    status : int
        Reason for algorithm termination:
    
            * -1: Integration step failed.
            *  0: The solver successfully reached the end of `tspan`.
            *  1: A termination event occurred.
    
    message : string
        Human-readable description of the termination reason.
    success : bool
        True if the solver reached the interval end or a termination event
        occurred (``status >= 0``).
    
    References
    ----------
    .. [1] J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta
           formulae", Journal of Computational and Applied Mathematics, Vol. 6,
           No. 1, pp. 19-26, 1980.
    .. [2] L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics
           of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
    .. [3] P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
           Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
    .. [4] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:
           Stiff and Differential-Algebraic Problems", Sec. IV.8.
    .. [5] `Backward Differentiation Formula
            <https://en.wikipedia.org/wiki/Backward_differentiation_formula>`_
            on Wikipedia.
    .. [6] L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.
           COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
    .. [7] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
           Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
           pp. 55-64, 1983.
    .. [8] L. Petzold, "Automatic selection of methods for solving stiff and
           nonstiff systems of ordinary differential equations", SIAM Journal
           on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
           1983.
    .. [9] `Stiff equation <https://en.wikipedia.org/wiki/Stiff_equation>`_ on
           Wikipedia.
    .. [10] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
            sparse Jacobian matrices", Journal of the Institute of Mathematics
            and its Applications, 13, pp. 117-120, 1974.
    .. [11] `Cauchy-Riemann equations
             <https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>`_ on
             Wikipedia.
    .. [12] `Lotka-Volterra equations
            <https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations>`_
            on Wikipedia.
    .. [13] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
            Equations I: Nonstiff Problems", Sec. II.
    .. [14] `Page with original Fortran code of DOP853
            <http://www.unige.ch/~hairer/software.html>`_.
    
    Examples
    --------
    Basic exponential decay showing automatically chosen time points.
    
    >>> from scipy.integrate import solve_ivp
    >>> def exponential_decay(t, y): return -0.5 * y
    >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])
    >>> print(sol.t)
    [ 0.          0.11487653  1.26364188  3.06061781  4.81611105  6.57445806
      8.33328988 10.        ]
    >>> print(sol.y)
    [[2.         1.88836035 1.06327177 0.43319312 0.18017253 0.07483045
      0.03107158 0.01350781]
     [4.         3.7767207  2.12654355 0.86638624 0.36034507 0.14966091
      0.06214316 0.02701561]
     [8.         7.5534414  4.25308709 1.73277247 0.72069014 0.29932181
      0.12428631 0.05403123]]
    
    Specifying points where the solution is desired.
    
    >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8],
    ...                 t_eval=[0, 1, 2, 4, 10])
    >>> print(sol.t)
    [ 0  1  2  4 10]
    >>> print(sol.y)
    [[2.         1.21305369 0.73534021 0.27066736 0.01350938]
     [4.         2.42610739 1.47068043 0.54133472 0.02701876]
     [8.         4.85221478 2.94136085 1.08266944 0.05403753]]
    
    Cannon fired upward with terminal event upon impact. The ``terminal`` and
    ``direction`` fields of an event are applied by monkey patching a function.
    Here ``y[0]`` is position and ``y[1]`` is velocity. The projectile starts
    at position 0 with velocity +10. Note that the integration never reaches
    t=100 because the event is terminal.
    
    >>> def upward_cannon(t, y): return [y[1], -0.5]
    >>> def hit_ground(t, y): return y[0]
    >>> hit_ground.terminal = True
    >>> hit_ground.direction = -1
    >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
    >>> print(sol.t_events)
    [array([40.])]
    >>> print(sol.t)
    [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
     1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
    
    Use `dense_output` and `events` to find position, which is 100, at the apex
    of the cannonball's trajectory. Apex is not defined as terminal, so both
    apex and hit_ground are found. There is no information at t=20, so the sol
    attribute is used to evaluate the solution. The sol attribute is returned
    by setting ``dense_output=True``. Alternatively, the `y_events` attribute
    can be used to access the solution at the time of the event.
    
    >>> def apex(t, y): return y[1]
    >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10],
    ...                 events=(hit_ground, apex), dense_output=True)
    >>> print(sol.t_events)
    [array([40.]), array([20.])]
    >>> print(sol.t)
    [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
     1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
    >>> print(sol.sol(sol.t_events[1][0]))
    [100.   0.]
    >>> print(sol.y_events)
    [array([[-5.68434189e-14, -1.00000000e+01]]), array([[1.00000000e+02, 1.77635684e-15]])]
    
    As an example of a system with additional parameters, we'll implement
    the Lotka-Volterra equations [12]_.
    
    >>> def lotkavolterra(t, z, a, b, c, d):
    ...     x, y = z
    ...     return [a*x - b*x*y, -c*y + d*x*y]
    ...
    
    We pass in the parameter values a=1.5, b=1, c=3 and d=1 with the `args`
    argument.
    
    >>> sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1),
    ...                 dense_output=True)
    
    Compute a dense solution and plot it.
    
    >>> t = np.linspace(0, 15, 300)
    >>> z = sol.sol(t)
    >>> import matplotlib.pyplot as plt
    >>> plt.plot(t, z.T)
    >>> plt.xlabel('t')
    >>> plt.legend(['x', 'y'], shadow=True)
    >>> plt.title('Lotka-Volterra System')
    >>> plt.show()

In [17]:

tf=100.
sol=solve_ivp(anHarmonic, [0,tf], [x,v], rtol=1e-14, atol=5e-7)
print(sol)

  message: 'The solver successfully reached the end of the integration interval.'
     nfev: 2630
     njev: 0
      nlu: 0
      sol: None
   status: 0
  success: True
        t: array([0.00000000e+00, 2.34367291e-02, 2.32365158e-01, 4.47372590e-01,
       6.87883729e-01, 9.61928130e-01, 1.26831360e+00, 1.59744802e+00,
       1.96014689e+00, 2.35100741e+00, 2.73837360e+00, 3.12707393e+00,
       3.51178920e+00, 3.87096546e+00, 4.20505645e+00, 4.54962871e+00,
       4.78908650e+00, 5.00018945e+00, 5.21129240e+00, 5.40503334e+00,
       5.58675814e+00, 5.76098202e+00, 5.93054034e+00, 6.09727981e+00,
       6.26240404e+00, 6.42670351e+00, 6.59080280e+00, 6.75541163e+00,
       6.92150115e+00, 7.09038495e+00, 7.26376557e+00, 7.44382418e+00,
       7.63341909e+00, 7.83647041e+00, 8.05863759e+00, 8.30811077e+00,
       8.59284643e+00, 8.90607918e+00, 9.24346658e+00, 9.61767274e+00,
       1.00087745e+01, 1.03957536e+01, 1.07849143e+01, 1.11642218e+01,
       1.15142969e+01, 1.18459430e+01, 1.21422326e+01, 1.23752485e+01,
       1.25822851e+01, 1.27893218e+01, 1.29805028e+01, 1.31606602e+01,
       1.33339061e+01, 1.35028646e+01, 1.36692523e+01, 1.38341849e+01,
       1.39984056e+01, 1.41625367e+01, 1.43273237e+01, 1.44937933e+01,
       1.46633229e+01, 1.48376894e+01, 1.50191730e+01, 1.52107819e+01,
       1.54166794e+01, 1.56429183e+01, 1.58981004e+01, 1.61888364e+01,
       1.65057985e+01, 1.68487205e+01, 1.72290004e+01, 1.76192308e+01,
       1.80064171e+01, 1.83953852e+01, 1.87704807e+01, 1.91156414e+01,
       1.94475159e+01, 1.97249054e+01, 1.99528431e+01, 2.01568430e+01,
       2.03608430e+01, 2.05500833e+01, 2.07290416e+01, 2.09015393e+01,
       2.10700415e+01, 2.12361639e+01, 2.14009549e+01, 2.15651242e+01,
       2.17292981e+01, 2.18942551e+01, 2.20610693e+01, 2.22311717e+01,
       2.24064026e+01, 2.25891289e+01, 2.27824918e+01, 2.29908719e+01,
       2.32206672e+01, 2.34807557e+01, 2.37761563e+01, 2.40961914e+01,
       2.44440623e+01, 2.48285572e+01, 2.52178865e+01, 2.56054149e+01,
       2.59938479e+01, 2.63649919e+01, 2.67064977e+01, 2.70397294e+01,
       2.73041712e+01, 2.75277137e+01, 2.77512561e+01, 2.79526153e+01,
       2.81388903e+01, 2.83159402e+01, 2.84872467e+01, 2.86550280e+01,
       2.88207364e+01, 2.89853130e+01, 2.91494223e+01, 2.93137092e+01,
       2.94790169e+01, 2.96465047e+01, 2.98176998e+01, 2.99945592e+01,
       3.01796138e+01, 3.03762599e+01, 3.05892980e+01, 3.08257683e+01,
       3.10947763e+01, 3.13975395e+01, 3.17232454e+01, 3.20804902e+01,
       3.24699831e+01, 3.28578891e+01, 3.32461855e+01, 3.36326475e+01,
       3.39961956e+01, 3.43324483e+01, 3.46713857e+01, 3.49188239e+01,
       3.51343365e+01, 3.53498491e+01, 3.55463294e+01, 3.57297119e+01,
       3.59049684e+01, 3.60751634e+01, 3.62422798e+01, 3.64076130e+01,
       3.65720067e+01, 3.67360952e+01, 3.69005586e+01, 3.70663161e+01,
       3.72346202e+01, 3.74071036e+01, 3.75858566e+01, 3.77736025e+01,
       3.79740439e+01, 3.81924849e+01, 3.84366742e+01, 3.87153342e+01,
       3.90247652e+01, 3.93572576e+01, 3.97249564e+01, 4.01162944e+01,
       4.05033707e+01, 4.08923596e+01, 4.12750462e+01, 4.16302593e+01,
       4.19629170e+01, 4.23147985e+01, 4.25473133e+01, 4.27544459e+01,
       4.29615784e+01, 4.31528127e+01, 4.33330025e+01, 4.35062686e+01,
       4.36752396e+01, 4.38416345e+01, 4.40065712e+01, 4.41707934e+01,
       4.43349237e+01, 4.44997067e+01, 4.46661679e+01, 4.48356832e+01,
       4.50100280e+01, 4.51914802e+01, 4.53830448e+01, 4.55888796e+01,
       4.58150287e+01, 4.60700858e+01, 4.63606971e+01, 4.66775802e+01,
       4.70203784e+01, 4.74005383e+01, 4.77907922e+01, 4.81779728e+01,
       4.85669509e+01, 4.89421434e+01, 4.92874035e+01, 4.96192583e+01,
       4.98970113e+01, 5.01250597e+01, 5.03291260e+01, 5.05331923e+01,
       5.07224752e+01, 5.09014600e+01, 5.10739743e+01, 5.12424866e+01,
       5.14086149e+01, 5.15734091e+01, 5.17375795e+01, 5.19017523e+01,
       5.20667054e+01, 5.22335116e+01, 5.24036006e+01, 5.25788114e+01,
       5.27615090e+01, 5.29548312e+01, 5.31631537e+01, 5.33928663e+01,
       5.36528417e+01, 5.39481401e+01, 5.42681050e+01, 5.46158600e+01,
       5.50002696e+01, 5.53896209e+01, 5.57771417e+01, 5.61655921e+01,
       5.65368292e+01, 5.68784126e+01, 5.72116009e+01, 5.74763049e+01,
       5.76999479e+01, 5.79235909e+01, 5.81250108e+01, 5.83113212e+01,
       5.84883931e+01, 5.86597133e+01, 5.88275029e+01, 5.89932160e+01,
       5.91577951e+01, 5.93219049e+01, 5.94861903e+01, 5.96514935e+01,
       5.98189727e+01, 5.99901539e+01, 6.01669928e+01, 6.03520181e+01,
       6.05486230e+01, 6.07616024e+01, 6.09979886e+01, 6.12668870e+01,
       6.15695680e+01, 6.18952020e+01, 6.22523294e+01, 6.26417821e+01,
       6.30297036e+01, 6.34179920e+01, 6.38044872e+01, 6.41681294e+01,
       6.45044353e+01, 6.48432742e+01, 6.50908899e+01, 6.53064983e+01,
       6.55221067e+01, 6.57186460e+01, 6.59020639e+01, 6.60773426e+01,
       6.62475513e+01, 6.64146759e+01, 6.65800137e+01, 6.67444097e+01,
       6.69084983e+01, 6.70729590e+01, 6.72387102e+01, 6.74070030e+01,
       6.75794688e+01, 6.77581960e+01, 6.79459053e+01, 6.81462950e+01,
       6.83646623e+01, 6.86087465e+01, 6.88872812e+01, 6.91966318e+01,
       6.95290292e+01, 6.98965919e+01, 7.02879263e+01, 7.06750094e+01,
       7.10639933e+01, 7.14467412e+01, 7.18020610e+01, 7.21347505e+01,
       7.24864018e+01, 7.27191012e+01, 7.29263408e+01, 7.31335805e+01,
       7.33248830e+01, 7.35051147e+01, 7.36784071e+01, 7.38473941e+01,
       7.40137985e+01, 7.41787403e+01, 7.43429646e+01, 7.45070939e+01,
       7.46718716e+01, 7.48383218e+01, 7.50078187e+01, 7.51821353e+01,
       7.53635470e+01, 7.55550543e+01, 7.57608080e+01, 7.59868409e+01,
       7.62417362e+01, 7.65321857e+01, 7.68489663e+01, 7.71916045e+01,
       7.75716081e+01, 7.79618920e+01, 7.83490649e+01, 7.87380553e+01,
       7.91133731e+01, 7.94587626e+01, 7.97905929e+01, 8.00688202e+01,
       8.02970121e+01, 8.05011644e+01, 8.07053166e+01, 8.08946549e+01,
       8.10736739e+01, 8.12462096e+01, 8.14147350e+01, 8.15808709e+01,
       8.17456692e+01, 8.19098410e+01, 8.20740125e+01, 8.22389604e+01,
       8.24057561e+01, 8.25758280e+01, 8.27510128e+01, 8.29336731e+01,
       8.31269427e+01, 8.33351906e+01, 8.35647963e+01, 8.38246255e+01,
       8.41197912e+01, 8.44396652e+01, 8.47872704e+01, 8.51715687e+01,
       8.55609481e+01, 8.59484587e+01, 8.63369311e+01, 8.67082887e+01,
       8.70499730e+01, 8.73831059e+01, 8.76481522e+01, 8.78719257e+01,
       8.80956992e+01, 8.82971977e+01, 8.84835541e+01, 8.86606543e+01,
       8.88319923e+01, 8.89997926e+01, 8.91655118e+01, 8.93300941e+01,
       8.94942047e+01, 8.96584880e+01, 8.98237852e+01, 8.99912533e+01,
       9.01624166e+01, 9.03392291e+01, 9.05242165e+01, 9.07207679e+01,
       9.09336713e+01, 9.11699488e+01, 9.14387055e+01, 9.17412796e+01,
       9.20668206e+01, 9.24237960e+01, 9.28131960e+01, 9.32011373e+01,
       9.35894148e+01, 9.39759524e+01, 9.43397167e+01, 9.46760918e+01,
       9.50148041e+01, 9.52626508e+01, 9.54783836e+01, 9.56941164e+01,
       9.58907323e+01, 9.60741963e+01, 9.62495037e+01, 9.64197302e+01,
       9.65868654e+01, 9.67522091e+01, 9.69166079e+01, 9.70806966e+01,
       9.72451540e+01, 9.74108970e+01, 9.75791752e+01, 9.77516182e+01,
       9.79303120e+01, 9.81179738e+01, 9.83182966e+01, 9.85365684e+01,
       9.87805167e+01, 9.90588886e+01, 9.93681347e+01, 9.97004092e+01,
       1.00000000e+02])
 t_events: None
        y: array([[ 0.        ,  0.0234346 ,  0.23040018,  0.43421579,  0.64358859,
         0.85146917,  1.04064699,  1.1919187 ,  1.29835755,  1.34590673,
         1.32639948,  1.23920694,  1.08277299,  0.87074982,  0.61711532,
         0.30672984,  0.07183174, -0.13874558, -0.34257795, -0.51414781,
        -0.65322639, -0.75985989, -0.83323102, -0.87243946, -0.87708576,
        -0.84770695, -0.78596996, -0.69453041, -0.57665919, -0.43586027,
        -0.27561643, -0.09926548,  0.09004217,  0.28933311,  0.49568653,
         0.70549494,  0.91036177,  1.08955243,  1.2283527 ,  1.31932021,
         1.34702514,  1.30794466,  1.20009623,  1.02513492,  0.80029204,
         0.53359558,  0.25920841,  0.02898161, -0.17706996, -0.37474747,
        -0.54057121, -0.67405497, -0.77499128, -0.84249673, -0.8757201 ,
        -0.87439798, -0.83925842, -0.77215894, -0.67589166, -0.55379257,
        -0.40937599, -0.24610554, -0.06728795,  0.12395933,  0.32468702,
         0.53194444,  0.74177072,  0.9439573 ,  1.11665493,  1.24790832,
         1.32901573,  1.34499562,  1.29427936,  1.17422735,  0.98916381,
         0.75717995,  0.48190421,  0.22084951, -0.00535607, -0.20775043,
        -0.40050433, -0.56163882, -0.69050301, -0.78671804, -0.84936405,
        -0.87764293, -0.87141291, -0.83156301, -0.76010336, -0.65992941,
        -0.53442233, -0.38709984, -0.2214074 , -0.04062546,  0.15215542,
         0.35400288,  0.56192429,  0.77156435,  0.97096813,  1.13807714,
         1.26294293,  1.33547292,  1.34182618,  1.28141129,  1.15134128,
         0.95843429,  0.72066485,  0.43754239,  0.18553572, -0.03696767,
        -0.25739168, -0.44380377, -0.5969029 , -0.71766915, -0.80557749,
        -0.85968267, -0.87930557, -0.86453589, -0.81654833, -0.73760709,
        -0.63076769, -0.49947306, -0.34723703, -0.17746909,  0.00659739,
         0.20191748,  0.40557943,  0.61446413,  0.8232228 ,  1.01647458,
         1.17347428,  1.28669417,  1.34337088,  1.33303863,  1.25545265,
         1.10820996,  0.90265477,  0.65476135,  0.3554111 ,  0.11491439,
        -0.10017819, -0.31022924, -0.48748144, -0.63200971, -0.74416487,
        -0.82321985, -0.86824315, -0.87871578, -0.8549946 , -0.79855663,
        -0.71191317, -0.59826039, -0.46108483, -0.30388455, -0.13002643,
         0.05730678,  0.25511534,  0.46049094,  0.67007047,  0.87688225,
         1.06198398,  1.20797115,  1.30798132,  1.34708728,  1.31922967,
         1.22323208,  1.05868387,  0.84104266,  0.58199883,  0.25962582,
         0.02991358, -0.17624707, -0.37406709, -0.54001383, -0.6736175 ,
        -0.77467618, -0.84230762, -0.8756594 , -0.87446503, -0.83944842,
        -0.77246307, -0.6762983 , -0.55428882, -0.40994881, -0.24674232,
        -0.06797679,  0.12322966,  0.32392725,  0.53116616,  0.74099439,
         0.94324588,  1.11608577,  1.24750229,  1.32882636,  1.34505584,
         1.29458841,  1.17479511,  0.98993917,  0.75810501,  0.48302073,
         0.22171099, -0.00458467, -0.20706142, -0.39992582, -0.56116634,
        -0.69013549, -0.78645795, -0.84921454, -0.87760581, -0.87148712,
        -0.83174382, -0.76038264, -0.66029676, -0.5348664 , -0.38760926,
        -0.22197125, -0.04123339,  0.15151312,  0.35333558,  0.56124251,
         0.77088881,  0.97036138,  1.13759866,  1.26261074,  1.3353389 ,
         1.34190882,  1.28171322,  1.15186712,  0.95913122,  0.72149124,
         0.43855421,  0.1863668 , -0.0362235 , -0.25677079, -0.44328975,
        -0.59648716, -0.71735155, -0.8053607 , -0.85956946, -0.87929707,
        -0.86463032, -0.81674041, -0.7378886 , -0.63112862, -0.49990278,
        -0.34772502, -0.17800531,  0.00602241,  0.20131267,  0.40495355,
         0.61382798,  0.82260198,  1.01593732,  1.17306066,  1.28642496,
         1.34329777,  1.3331648 ,  1.25578574,  1.10874424,  0.90333248,
         0.65556134,  0.35642953,  0.11581206, -0.09937433, -0.30955545,
        -0.48692495, -0.63156475, -0.7438326 , -0.82300357, -0.86814571,
        -0.87873769, -0.85513278, -0.79880407, -0.71225979, -0.59869447,
        -0.46159422, -0.30445735, -0.13065133,  0.05664041,  0.25441758,
         0.45977197,  0.66934466,  0.87618969,  1.06140682,  1.20753904,
         1.30772801,  1.34706532,  1.31943499,  1.22367684,  1.05934481,
         0.84185282,  0.58295847,  0.26093989,  0.03108556, -0.1751989 ,
        -0.3731866 , -0.53929177, -0.6730507 , -0.77426782, -0.84206241,
        -0.87558046, -0.87455153, -0.83969416, -0.77285665, -0.67682467,
        -0.55493127, -0.41069045, -0.24756683, -0.06886871,  0.12228485,
         0.32294347,  0.5301584 ,  0.73998907,  0.94232419,  1.11534823,
         1.24697607,  1.32858057,  1.34513368,  1.29498896,  1.17553148,
         0.99094546,  0.75930568,  0.48446927,  0.2228266 , -0.0035858 ,
        -0.20616929, -0.39917682, -0.56055463, -0.68965964, -0.78612116,
        -0.84902086, -0.87755759, -0.87158302, -0.83197773, -0.76074402,
        -0.66077215, -0.53544112, -0.3882686 , -0.22270103, -0.04202025,
         0.15068179,  0.35247191,  0.56036013,  0.77001441,  0.96957572,
         1.13697906,  1.26218056,  1.33516524,  1.34201604,  1.28210482,
         1.1525494 ,  0.96003588,  0.72256392,  0.43986688,  0.18744305,
        -0.03525988, -0.25596692, -0.44262433, -0.59594897, -0.7169404 ,
        -0.80508003, -0.85942285, -0.87928599, -0.86475245, -0.81698896,
        -0.73825291, -0.63159574, -0.50045896, -0.34835664, -0.17869935,
         0.00527822,  0.20052989,  0.40414354,  0.6130047 ,  0.82179842,
         1.01524169,  1.17252516,  1.28607642,  1.34320319,  1.33332862,
         1.25621799,  1.10943781,  0.90421249,  0.65659998,  0.35775074,
         0.11697619, -0.09833194, -0.30868187, -0.48620348, -0.63098785,
        -0.74340176, -0.82272305, -0.86801918, -0.87876586, -0.85531165,
        -0.79912457, -0.71270887, -0.59925694, -0.4622543 , -0.30519963,
        -0.13146116,  0.05577684,  0.25351334,  0.45884027,  0.66840404,
         0.87529192,  1.0606584 ,  1.20697872,  1.29356321],
       [ 1.        ,  0.99972752,  0.97518844,  0.91583072,  0.82137255,
         0.69337687,  0.54102523,  0.37917082,  0.20935221,  0.03496053,
        -0.13610495, -0.31414713, -0.50074846, -0.67960538, -0.83516902,
        -0.95683637, -0.99747854, -0.98987833, -0.93232963, -0.83087136,
        -0.69309704, -0.52570744, -0.33599942, -0.1322742 ,  0.07623757,
         0.27988644,  0.46947496,  0.63714429,  0.77688071,  0.88454759,
         0.95763033,  0.99489674,  0.99605967,  0.96143594,  0.89156679,
         0.78697061,  0.65017708,  0.49398284,  0.33007353,  0.15768663,
        -0.01540291, -0.18740732, -0.36865481, -0.55515993, -0.72827547,
        -0.87516691, -0.96882154, -0.99958345, -0.98325712, -0.91761685,
        -0.80939918, -0.66600544, -0.49419758, -0.3014353 , -0.09619808,
         0.11216176,  0.31400451,  0.50030786,  0.66351086,  0.79793411,
         0.8997449 ,  0.96667565,  0.99768188,  0.9926063 ,  0.95183307,
         0.87590356,  0.76538932,  0.62395134,  0.4660531 ,  0.30083491,
         0.12723074, -0.04497735, -0.21805221, -0.40108206, -0.58665759,
        -0.75589621, -0.89725738, -0.97714591, -0.99998467, -0.97665194,
        -0.9045332 , -0.79088931, -0.6430094 , -0.46772459, -0.27263448,
        -0.06636246,  0.14165211,  0.34179877,  0.52521964,  0.68460916,
         0.81456259,  0.91149121,  0.97332342,  0.99916202,  0.98894911,
         0.94311555,  0.86219633,  0.74685516,  0.60195858,  0.44286655,
         0.27649118,  0.10218824, -0.06941034, -0.24361372, -0.427954  ,
        -0.61226817, -0.77807478, -0.91459026, -0.98371877, -0.99930665,
        -0.96336013, -0.8797148 , -0.75684444, -0.60144555, -0.42045334,
        -0.22171805, -0.01410382,  0.19283073,  0.38957376,  0.56759515,
         0.7200521 ,  0.84201455,  0.93030502,  0.98317768,  0.99997663,
         0.98080071,  0.92614274,  0.83651889,  0.71292735,  0.5629001 ,
         0.40198941,  0.2334315 ,  0.05876747, -0.1122672 , -0.28884868,
        -0.47491271, -0.656024  , -0.81538462, -0.94267924, -0.99362816,
        -0.9947993 , -0.945411  , -0.85073125, -0.71861747, -0.55574   ,
        -0.36924608, -0.16726119,  0.04111735,  0.24625665,  0.43881542,
         0.61065959,  0.75545026,  0.86874705,  0.9477848 ,  0.99113888,
         0.99838595,  0.96976623,  0.90581072,  0.80702359,  0.67514775,
         0.52098388,  0.35826812,  0.18731013,  0.01339107, -0.15793375,
        -0.33734965, -0.52413012, -0.70067643, -0.85263414, -0.96872258,
        -0.99955476, -0.98341635, -0.91794459, -0.8098696 , -0.66659446,
        -0.49487924, -0.30218004, -0.09697255,  0.11139334,  0.31327744,
         0.4996534 ,  0.66295374,  0.79749193,  0.8994288 ,  0.96649161,
         0.99763197,  0.99268989,  0.95204822,  0.87624841,  0.76586027,
         0.62451675,  0.46665156,  0.30146226,  0.12787993, -0.04434694,
        -0.21739651, -0.40039117, -0.58599367, -0.75531787, -0.89679991,
        -0.97697072, -0.99998688, -0.97681175, -0.90483853, -0.79131666,
        -0.64353739, -0.46833013, -0.27329122, -0.06704083,  0.14098352,
         0.34117051,  0.52465832,  0.6841355 ,  0.81419112,  0.91123096,
         0.97317907,  0.99913497,  0.98903859,  0.94331978,  0.86251376,
         0.74728145,  0.60245921,  0.4433925 ,  0.27704397,  0.10275267,
        -0.06885889, -0.24303577, -0.42734951, -0.61169706, -0.77758279,
        -0.91420992, -0.98357461, -0.99933277, -0.96354541, -0.88002952,
        -0.75726773, -0.60195743, -0.42103162, -0.22233739, -0.01473609,
         0.19221487,  0.38900208,  0.56709118,  0.71963367,  0.84169382,
         0.9300893 ,  0.98307046,  0.99997862,  0.98091083,  0.92635933,
         0.8368405 ,  0.71334668,  0.56337411,  0.40248407,  0.23395321,
         0.05928609, -0.1117519 , -0.2883036 , -0.47435305, -0.65550964,
        -0.81495028, -0.94236136, -0.9935283 , -0.99488231, -0.94566432,
        -0.85112557, -0.71912979, -0.55634703, -0.36992158, -0.1679753 ,
         0.04039736,  0.24556413,  0.43818103,  0.61010856,  0.75500117,
         0.8684122 ,  0.94757127,  0.9910497 ,  0.99842126,  0.96992449,
         0.90609007,  0.8074221 ,  0.67565023,  0.52153226,  0.35884013,
         0.18791254,  0.01397748, -0.15733857, -0.33671793, -0.5234983 ,
        -0.70011028, -0.85216745, -0.96841372, -0.99951791, -0.98361842,
        -0.91836805, -0.81047833, -0.66735687, -0.49576169, -0.30314427,
        -0.09797539,  0.11039823,  0.31233576,  0.49880566,  0.66223201,
         0.79691899,  0.89901914,  0.96625299,  0.99756707,  0.99279791,
         0.95232663,  0.87669481,  0.76647002,  0.62524921,  0.46742708,
         0.30227517,  0.12872155, -0.04352945, -0.21654613, -0.39949474,
        -0.58513172, -0.7545668 , -0.89620552, -0.97674342, -0.99998933,
        -0.97701826, -0.90523345, -0.79186954, -0.64422059, -0.46911374,
        -0.27414116, -0.06791884,  0.14011811,  0.34035722,  0.52393162,
         0.68352226,  0.81371013,  0.91089396,  0.97299211,  0.99909987,
         0.98915434,  0.9435841 ,  0.86292461,  0.74783328,  0.60310761,
         0.44407382,  0.27776007,  0.10348426, -0.06814395, -0.24228628,
        -0.4265652 , -0.61095568, -0.77694397, -0.91371585, -0.98338753,
        -0.99936622, -0.96378498, -0.88043669, -0.75781542, -0.60261981,
        -0.42177997, -0.22313893, -0.0155544 ,  0.19141774,  0.38826209,
         0.56643879,  0.719092  ,  0.84127861,  0.92981002,  0.98293165,
         0.99998118,  0.98105339,  0.92663971,  0.83725683,  0.71388962,
         0.56398811,  0.40312489,  0.23462911,  0.05995841, -0.11108373,
        -0.28759657, -0.47362669, -0.65484179, -0.81438621, -0.94194828,
        -0.9933982 , -0.99498935, -0.94599223, -0.85163628, -0.71979348,
        -0.55713352, -0.37079687, -0.1689007 ,  0.03946425,  0.24466651,
         0.43735868,  0.6093942 ,  0.75441891,  0.86797801,  0.94729433,
         0.99093393,  0.99846688,  0.97012949,  0.90645203,  0.80793853,
         0.67630161,  0.52224349,  0.35958202,  0.21952698]])
 y_events: None

/Users/haule/opt/anaconda3/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:47: UserWarning: At least one element of `rtol` is too small. Setting `rtol = np.maximum(rtol, 2.220446049250313e-14)`.
  warn("At least one element of `rtol` is too small. "

We choose absolute tolerance atol=5e-7 so that it uses approximately the same function evaluations, namely, around 2000.

In [18]:

plt.plot(sol.t, sol.y[0], label='x:scipy.solve_ivp')
plt.plot(t, data2R[:,0], label='x:RK4')
plt.legend(loc='best');

Solution is very similar to our RK4. What about total energy and stability?

In [19]:

Enr = E_anharmonic(sol.y)

plt.plot(sol.t, Enr[0]-E0)
print(Enr[0,-1]-E0)

-5.527636871005548e-06

Dissapointingly, the error is even a bit larger, and definitely not stable.

Verlet algorithm¶

Verlet algorithm is used in molecular dynamics simulations (when there is no energy loss), because it is stable and does not lead to energy loss with time. The precision of Verlet algorithm is only $O(h^3)$, which is worse than RK4, but stability is here more important than precision. Below we derive Verlet's algorithm.

If there is no friction, we can explore the time-symmetry of Newton's Eq: they are time invariant when there is no friction.

First we recall the Taylor expansion of the soltuion, which should exist:

$$r(t+h)=r(t)+\dot{r}(t) h + \frac{1}{2}\ddot{r}(t) h^2 + O(h^3)$$

and because of the symmetry of Newton's Eq, we also have:

$$r(t-h)=r(t)-\dot{r}(t) h + \frac{1}{2}\ddot{r}(t) h^2 - O(h^3)$$

Next we assume that the explicit form for the acceleration $\ddot{r}(t)=a(t)$ exists in the form of the Newton's Eq. We will also replace $\dot{r}(t)$ with symbol for velocity $v$.

We will now rewrite the last equation, but shift the time variable from $t\rightarrow t+h$: $$r(t)=r(t+h)-v(t+h) h + \frac{1}{2}a(t+h) h^2 + O(h^3)$$ The very first Eq. above can be simply rewritten with $v$ and $a$ as: $$r(t+h)=r(t)+ v(t) h + \frac{1}{2}a(t) h^2 + O(h^3)$$

If we sum these two Eq. together we get: $$ v(t+h) = v(t) + \frac{1}{2}(a(t)+a(t+h)) h + O(h^3)$$

The velocity Verlet algorithm uses the last two Eq. (for $r(t+h)$ and $v(t+h)$).

We notice that both Equations contain a common term, which we will call $v(t+h/2)$, namely, $$v(t+\frac{1}{2}h) \equiv v(t) + \frac{h}{2} a(t)$$ In terms of this quantity, we can rewrite the above two Equations as \begin{eqnarray} && v(t+h) = v(t+\frac{1}{2}h) + \frac{h}{2}a(t+h) \\ && r(t+h) = r(t) + h\; v(t+\frac{1}{2}h) \end{eqnarray}

Below we implement velocity Verlet algorithm:

In [20]:

def velocity_verlet(y, f, t, h):
    """Velocity Verlet

    Low-performance implementation because the force is calculated
    twice; should remember the second force calculation and use as
    input for the next step.

    For comparing different algorithms it is ok to use this
    implementation for convenience. For production level code you
    should use a better implementation that saves the second force
    evaluation.

    """
    F = f(t, y)   # [dy/dt, a]
    # a = F[1]      # this is a(t)
    # half step velocity
    y[1] += 0.5*h*F[1] # y[1] <- v[t+h/2] = v[t] + h/2 * a
    # full step position
    y[0] += h*y[1]     # r[t+h] = r[t] + h * v[t+h/2]
    # full step velocity (updated positions!)
    F = f(t+h, y)   # this is a(t+h), which should be reused in state of the art implementation
    y[1] += 0.5*h*F[1] # v[t+h] = v[t+h/2] + h/2 * a
    return y

Now we can just use it with previously developed fixed step Solve routine:

In [21]:

data3R = Solve(t, [x,v], velocity_verlet, anHarmonic)
Env = E_anharmonic(data3R.T)

In [22]:

plt.plot(t,data3R[:,0], label='verlet')
plt.plot(t,data2R[:,0], label='RK4')
plt.legend(loc='best')

Out[22]:

<matplotlib.legend.Legend at 0x7fba981d09a0>

The verlet algorithm gives similar result as RK4 (but slightly less precise).

What about total energy?

In [23]:

plt.plot(t,Env[0]-E0)

Out[23]:

[<matplotlib.lines.Line2D at 0x7fbaa80870a0>]

The total energy is substantially less precise, however, it is stable, in the sense that no matter how long we evolve the system, the error of the total energy is bounded. It can be shown that the cummulative error of Verlet algorithm is $O((\Delta t)^2)$, where $\Delta t$ is time-step. While RK4 is locally $O(h^5)$ for each step, but cummulative error is unbounded, hence diverging. Only Verlet has finite cummulative error.

Homework: The Keppler problem¶

Simulate the motion of Earth in the solar system as a two body problem (taking into account only the Sun and the Earth). Derive the Newton's equation for the relative coordinate:

$$m \ddot{\vec{r}}=-G m M \frac{\vec{r}}{r}.$$

Turn the above equation into atronomical units (AU). The astronomical units are: - length is meassured in units of distance between Earth and Sun $R\approx 1.5\;10^{11}$m - time is meassured in years.

Plot the Earth's orbit $(x,y)$ for 5000 years, and verify it is stable (the orbit is a circle with no time dependence). Also plot $x(t)$ and $y(t)$ for the last 5 out of 5000 years to see that period is what is expected.

Simulate the three body problem, which consists of Sun, Earth and Jupiter. The latter has mass $m_J/m_S\approx 9.55\; 10^{-4}$, and distance $5.2\; AU$. Note that the mass of Earth is $m_E/m_S\approx 3 10^{-6}$ and distance $1\; AU$.

Simulate 5000 Earth years and plot the orbits of Sun, Eart, and Jupiter ($x(t)$ versus $y(t)$ for all three objects). Are the orbits stable?
Check how strong is the influence of the Jupiter on motion of the Earth. Plot $x(t)$ for the last 5 years of 5000 years for the case with and withouth Jupiter.

It turns out that our solar system has uneven distribution of asteroids, as demonstrated by the picture below. We are plotting the number of asteroids as a function of the distance from the Sun. There are many gaps in the distribution plot, which are now named Kirkwood gaps, after Daniel Kirkwood, who discovered them. He showed that gaps are associated with Jupiter, because the orbits are in resonance with Jupiter's motion. For example, the 2/1 gap is such that an asteroid placed there would complete two orbits every time Jupiter completes one. Similarly there are 3/1, 5/2, and 7/3 resonance, all related to Jupiter.

We would like to simulate Kirkwood gaps. To simulate 2/1 gap, we need to find the distance from the Sun that such asteroid would be placed at, and his initial velocity. Than we need to check the long term stability of such orbit.

We first recall that all orbits in the solar system satisfy $R_i^3/T_i^2=const$. This is because centrifugal force ($m_i\omega_i^2 R_i$) and gravitational force ($G M m_i/R_i^2$) have to be balanced, hence $m_i (2\pi/T_i)^2 R_i = G M m_i/R_i^2$, where $M$ is the Solar's mass. We hence see that $R_i^3/T_i^2=G M/(4\pi^2)$. In AU units this equation is $R_i^3/T_i^2=1$

For Asteroid that completes the orbit in half the Saturn's year, it must satisfy $T_{asteroid}=R_{Saturn}^{3/2}/2$, and hence $R_{asteroid}=R_{Saturn}/2^{2/3}$, which is $R_{asteroid}=5.2/2^{2/3}\approx 3.2758$. It's starting velocity should be $v_{asteroid}=2\pi/\sqrt{R_{asteroid}}\approx 3.4715$.

For 3/1 gap, we should similarly have $R_{asteroid}=5.2/3^{2/3}\approx 2.5$ and $v_{asteroid}=2\pi/\sqrt{R_{asteroid}}\approx 3.974$.

For homework, simulate the three body problem: Sun, Jupiter and Asteroid. You can assume that the mass of Asteroid is vanishingly small. Here it is also safe to ignore Earth and its influence on Asteroid. Simulate 5000 Earth years of motion, and plot the orbits of the three objects $x(t)$ versus $y(t)$. Make sure you subtract the center of motion movement $\vec{R}_{cm}=m_1\vec{r}_1+m_2\vec{r}_2+m_3\vec{r_3}$ when plotting the orbits.

Do you see any change of the orbit of Asteroid over 5000 years?

Sketch of the solution¶

For Keppler problem we use $\vec{r}_1$ and $\vec{r}_2$ for vector positions of the two bodies, like the Earth and the Sun. The Newton's Eq require:

\begin{eqnarray} m_1 \ddot{\vec{r}_1} = -G m_1 m_2 \frac{\vec{r}_1-\vec{r}_2}{|\vec{r}_1-\vec{r}_2|^3}\\ m_2 \ddot{\vec{r}_2} = -G m_1 m_2 \frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|^3} \end{eqnarray}

The center of mass can be fixed at the origin $m_1 \vec{r}_1+m_2 \vec{r}_2=0$.

The relative vector $$\vec{r}=\vec{r}_1-\vec{r}_2$$ then satisfies the equation \begin{equation} \ddot{\vec{r}}=G(m_1+m_2)\frac{\vec{r}}{|\vec{r}|^3} \end{equation} The gravitational constant $G$ can be obtained from information about the Earth's orbit. We know that Earth's orbit is almost circular. For circular orbits we know that centrifugal force and gravitational force have to be balanced, which means $m_E \omega^2 R= G m_S m_E/R^2$, where $R$ is Sun-Earth distance $m_S$ and $m_E$ are Sun and Earth mass, and $\omega$ frequency of the Earth's rotation, which is $\omega=2\pi/T$, with $T$ being one year. We thus have $$G=\frac{R^3}{m_S} \left(\frac{2\pi}{T}\right)^2 .$$

We want to meassure time in years $T$, and distance in Sun-Earth distance $R$. We thus define $$\vec{r}=R \vec{r_d}$$ $$t = T t_d $$ and in this AU units the above equation for relative vector is $$\ddot{\vec{r}_d} = -4\pi^2 (1+\frac{m_E}{m_S})\frac{\vec{r}_d}{r_d^3}$$ The ratio of the mass $m_E/m_S= 3\, 10^{-6}$ is safe to neglect.

The initial conditions will be choosen so that we have a circular orbit. For example $\vec{r}_d=[1,0]$ and $\dot{\vec{r}_d}=[0,2\pi]$. The velocity for a circular motion is $v=\omega R$, which is $v=2\pi R/T$, and hence in AU units is just $\dot{y}_0=2\pi$.

We can choose the following set of variables $[x,y,\dot{x},\dot{y}]$, where $\vec{r}_d=[x,y]$. With these we can solve the Kepler equations by: \begin{eqnarray} \begin{bmatrix}\frac{d x}{dt}\\\frac{d y}{dt}\\ \frac{d \dot{x}}{dt}\\ \frac{d \dot{y}}{dt}\\ \end{bmatrix} =\begin{bmatrix} \dot{x}\\\dot{y}\\ -4\pi^2 \frac{x}{(x^2+y^2)^{3/2}}\\ -4\pi^2 \frac{x}{(x^2+y^2)^{3/2}} \end{bmatrix} \end{eqnarray} with initial conditions $y_0=[1,0,0,2\pi]$

In the Keppler's three body problem, we have Sun, Earth and Jupiter with $$\frac{m_2}{m_1}=3\; 10^{-6}$$ and $$\frac{m_3}{m_1}=9.55\; 10^{-4}.$$

The Newton's equations are \begin{eqnarray} && \ddot{\vec{r}}_1 = -G m_1\left( \frac{m_2}{m_1}\frac{\vec{r}_1-\vec{r}_2}{|\vec{r}_1-\vec{r}_2|^3} +\frac{m_3}{m_1}\frac{\vec{r}_1-\vec{r}_3}{|\vec{r}_1-\vec{r}_3|^3}\right) \\ && \ddot{\vec{r}}_2 = -G m_1\left(\frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|^3} +\frac{m_3}{m_1}\frac{\vec{r}_2-\vec{r}_3}{|\vec{r}_2-\vec{r}_3|^3}\right) \\ && \ddot{\vec{r}}_3 = -G m_1\left(\frac{\vec{r}_3-\vec{r}_1}{|\vec{r}_3-\vec{r}_1|^3} +\frac{m_2}{m_1}\frac{\vec{r}_3-\vec{r}_2}{|\vec{r}_3-\vec{r}_2|^3}\right) \\ \end{eqnarray} and in AU units become \begin{eqnarray} && \ddot{\vec{r}}_1 = -4\pi^2\left( \frac{m_2}{m_1}\frac{\vec{r}_1-\vec{r}_2}{|\vec{r}_1-\vec{r}_2|^3} +\frac{m_3}{m_1}\frac{\vec{r}_1-\vec{r}_3}{|\vec{r}_1-\vec{r}_3|^3}\right) \\ && \ddot{\vec{r}}_2 = -4\pi^2\left(\frac{\vec{r}_2-\vec{r}_1}{|\vec{r}_2-\vec{r}_1|^3} +\frac{m_3}{m_1}\frac{\vec{r}_2-\vec{r}_3}{|\vec{r}_2-\vec{r}_3|^3}\right) \\ && \ddot{\vec{r}}_3 = -4\pi^2\left(\frac{\vec{r}_3-\vec{r}_1}{|\vec{r}_3-\vec{r}_1|^3} +\frac{m_2}{m_1}\frac{\vec{r}_3-\vec{r}_2}{|\vec{r}_3-\vec{r}_2|^3}\right) \\ \end{eqnarray}

The initial condistions for Sun can be $\vec{r}_1=0$ and $\dot{\vec{r}}_1=0$, for Earth as before $\vec{r}_2=[1,0]$ and $\dot{\vec{r}}_2=[0,2\pi]$, and for Jupiter $\vec{r}_3=[R_J/R,0]$ and $\dot{\vec{r}}_3=[0,\frac{2\pi}{\sqrt{R_J/R}}]$, where $R_J/R\approx 5.2$. The velocity folows from circular motion of planets, for which we know that gravitational force and centrigunal force are balanced $\omega^2 R_i = G m_S/R_i^2$ or $v_i^2/R_i=G m_S/R_i^2$ and hence $v_i^2 R_i=G m_S = R^3 (2\pi/T)^2$. In AU units each planet with circular orbit should satisfy $v_i=\frac{2\pi}{\sqrt{R_i}}$

Because asteroids are very light, we can neglect their mass compared to mass of Sun and Jupiter. In this case we can treat each Asteroid as an independent problem. We need to simulate the three bodies: Sun, Jupiter, and Asteroid. The masses are $m_2/m_1=9.55\; 10^{-4}$ and $m_3/m_1=0$.

The equations are identical to previous equations simulating Sun, Eart and Jupiter.

The initial condistions for circular orbits are also $R_{asteroid\;(2/1)}=R_{Saturn}/2^{2/3}$, $R_{asteroid\;(3/1)}=R_{Saturn}/3^{2/3}$ and $v_{asteroid}=2\pi/\sqrt{R_{asteroid}}$.

In [26]:

from math import *
from numba import njit
@njit
def Keppler(t, y):
    """anharmonic EOM"""
    (x1,y1,vx1,vy1)=y
    c=4*np.pi**2
    t = np.array([x1,y1])/(x1**2+y1**2)**(3/2.)
    a = -c*t
    return np.array([vx1,vy1,a[0],a[1]])

@njit
def ThreeBody(t, y, m2_m1=9.55e-4, m3_m1=0):
    """anharmonic EOM"""
    (x1,y1,vx1,vy1,x2,y2,vx2,vy2,x3,y3,vx3,vy3)=y
    c=4*np.pi**2
    r12=np.array([x1-x2,y1-y2])
    r13=np.array([x1-x3,y1-y3])
    r23=np.array([x2-x3,y2-y3])
    t12=r12/(r12[0]**2+r12[1]**2)**(3/2.)
    t13=r13/(r13[0]**2+r13[1]**2)**(3/2.)
    t23=r23/(r23[0]**2+r23[1]**2)**(3/2.)
    a1 = -c*(m2_m1*t12+m3_m1*t13)
    a2 = -c*(-t12+m3_m1*t23)
    a3 = -c*(-t13-m2_m1*t23)
    #print(a1,a2,a3)
    # center of mass
    xCM = x1+m2_m1*x2
    return np.array([vx1,vy1,a1[0],a1[1],vx2,vy2,a2[0],a2[1],vx3,vy3,a3[0],a3[1]])
    

In [27]:

# Keppler problem for 5000 years
tf=5000.
sol=solve_ivp(Keppler, [0,tf], [1,0,0,2*pi], rtol=1e-14, atol=5e-7)
print(sol)

  message: 'The solver successfully reached the end of the integration interval.'
     nfev: 956570
     njev: 0
      nlu: 0
      sol: None
   status: 0
  success: True
        t: array([0.00000000e+00, 8.31780817e-03, 3.87543305e-02, ...,
       4.99996119e+03, 4.99999252e+03, 5.00000000e+03])
 t_events: None
        y: array([[ 1.        ,  0.99863464,  0.97049986, ...,  0.23956125,
         0.42491401,  0.46699369],
       [ 0.        ,  0.05223854,  0.24110146, ..., -0.96955822,
        -0.90381473, -0.882808  ],
       [ 0.        , -0.32822444, -1.51488601, ...,  6.10367325,
         5.68979695,  5.55755286],
       [ 6.28318531,  6.27460647,  6.09783093, ...,  1.50810725,
         2.67496071,  2.93986513]])
 y_events: None

In [28]:

# last 5 of 5000 years
plt.plot(sol.t, sol.y[0])
plt.xlim([4995,5000]);
plt.xlabel('t(years)')
plt.ylabel('x(Earth-Sun distance)');

In [29]:

# Eart's orbit for 5000 years
fig, ax = plt.subplots(1,1)
ax.set_title('Earts orbit over 5000 years')
ax.plot(sol.y[0],sol.y[1],lw=0.1)
ax.set_aspect('equal', adjustable='box')

In [30]:

# Sun, Earth, Jupyter
# Initial condistions
rS=0
vS=0
rE=[1,0]
vE=[0,2*pi]
rJ=[5.2,0]
vJ=[0,2*pi/sqrt(5.2)]
# combined into initial condition
y0=[rS,rS,vS,vS,rE[0],rE[1],vE[0],vE[1],rJ[0],rJ[1],vJ[0],vJ[1]]
print('y0=', y0)
m2m1=3.0e-6
m3m1=9.55e-4

tf=5000. # total time in Earth years

sol1=solve_ivp(ThreeBody, [0,tf], y0, rtol=1e-14, atol=5e-7, args=(m2m1,m3m1))

print(sol1)
Rcm1=(sol1.y[0:2]+sol1.y[4:6]*m2m1+sol1.y[8:10]*m3m1)/(1+m2m1+m3m1)

y0= [0, 0, 0, 0, 1, 0, 0, 6.283185307179586, 5.2, 0, 0, 2.7553590302269777]

/Users/haule/opt/anaconda3/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:47: UserWarning: At least one element of `rtol` is too small. Setting `rtol = np.maximum(rtol, 2.220446049250313e-14)`.
  warn("At least one element of `rtol` is too small. "

  message: 'The solver successfully reached the end of the integration interval.'
     nfev: 851918
     njev: 0
      nlu: 0
      sol: None
   status: 0
  success: True
        t: array([0.00000000e+00, 9.28370501e-03, 4.32735233e-02, ...,
       4.99996166e+03, 4.99999673e+03, 5.00000000e+03])
 t_events: None
        y: array([[ 0.00000000e+00,  6.51876508e-08,  1.41563469e-06, ...,
         9.86445821e-03,  9.87720317e-03,  9.87829715e-03],
       [ 0.00000000e+00,  1.97744447e-10,  1.99910273e-08, ...,
         1.32375709e+01,  1.32377559e+01,  1.32377732e+01],
       [ 0.00000000e+00,  1.40431181e-05,  6.53932239e-05, ...,
         3.89817009e-04,  3.37026763e-04,  3.32103772e-04],
       ...,
       [ 0.00000000e+00,  2.55798373e-02,  1.19223645e-01, ...,
         1.40050170e+01,  1.39093223e+01,  1.39003871e+01],
       [ 0.00000000e+00, -1.35542231e-02, -6.31741362e-02, ...,
        -4.07442073e-01, -3.56540045e-01, -3.51787252e-01],
       [ 2.75535903e+00,  2.75532569e+00,  2.75463473e+00, ...,
        -2.72512576e+00, -2.73225166e+00, -2.73286767e+00]])
 y_events: None

In [31]:

fig, ax = plt.subplots(1,1)
ax.plot(sol1.y[0]-Rcm1[0], sol1.y[1]-Rcm1[1],lw=0.1)
ax.plot(sol1.y[4]-Rcm1[0], sol1.y[5]-Rcm1[1],lw=0.1)
ax.plot(sol1.y[8]-Rcm1[0], sol1.y[9]-Rcm1[1],lw=0.1)
ax.set_title('orbits of Sun, Earth and Jupiter over 5000 years')
ax.set_aspect('equal', adjustable='box')

In [32]:

fig, ax = plt.subplots(2,1)

ax[0].plot(sol1.y[4]-Rcm1[0], sol1.y[5]-Rcm1[1],lw=0.1)
ax[0].plot(sol.y[0], sol.y[1], lw=0.1)
ax[0].set_aspect('equal', adjustable='box')
ax[0].set_title('Comparison of Earth\'s orbit with and withouth Jupiter')
ax[0].axes.get_xaxis().set_visible(False)
ax[1].plot(sol1.t, sol1.y[4]-Rcm1[0], label='w Jupiter')
ax[1].plot(sol.t, sol.y[0], label='no Jupiter')
ax[1].set_xlim([4990,5000])
ax[1].set_title('Comparison of x(t) with and withouth Jupiter')
ax[1].set_ylabel('x(Earth-Sun distance)')
ax[1].set_xlabel('time(years)')
ax[1].legend(loc='best');

In [33]:

# distance from Sun to asteroid
ra=5.2/2**(2/3.)
# Initial condistions
rS=0
vS=0
rJ=[5.2,0]
vJ=[0,2*pi/sqrt(5.2)]
rA=[ra,0]
vA=[0,2*pi/sqrt(ra)]
# combined into initial condition
y0=[rS,rS,vS,vS,rJ[0],rJ[1],vJ[0],vJ[1],rA[0],rA[1],vA[0],vA[1]]
m2m1=9.55e-4
m3m1=0

tf=5000. # total time in Earth years

sol2=solve_ivp(ThreeBody, [0,tf], y0, rtol=1e-14, atol=5e-7, args=(m2m1,m3m1))

print(sol2)
Rcm2=(sol2.y[0:2]+sol2.y[4:6]*m2m1+sol2.y[8:10]*m3m1)/(1+m2m1+m3m1)

/Users/haule/opt/anaconda3/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:47: UserWarning: At least one element of `rtol` is too small. Setting `rtol = np.maximum(rtol, 2.220446049250313e-14)`.
  warn("At least one element of `rtol` is too small. "

  message: 'The solver successfully reached the end of the integration interval.'
     nfev: 157676
     njev: 0
      nlu: 0
      sol: None
   status: 0
  success: True
        t: array([0.00000000e+00, 1.96342578e-02, 2.11350743e-01, ...,
       4.99976071e+03, 4.99995687e+03, 5.00000000e+03])
 t_events: None
        y: array([[ 0.00000000e+00,  2.68751868e-07,  3.11086368e-05, ...,
         9.60937759e-03,  9.76187890e-03,  9.78846245e-03],
       [ 0.00000000e+00,  9.32011252e-10,  1.16177293e-06, ...,
         1.31419491e+01,  1.31429581e+01,  1.31431819e+01],
       [ 0.00000000e+00,  2.73755649e-05,  2.94072052e-04, ...,
         9.07792146e-04,  6.45653987e-04,  5.86971088e-04],
       ...,
       [ 0.00000000e+00,  6.81560868e-02,  7.27586711e-01, ...,
         9.91709418e+00,  9.94662118e+00,  9.97276141e+00],
       [ 0.00000000e+00, -7.20287362e-02, -7.68929316e-01, ...,
         3.44853713e+00,  3.49254338e+00,  3.48093175e+00],
       [ 3.47153484e+00,  3.47078264e+00,  3.38474843e+00, ...,
        -2.14889946e-01,  5.22784274e-01,  6.89514721e-01]])
 y_events: None

In [34]:

fig, ax = plt.subplots(1,1)
ax.plot(sol2.y[0]-Rcm2[0], sol2.y[1]-Rcm2[1],lw=0.1)
ax.plot(sol2.y[4]-Rcm2[0], sol2.y[5]-Rcm2[1],lw=0.1)
ax.plot(sol2.y[8]-Rcm2[0], sol2.y[9]-Rcm2[1],lw=0.1)
ax.set_title('Sun, Asteroid in 2/1 gap, and Jupiter orbit over 5000 years')
ax.set_aspect('equal', adjustable='box')

In [35]:

# distance from Sun to asteroid
ra=5.2/3**(2/3.)
# Initial condistions
rS=0
vS=0
rJ=[5.2,0]
vJ=[0,2*pi/sqrt(5.2)]
rA=[ra,0]
vA=[0,2*pi/sqrt(ra)]
# combined into initial condition
y0=[rS,rS,vS,vS,rJ[0],rJ[1],vJ[0],vJ[1],rA[0],rA[1],vA[0],vA[1]]
# masses
m2m1=9.55e-4
m3m1=0

tf=5000.  # total time in Earth years

sol3=solve_ivp(ThreeBody, [0,tf], y0, rtol=1e-14, atol=5e-7, args=(m2m1,m3m1))
print(sol3)
Rcm3=(sol3.y[0:2]+sol3.y[4:6]*m2m1+sol3.y[8:10]*m3m1)/(1+m2m1+m3m1)

/Users/haule/opt/anaconda3/lib/python3.8/site-packages/scipy/integrate/_ivp/common.py:47: UserWarning: At least one element of `rtol` is too small. Setting `rtol = np.maximum(rtol, 2.220446049250313e-14)`.
  warn("At least one element of `rtol` is too small. "

  message: 'The solver successfully reached the end of the integration interval.'
     nfev: 223622
     njev: 0
      nlu: 0
      sol: None
   status: 0
  success: True
        t: array([0.00000000e+00, 1.70626240e-02, 1.49703289e-01, ...,
       4.99984203e+03, 4.99997571e+03, 5.00000000e+03])
 t_events: None
        y: array([[ 0.00000000e+00,  2.02962178e-07,  1.56157131e-05, ...,
         9.70380591e-03,  9.79297876e-03,  9.80657205e-03],
       [ 0.00000000e+00,  6.11667949e-10,  4.12989005e-07, ...,
         1.31424463e+01,  1.31431382e+01,  1.31432646e+01],
       [ 0.00000000e+00,  2.37901038e-05,  2.08512973e-04, ...,
         7.57225039e-04,  5.76278679e-04,  5.43064607e-04],
       ...,
       [ 0.00000000e+00,  6.77971086e-02,  5.89307810e-01, ...,
         1.15109506e+01,  1.11483526e+01,  1.10917889e+01],
       [ 0.00000000e+00, -1.07684101e-01, -9.36012023e-01, ...,
         2.60430499e+00,  3.18245218e+00,  3.27253268e+00],
       [ 3.97391538e+00,  3.97245338e+00,  3.86189991e+00, ...,
        -3.01281930e+00, -2.39131634e+00, -2.26592707e+00]])
 y_events: None

In [36]:

fig, ax = plt.subplots(1,1)
ax.plot(sol3.y[0]-Rcm3[0], sol3.y[1]-Rcm3[1],lw=0.1)
ax.plot(sol3.y[4]-Rcm3[0], sol3.y[5]-Rcm3[1],lw=0.1)
ax.plot(sol3.y[8]-Rcm3[0], sol3.y[9]-Rcm3[1],lw=0.1)
ax.set_title('Sun, Asteroid in 3/1 gap, and Jupiter orbit over 5000 years')
ax.set_aspect('equal', adjustable='box')

In [ ]: