Day 21: Differential Equations

Yesterday we had a look at Julia’s support for Calculus. The next logical step is to solve some differential equations. We’ll look at two packages today: Sundials and ODE.

Sundials

The Sundials package is based on a library which implements a number of solvers for differential equations. First off you’ll need to install that library. In Ubuntu this is straightforward using the package manager. Alternatively you can download the source distribution.

sudo apt-get install libsundials-serial-dev

Next install the Julia package and load it.

julia> Pkg.add("Sundials")
julia> using Sundials

To demonstrate we’ll look at a standard “textbook” problem: a damped harmonic oscillator (mass on a spring with friction). This is a second order differential equation with general form

$$ \ddot{x} + a \dot{x} + b x = 0 $$

where \(x\) is the displacement of the oscillator, while \(a\) and \(b\) characterise the damping coefficient and spring stiffness respectively. To solve this numerically we need to convert it into a system of first order equations:

$$ \begin{aligned} \dot{x} &= v \ \dot{v} &= - a v - b x
\end{aligned} $$

We’ll write a function for those relationships and assign specific values to \(a\) and \(b\).

function oscillator(t, y, ydot)
  ydot[1] = y[2]
  ydot[2] = - 3 * y[1] - y[2] / 10
end

oscillator (generic function with 2 methods)

Next the initial conditions and time steps for the solution.

initial = [1.0, 0.0];                        # Initial conditions
t = float([0:0.125:30]);                     # Time steps

And finally use cvode() to integrate the system.

xv = Sundials.cvode(oscillator, initial, t);
xv[1:5,:]

5x2 Array{Float64,2}:
 1.0 0.0
 0.97676 -0.369762
 0.908531 -0.717827
 0.799076 -1.02841
 0.65381 -1.28741

The results for the first few time steps look reasonable: the displacement (left column) is decreasing and the velocity (right column) is becoming progressively more negative. To be sure that the solution has the correct form, have a look at the Gadfly plot below. The displacement (black) and velocity (blue) curves are 90° out of phase, as expected, and both gradually decay with time due to damping. Looks about right to me!

ODE

The ODE package provides a selection of solvers, all of which are implemented with a consistent interface (which differs a bit from Sundials).

using ODE

Again we need to define a function to characterise our differential equations. The form of the function is a little different with the ODE package: rather than passing the derivative vector by reference, it’s simply returned as the result. I’ve consider the same problem as above, but to spice things up I added a sinusoidal driving force.

function oscillator(t, y)
  [y[2]; - a * y[1] - y[2] / 10 + sin(t)]
end

oscillator (generic function with 2 methods)

We’ll solve this with ode23(), which is a second order adaptive solver with third order error control. Because it’s adaptive we don’t need to explicitly specify all of the time steps, just the minimum and maximum.

a = 1; # Resonant
T, xv = ode23(oscillator, initial, [0.; 40]);
xv = hcat(xv...).'; # Vector{Vector{Float}} -> Matrix{Float}

The results are plotted below. Driving the oscillator at the resonant frequency causes the amplitude of oscillation to grow with time as energy is transferred to the oscillating mass.

If we move the oscillator away from resonance the behavior becomes rather interesting.

a = 3; # Far from resonance

Now, because the oscillation and the driving force aren’t synchronised (and there’s a non-rational relationship between their frequencies) the displacement and velocity appear to change irregularly with time.

How about a double pendulum (a pendulum with a second pendulum suspended from its end)? This seemingly simple system exhibits a rich range of dynamics. It’s behaviour is sensitive to initial conditions, one of the characteristics of chaotic systems.

First we set up the first order equations of motion. The details of this system are explained in the video below.

function pendulum(t, y)
  Y = [
  6 * (2 * y[3] - 3 * cos(y[1] - y[2]) * y[4]) / (16 - 9 * cos(y[1] - y[2])^2);
  6 * (8 * y[4] - 3 * cos(y[1] - y[2]) * y[3]) / (16 - 9 * cos(y[1] - y[2])^2)
  ]
  [
  Y[1];
  Y[2];
  - (Y[1] * Y[2] * sin(y[1] - y[2]) + 3 * sin(y[1])) / 2;
  - (sin(y[2]) - Y[1] * Y[2] * sin(y[1] - y[2])) / 2;
  ]
end

pendulum (generic function with 1 method)

Define initial conditions and let it run…

initial = [pi / 4, 0, 0, 0]; # Initial conditions -> deterministic behaviour

T, xv = ode23(pendulum, initial, [0.; 40]);

Below are two plots which show the results. The first is a time series showing the angular displacement of the first (black) and second (blue) mass. Next is a phase space plot which shows a different view of the same variables. It’s clear to see that there is a regular systematic relationship between them.

Next we’ll look at a different set of initial conditions. This time both masses are initially located above the primary vertex of the pendulum. This represents an initial configuration with much more potential energy.

initial = [3/4 * pi, pi, 0, 0]; # Initial conditions -> chaotic behaviour

The same pair of plots now illustrate much more interesting behaviour. Note the larger range of angles, θ₂, achieved by the second bob. With these initial conditions the pendulum is sufficiently energetic for it to “flip over”. Look at the video below to get an idea of what this looks like with a real pendulum.

It’s been a while since I’ve played with any Physics problems. That was fun. The full code for today is available at GitHub. Come back tomorrow when I’ll take a look at Optimisation in Julia.