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, the seed will follow this wind, moving however it blows at each place and time. We might imagine, then, that by knowing how the wind will be blowing, we can predict where the dandelion seed 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 with which we usually work. This is a great advantage for modelers.

Dimensional analysis

The dimensional analysis techniques developed in the earlier lecture can be applied directly do differential equations. This is often very useful because it reduces the number of parameters that need to be considered.

One elementary example of this can be found in the logistic growth equation describing the growth of bacteria in a test tube. Let \(N(t)\) be the number of cells at time \(t\). The population increases in size according to \[\frac{dN}{dt}=rN(1-N/K)\] where \(r\) is the growth rate with dimensions of 1/time while \(K\) is called the carrying capacity and has dimensions of population size. To solve this differential equation, we also need an initial condition specifying the population size \(N_0\) at some time \(t_0\). So, in the abstract, if the population starts at size \(N_0\) at time \(t_0\) and ends at size \(N _ 1\) at time \(t_1 > t_0\), then the logistic growth model specifies a relationship \[f(N_1,t_1,r,K,N_0,t_0) = 0.\]

The first observation we can make is that the logistic growth equation is autonomous – it does not depend explicitly on the absolute time. This means that starting from initial population \(N_0\) at time \(t_0\) and running to time \(t_1\) will give the same answer as starting with the same population size at time \(0\) and running to time \(t = t_1 - t_0\). This is a common property of all autonomous ODE systems, and the reason we so commonly get away with assuming that the initial condition is applied at time \(t=0\).

The second observation we can make is that we can perform dimensional analysis on the logistic growth equation to reveal scaling symmetries in its solutions. Introducing new variables \(\hat{N}\), \(\hat{N} _ 0\) and \(\hat{t}\) such that \(N = K \hat{N}\), \(t = (1/r) \hat{t}\), and \(N_0 = K \hat{N}_0\), \[\frac{d\hat{N}}{d\hat{t}}=\hat{N}(1-\hat{N}), \quad \hat{N}(0) = \hat{N}_0.\]

When converted back to dimensional variables and our original calendar timescale, \(N_1(t) = K \hat{N}( r (t_1-t _ 0) )\).

This is a good way to show that

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 replace the derivative with it’s approximate form and derive an update formula \[\begin{gather} \frac{x(t+\Delta t)-x(t)}{\Delta t}=V(x(t),t) \\ x(t+\Delta t)=x(t)+V(x(t),t) \Delta t . \end{gather}\] 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)+ V(x(t_0),t_0) \Delta t\\ x(t_0+2\Delta t)&=x(t_0+\Delta t)+ V(x(t_0+\Delta t),t_0+\Delta t) \Delta t\\ x(t_0+3\Delta t)&=x(t_0+2\Delta t)+ V(x(t_0+2\Delta t),t_0+2\Delta t) \Delta t\\ \vdots&\quad\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.

There are actually many different ways to construct finite-difference methods. This method is generally called Euler’s method. Euler’s method is very simple to explain and implement, but is often not as accurate as some other methods. The most commonly used methods are often Runge-Kutta methods, which you can look up and read more about if you are interested.

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+ rN_n(1-N_n/K)\Delta t.\]

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 efficiency 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 our “numerical” solution of the logistic equation when N(0) = 10, K=30, and r = 0.3.
Result of our “numerical” solution of the logistic equation when N(0) = 10, K=30, and r = 0.3.

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 Van der Pol system, \[\begin{align*} \dot{X} &= Y, \\ \dot{Y} &= -X + a(1-X^2)Y. \end{align*}\]

Unlike the logistic equation, the solution of this system can not be written in terms of elementary functions and algebraic operations. However, we can integrate it numerically.

[Show code]
#!/usr/bin/env python
"""
Numerical integration of Vander Pol's equation

    X' = Y
    Y' = -X +
"""

from pylab import *
from scipy import integrate

def VectorField(state, t, a, b, w):
    X, Y = state
    vec = array([ \
        Y, \
        -X + a*(1-X**2)*Y + b*sin(w*t) \
    ])
    return vec

def main():
    starttime = 0.
    endtime = 50.
    numsteps = 10000
    time = linspace(starttime,endtime,numsteps)

    a = 1. # set the value of parameter "a"
    b = 0.
    w = 1.

    S0 = array([0,0.1]) # set our initial condition
    # Now we call the numerical integration function
    S = integrate.odeint(VectorField, S0, time, args=(a,b,w))

    ############################
    figure(1)
    subplot(2,1,1)
    plot(time, S[:,0])
    ylabel('X(t)')

    subplot(2,1,2)
    plot(time, S[:,1])
    ylabel('Y(t)')

    xlabel('Time t')
    savefig('VanderPolTimeseries.png',transparent=True)
    savefig('VanderPolTimeseries.pdf')

    ############################
    figure(2)
    plot(S[:,0], S[:,1], '-')
    xlabel('X(t)')
    ylabel('Y(t)')

    savefig('VanderPolPhaseplane.png',transparent=True)
    savefig('VanderPolPhaseplane.pdf')

main()

We can then use the approximate solutions to then investigate some scientific questions.

Time series plots of numerical solutions of the Van der Pol equation when a = 1 for initial conditon X(0)=0, Y(0)=0.1. Phaseplane plot of numerical solution of the Van der Pol equation when m = 1 for initial conditon X(0)=0, Y(0)=0.1.

Try it for yourself.

Exercises

  1. Use dimensional analysis to express the solutions \(y(t)\) of the dimensional equation \(dy/dt = m - r y\) in terms of solutions \(u(s)\) of the dimensionless equation \(du/ds = 1 - u\).

  2. Use algebra and dimensional analysis to show that the 4 parameter system \[\begin{align*} \dot{p} &= \alpha q, \\ \dot{q} &= -\beta p + \gamma q - \delta p^2 q, \end{align*}\] is equivalent to the Van der Pol equation (see above) under a rescaling of the states and time. What is the formula for \(a\) in terms of \(\alpha, \beta, \gamma, \delta\)?

  3. 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.
  4. 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?

  5. One day it started snowing at a 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.)

  6. (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?

  7. (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.

  8. 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.

  9. 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?

  10. (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 )