Differential equations

( previous, home, next )

calculus spiral diagram from the archimedes palimpsest

Introduction to differential equations

Imagine the winds blowing across the country. (Well, you don’t have to imagine it – some artists have already visualized these winds for us.) At each place in the country, the wind is blowing in some direction with some speed. The speed and direction differ from one place to another, and they can change over time as well, becoming stronger or weaker and switching direction. And if we release a dandelion seed or weather balloon, the motion tracks this wind, moving however it blows at each place and time. We might imagine, then, that just by knowing how the wind is blowing, we can predict where the weather balloon will be an hour or a day after we let it go.

This is the kind of thinking, that, together with the recreation of ancient Greek mathematics, leads to the ideas we know call “calculus”. Calculus was a product of the climax of the scientific revolution in the 17th century, the simultaneous synthesis by Godried Liebnitz and Isaac Newton of the related works of Galileo, Descartes, Pascal, Fermat, Wallis, Barrow, Roberval, Torricelli, Huygens, and others. Calculus itself opened the flood-gates for exacting descriptions of the laws of nature by successors like Euler, the Bernoulli’s, Lagrange, and Fourier.

One of the most important realizations of this concept is a differential equation. A differential equation is a statement that proposes to relate the expected change of the world moment by moment to the world’s current state. We let \(x(t)\) represent the world’s state at time \(t\), and then find some function \(V(x,t)\) that tells us the infinitessimal amount \(dx\) that state will change in the next instant of time \(dt\) as some function of the current state and time. You can think of this function \(V(x,t)\) as the wind field that is carrying our dandelion seed. It leads to a very general infinitessimal equation of the form \[dx = V(x,t) dt \quad \text{or} \quad \frac{dx}{dt} = V(x,t).\] Then, knowing we start at \(x(0)\) at time \(t=0\), we can sum all these little changes for each instant of time together to find where the world will be at the future instant \(T\) that we care about. \[x(0) + \int_{0}^{T} dx = x(0) + \int_{0}^{T} V(x(t), t) \mathrm{d} t = x(T).\]

Sometimes, we aren’t fortunate enough to have our equation in this form. Sometimes, all we have is a more general form \[F\left(t,x(t), \frac{dx}{dt}\right) = 0,\] but these general differential equations can still be analyzed to tell us something about how our system is changing.

The solutions of very few differential equations can be expressed in terms of elementary functions like polynomials and exponentials. Most differential equations do not have “solutions” that can be written down by a quick hand. At first, this seems a liability – it prevents us from giving a nice closed-form solution for the future. But on the other hand, it means differential equations have much greater expressive power than the standard kinds of functions we are used to working with – and that is a great advantage for modelers.

Numerically solution

One of the most straightforward methods of numerically solving first-order ordinary differential equations is the finite difference method using forward differences to represent the derivative. \[\frac{dx}{dt}=\lim_{\Delta t\rightarrow 0}\frac{x(t+\Delta t)-x(t)}{\Delta t}\Rightarrow \frac{dx}{dt}\approx \frac{x(t+\Delta t)-x(t)}{\Delta t},\] with the approximation improving as \(\Delta t\) gets small.

To solve an ordinary differential equation in the form \(\dfrac{dx}{dt}=V(x(t),t)\), we can derive an update formula \[\frac{x(t+\Delta t)-x(t)}{\Delta t}=V(x(t),t)\Rightarrow x(t+\Delta t)=x(t)+\Delta t V(x(t),t).\] Starting with an initial condition \(x(t_0)=x_0\), we can use this update formula to build the solution time step by time step: \[\begin{align*} x(t_0+\Delta t)&=x(t_0)+\Delta t V(x(t_0),t_0)\\ x(t_0+2\Delta t)&=x(t_0+\Delta t)+\Delta t V(x(t_0+\Delta t),t_0+\Delta t)\\ x(t_0+3\Delta t)&=x(t_0+2\Delta t)+\Delta t V(x(t_0+2\Delta t),t_0+2\Delta t)\\ \vdots&=\vdots \end{align*}\] Note: Another way of writing the update rule is as a difference equation. Let \(\Delta t=k\), \(t_n=t_0+nk\), \(x_n=x(t_n)\). Then \[x_{n+1}=x_n + kV(x_n,t_n).\] That’s why this approach is called a finite difference method. (finite difference methods are actually a large toolbox!)

Example: Numerical solution of the logistic equation

Consider the logistic equation, frequently used to model population growth: \[\frac{dN}{dt}=rN(1-N/K)\] with \(N(0)=5\). Set the growth rate \(r=0.3\) per day and the carrying capacity \(K=30\). Solve numerically using the forward difference method described above.

The update formula difference equation is \[N_{n+1}=N_n+\Delta t rN_n(1-N_n/K).\]

See spreadsheet for worked-out solution.

Fun fact: as the time step \(\Delta t\) grows large, this difference equation gives period-doubling bifurcations and chaos!

Error in the method outlined above is \(O(\Delta t)\). More sophisticated methods will converge to the analytic solution more rapidly (e.g. Fourth order Runge–Kutta is \(O(\Delta t^4)\)). More sophisticated methods may also use adaptive time stepping and other techniques to improve accuracy and speed of calculation.

We are not in a numerics class so we’ll take advantage of what’s been coded up for us! In Python, scipy has an integrate toolbox full of tools to numerically solve ordinary differential equations. Depending on the nature of the ordinary differential equations, some methods may be better than others (e.g. if your equation is stiff).

For the Logistic equation, I’ll use integrate.odeint(). See code:

[Show code]
#!/usr/bin/env python3

from scipy import array, linspace
from scipy import integrate
from matplotlib.pyplot import *

def vector_field(x, t, r, K):
    # Pearl-Reed logistic growth differential equation:
    #
    #    dN
    #    -- = r N (1 - N/K)
    #    dt
    #
    # r = the per-capita growth rate at low numbers
    #
    # K = the carrying-capacity
    #
    N = x[0] # Rabbits density
    return array([r*N*(1-N/K)])

# Parameters
c_r = .3  # rabbit growth rate, units of per-day
c_K = 30. # carrying capacity of rabbits, units of rabbits

# set up our initial conditions
N0 = 10. # initial number of rabbits
x0 = array([N0])

# choose the time's we'd like to know the approximate solution
t = linspace(0., 30., 30)

# and integrate
x = integrate.odeint(vector_field, x0, t, args=(c_r, c_K))
# x is an array
# Each row of x is for the corresponding time in t
# Each column of x is the dynamic of the corresponding entry in x0

# now, plot the solution curve
figure(1)
plot(t, x[:,0], 'bx-', linewidth=1)
xlabel('Time (months)', fontsize=24)
ylabel('Number of rabbits', fontsize=24)
tight_layout()
savefig('LogisticGrowth.png', transparent=True, bbox_inches = 'tight')
show()
Result of numerical solution of the logistic equation
Result of numerical solution of the logistic equation

Example: Numerical solution of a system of equations

The same approach can be used to solve systems of first order ordinary differential equations. For example consider the “Competing species” model, \[\begin{align*} \frac{dR}{dt}&=r_1R(1-R/K_1)-c_1RS \\ \frac{dS}{dt}&=r_2S(1-S/K_2)-c_2RS \end{align*}\] \(R\) indicate a population of rabbits, \(S\) a population of sheep. Note that in the absence of the other, both populations experience logistic growth with different growth rates and carrying capacities. However they are competing for the same resource: both eat grass! So the rabbit population growth is impinged upon by the presence of sheep, hence the \(-c_1RS\) term (negative b/c competition slows growth; \(c_1\) is the magnitude of that competition). Similarly the sheep population growth is negatively impacted by the presence of rabbits.

Unlike the logistic equation, the solution of the model equations has no closed-form solution. Solving numerically,

[Show code]
#!/usr/bin/env python3

from scipy import array, linspace
from scipy import integrate
from matplotlib.pyplot import *

def vector_field(x, t, r1, K1, c1, r2, K2, c2):
    # Competing species differential equations model
    # from Section 9.4 of Boyce & DiPrima
    #
    #    dR
    #    -- = r1*R*(1-R/K1) - c1*R*S
    #    dt
    #
    #    dS
    #    -- = r2*S*(1-S/K2) - c2*R*S
    #    dt
    #
    R = x[0] # Rabbits density
    S = x[1] # Sheep density
    V = array([r1*R*(1-R/K1) - c1*R*S,  r2*S*(1-S/K2) - c2*R*S])
    return V


def main():
    # set up our initial conditions
    R0 = 10.
    S0 = 6.
    x0 = array([R0, S0])

    # Parameters
    r1 = .3    # rabbit growth rate
    r2 = .2    # sheep growth rate
    c1 = .2    # inhibition of rabbits due to competition
    c2 = .1    # inhibition of sheep due to competition
    K1 = 30.   # carrying capacity of rabbits
    K2 = 20.   # carrying capacity of sheep

    # choose the time's we'd like to know the approximate solution
    t = linspace(0., 40., 100)

    # and solve
    x = integrate.odeint(vector_field, x0, t, args=(r1,K1,c1,r2,K2,c2))

    # now, plot the solution curves
    figure(1)
    plot(t, x[:,0], 'b-', linewidth=2)
    plot(t, x[:,1], 'g-', linewidth=2)
    axis([0,40,0,31])
    xlabel('Time (days)')
    ylabel('Number')
    title('Initial conditions: 10 rabbits, 5 sheep')
    legend(['Rabbits', 'Sheep'], loc=2, framealpha=0.)
    tight_layout()
    savefig('CompetingSpecies2.png',transparent=True)
    show()

main()

We’d use the code to answer whatever specific scientific question we want to understand. For example, given certain conditions, which population will prevail, will they coexist?

Both the competing species and logistic models are nonlinear. Nonlinearities create beautiful and interesting outcomes. For example, in the competing species model, \[\begin{align*} \frac{dR}{dt}&=0.3R(1-R/30)-0.2RS \\ \frac{dS}{dt}&=0.2S(1-S/20)-0.1RS \end{align*}\] If you start with 10 rabbits, and 5 sheep, the competition offered by the sheep is not strong enough and the rabbit population prevails and reaches carrying capacity. But with just one more sheep, the model predicts that the competition will kill off the rabbit population and only sheep will survive!

Try it for yourself and see…

Exercises

  1. The end of a pocket watch chain is very slowly pulled in a straight line, and the pocket watch, at the end of the chain, follows behind.
    1. Find three parameters that the path of pursuit depends on.
    2. What relationship should these parameters satisfy if pulling the chain end also moves the watch.
    3. Use dimensional analysis to show we may assume, without loss of generality, that the chain length is 1 unit.
    4. Find a differential equation for the path of the center of the watch.
    5. Transform your differential equation back to dimensional form.
  2. Suppose that instead of moving in a straight line, the watch chain end is pulled in a circle with radius \(r\). What path does the watch follow?

  3. One day it started snowing at heavy and steady rate. A snowplow started out at noon, going 2 miles the first hour and 1 mile the second hour. If the speed of the snowplow is inversely proportional to the depth of the snow, at what time did it start snowing? (Agnew’s Snowplow problem – Agnew was a famous leader of the mathematics department at Cornell.)

  4. (Klamkin’s Great Snowplow Chase) One winter day, snow accumulates at steady rate. Three identical snowplows start plowing the same road, the first leaving at 12 noon, the second leaving at 1 pm, and the third leaving at 2 pm. At some time later, they all collide. At what time did the snow start to fall?

  5. (Golbart’s ladies) Four ladybugs are found on a unit circle at angles \(\pi/4\), \(3 \pi/4\), \(5 \pi/4\), and \(7\pi/4\). Each ladybug walks directly towards the next bug counterclockwise from itself. All lady bugs walk at the same speed. Find a differential equation for \(y(x)\) where \((x,y)\) are the coordinates of the top-right bug.

  6. A tank shaped like a hemisphere is filled to the top. There is a small hole in the bottom of the tank, out of which water drains at a rate proportional to the square root of the height of the water. Model the change in the height of the water in the tank over time.

  7. An ant walks along a elastic cord attached to the bumper of a truck at one end and a pole at the other. The truck moves at 1 foot per second. The ant walks at 1 centimeter per second. If the bungie cord is infinitely elastic and never breaks, will the ant ever reach the truck?

  8. (Letter by Thomas Godfrey to the Ben Franklin’s Pennsylvania Gazette, 1743) Mr. Franklin, As privateering is now so much in Fashion, the printing of the following question may be an amusement, if not to the privateers, yet to some of your correspondents or readers. Suppose a privateer, in the latitude of 10 degrees North, should at 6 in the morning spy a ship due south of her, distant 20 miles; upon which she steers directly for her, and runs at the rate of 8 miles an hour. The ship at the same time sees the privateer, but not being much afraid of her, keeps on her course due west, and sails at the rate of 6 miles an hour; how many hours will it be before the privateer overtakes the ship?

( previous, home, next )