Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 Sigurd Angenent.

The logistic equation and the leap frog method

The logistic differential equation is a diffeq of the form $$ \frac{dx}{dt} = f(x) $$ where $$ f(x) = x(1-x). $$ The leap frog method for approximating solutions to this diffeq begins by replacing the derivative $dx/dt$ by a centered difference: $$ \frac{dx}{dt} \approx \frac{x(t+\Delta t) - x(t-\Delta t)}{2\Delta t}. \tag{1} $$ If we write $x_n = x(n\Delta t)$ for the successive approximations, then the $x_n$ satisfy

\begin{equation} \frac{x_{n+1} - x _ {n-1}}{2\Delta t} = f(x_n) \implies x _ {n+1} = x _ {n-1} + 2f(x_n). \tag{2} \label{eq:leapfrog} \end{equation}

If we know $x_0$ and $x_1$ then this method tells us how to find $ x_2, x_3, \dots$. Here $x_0$ should be a given as initial data. We have to somehow come up with $x_1$. Below we will use the midpoint rule to approximate $x_1$.

The symmetric difference in (1) is a better approximation of the derivative than the forward difference $$ \frac{dx}{dt} \approx \frac{ x(t+\Delta t) - x(t)}{\Delta t} $$ which is used in the derivation of Euler's method. So you would think that the leap frog method is more accurate than the Euler method. But the leap frog method turns out to be problematic. In this notebook we compare the solutions of the leap frog method with those provided by the mid point rule.

The simulation


In [2]:
#  To begin the simulation we import the necessary pyplot stuff
import numpy
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16

Next, we define a function that returns the right hand side of the logistic ODE, and also a function that performs one step of the leapfrog method.


In [3]:
def f(x):
    """the right hand side in the logistic differential equation"""
    return 3*x*(1.0-x)
function_formula="$3x-x^2$"  #  this string goes in the title of our figure

def leapfrog(u0,u1,f,dt):
    """ Compute the solution at the next time step
    by applying the leap frog method to the previous two known solution points. """
    return u0 + 2*dt*f(u1)

We need the midpoint method to estimate $x_1$, and we will also use it to compare the leapfrog solution with a solution we trust.


In [4]:
def midpoint_step(u, f, dt):
    """Returns the solution at the next time-step using the midpoint method.
    
    Parameters
    ----------
    u : array of float
        solution at the previous time-step.
    f : function
        function to compute the right hand-side of the system of equation.
    dt : float
        time-increment.
    
    Returns
    -------
    u_n_plus_1 : array of float
        approximate solution at the next time step.
    """
    umid = u+dt*f(u)/2
    return u + dt * f(umid)

At this point we begin the computation. We will compute the solution over a time interval $0\lt t\lt T$, taking $N$ steps.

We store the results of the leap frog and midpoint methods in two arrays ulp and umdp


In [5]:
# the parameters
N=200
T=20.0
dt = T/N
u0 = 0.02  # initial value of the solution

# create arrays with indices 0, 1, 2, ..., N: 
# they must have N+1 entries since we include 0.
times = numpy.linspace(0.0, T, N+1)
ulp = numpy.zeros(N+1)
umdp = numpy.zeros(N+1)

# assign initial data
ulp[0] = u0
umdp[0] = u0

In [6]:
# Compute umdp[1] and use it as u1 for the leap frog method

umdp[1] = midpoint_step(u0, f, dt)
ulp[1] = umdp[1]

# Compute the remaining solution points
for i in range(2,N+1):
    umdp[i] = midpoint_step(umdp[i-1], f, dt)
    ulp[i] = leapfrog(ulp[i-2], ulp[i-1], f, dt)

In [7]:
# now we plot
# visualization of the path
plt.figure(figsize=(12,6))
plt.grid(True)
plt.xlabel(r't', fontsize=18)
plt.ylabel(r'u', fontsize=18)
plt.title("Approximate solutions of $3x-x^2$.     Time step: $\Delta t=$%.3f "%(dt))
plt.plot(times, ulp, 'r-', lw=1, label='Leap frog')
plt.plot(times, umdp, 'g-', lw=2, label='Mid point rule')
# uncomment the following two lines to plot the
# even entries and the odd entries in the leapfrog approximation
#plt.plot(times[0::2], ulp[0::2], 'y-', lw=3)
#plt.plot(times[1::2], ulp[1::2], 'm-', lw=3)
plt.legend(loc='upper left');


What is going on?

While the solution of the leapfrog method oscillates a lot, it looks like it oscillates between two nice graphs (you can see them by uncommenting two lines in the code above.)

To see what is happening, take the sequence of leap-frog values $$ \dots, x _ {n-2}, x _ {n-1}, x _ {n}, x _ {n+1}, x _ {n+2}, \dots $$ and separate the even from the odd numbered entries: let $$ p _ k = x _ {2k}, \qquad q _ k = x _ {2k+1}. $$ So, $$ x_0, x_1, x_2, x_3, x_4, \dots = p_0, q_0, p_1, q_2, p_3, \dots $$ Write the leap frog equation \eqref{eq:leapfrog} twice: once for $n$ even, and once for $n$ odd. You get \begin{align*} x _ {2k} &= x _ {2k-2} + 2\Delta t f(x _ {2k-1}) & \implies p _ k &= p _ {k-1} + 2\Delta t f(q _ {k-1})\\ x _ {2k+1} &= x _ {2k-1} + 2\Delta t f(x _ {2k}) & \implies q _ k &= q _ {k-1} + 2\Delta t f(p _ k) \end{align*} The equations for $p _ k$ and $q_k$ are almost the Euler method for the system of differential equations \begin{align*} p' &= f(q)\\ q' &= f(p) \end{align*} where $p_k\approx p(2k\Delta t)$ and $q_k\approx q(2k\Delta t)$. The leap frog method isn't solving one differential equation $x'=f(x)$, but a system of two diffeqs $p'=f(q)$, $q'=f(p)$! As long as $p=q$ the system of diffeqs reduces to the original single equation, but apparently there are points where a small difference between $p$ and $q$ gets amplified, and the leap frog method shows us a solution of $p'=f(q), q'=f(p)$ with $p\neq q$.

Suggestions

Play with the parameters: what happens when you double $N$? If you change the initial value?

Try different functions $f$.


The cell below loads the style of the notebook.

In [1]:
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())


Out[1]: