( previous, home, next )

calculus spiral diagram from the archimedes palimpsest

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 leads to the ideas we know call "calculus". Calculus was a climactic product 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 De Moivre, Euler, the Bernoulli's, Lagrange, and Laplace.

One of the most important creations of calculus is the ordinary differential equation. An ordinary differential equation relates condition of some part of the world to the rate of change of that condition, moment by moment. We let \(x(t)\) represent the condition, or state, at time \(t\). The state \(x\) is sometimes a single number and sometimes a vector with 2 or more entries. Sometimes, the independent variable of an ordinary differential equation represents distance or age instead of time. We then try to find some function \(V(x,t)\) that tells us the rate of change of that state, potentially depending on the state and time. You can think of this function \(V(x,t)\) as the wind field that is carrying our dandelion seed. Newton suggested we write the rate of change of the state as \(\dot{x}\), while Liebnitz suggested we write the rate of change as \(dx/dt\), where \(dx\) is the infinitessimal amount that state will change in the next instant, while \(dt\) is the infinitessimal amount time will change in that instant. So, our ordinary differential equation may be written as \[\dot{x} = V(x,t) \quad\text{or}\quad \frac{dx}{dt} = V(x,t).\] Liebnitz's notation is particularly intuitive, because it suggests rewriting the differential equation in infinitessimal form \[dx = V(x,t) dt.\] 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 time \(T\) that we care about. \[\begin{gather} \int_{0}^{T} dx = \int_{0}^{T} V(x,t) dt \\ x(T) - x(0) = \int _ 0^T V(x(t), t) \mathrm{d} t. \end{gather}\] So, we have four roughly interchangable ways of writting an ordinary differential equation. With further study, you'll learn that the last two of these four ways are actually the most useful, but it is standard to use the first two.

Sometimes, what we can say about a state and its rate of change is captured by a less specific form \[F\left(t,x(t), \frac{dx}{dt}\right) = 0,\] but these are still differential equations that can be analyzed to tell us something about how our system is changing.

Closed-form general solutions

A "solution" of an ordinary differential equation is a curve, expressed implicitly or explicitly, for which the given relationship between the state and its rate of change can be shown to hold. A general solution of a differential equation represents all possible solutions for all possible initial conditions and parameter values of interest. A "closed form" general solution is one that can be written exactly in a finite number of terms, possibly including common functions. A closed-form general solution allows us to quickly calculate all sorts of useful relationships that we want to know. But very few differential equations have closed-form solutions. At first, this seems a liability -- it prevents us from using our ordinary differential equation to predict the future. But on the other hand, it means differential equations have much greater expressive power than our standard library of functions. This is a great advantage for modeling.

If you have studied differential equations before, this may sound shocking, since the whole course usually revolves around techniques for finding closed-form solutions. But the equations you primarily study in that course are autonomous linear systems. Autonomous linear systems are special because they have translational symmetry with respect to time and linear scaling symmetry with respect to state. Without both of these symmetries, closed-form solutions are much harder to find. For example, the Airy equation \(d^2x/dt^2 = t x\) is a simple linear differential equation but is not autonomous and does not have a closed-form general solution.

Even though ordinary differential equations often have no closed-form solutions, we have come up with a variety of alternative "solution" methods that are useful. One classic technique was to create new "special" functions that let us express solutions of a particular equation in closed form. For example, the Airy equation's general solution is \(x(t) = C_A Ai(x) + C_B Bi(x)\), where \(Ai(x)\) and \(Bi(x)\) are linearly independent Airy functions. Other important special functions include Bessel functions, and Lambert's W function.

Numerically solution

Another approach is to numerically approximate solutions step-by-step. One of the most straightforward methods to numerically approximate first-order ordinary differential equations is the finite difference method. Since
\[\frac{dx}{dt}=\lim _ {\Delta t\rightarrow 0}\frac{x(t+\Delta t)-x(t)}{\Delta t},\] then for small time steps \(\Delta t\), \[\frac{dx}{dt}\approx \frac{x(t+\Delta t)-x(t)}{\Delta t},\] with the approximation improving as \(\Delta t\) gets small.

To approximate the solution of an ordinary differential equation \(\dfrac{dx}{dt}=V(x(t),t)\) from an initial condition \(x(t_0)=x_0\), 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 the initial condition, 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*}\]

Example 1: 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, \quad t_{n+1} = t_n + \Delta t.\]

The table below shows the results \(N_n\) of applying this formula 10 times, compared to the values calculated using the exact solution, \(N(t_n)\). The numerically calculated values are close to the exact values, but not perfect.

\(n\) \(t _ n\) \(N_n\) \(N(t_n)\) Error

0

0.0

5.000

5.000

0.000

1

0.5

5.625

5.656

0.031

2

1.0

6.311

6.377

0.066

3

1.5

7.058

7.163

0.105

4

2.0

7.868

8.012

0.144

5

2.5

8.738

8.923

0.185

6

3.0

9.667

9.891

0.224

7

3.5

10.650

10.910

0.260

8

4.0

11.680

11.971

0.291

9

4.5

12.750

13.065

0.315

10

5.0

13.850

14.180

0.330

Note that the error gets larger as the time from the initial condition increases. This is a common (but not universal) feature of numerical approximations to the solutions of differential equation.

There are actually many different ways to construct finite-difference methods. The method above is called Euler's method. Euler's method is very simple to explain and implement, but is less accurate than 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. 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, we can use integrate.odeint(), as illustrated in the following code.

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#!/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 numerically approximated solution of the logistic equation when N(0) = 10, K=30, and r = 0.3.

Result of our numerically approximated solution of the logistic equation when N(0) = 10, K=30, and r = 0.3.

Example 2: Van der Pol's system

The same approach can be used to approximate the solutions of 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 this system has no closed-form solution. However, we can numerically approximate the solutions.

[Show code]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#!/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()
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.

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.

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. One thing you may notice is that if your step size is too large, your plots stop looking like smooth curves. And when the step size is large, your approximate solutions may have so much error that they no longer resemble the true solutions. On the other hand, if your step size is too small, the solution algorithm becomes slower and very small step sizes can actually increase the over-all error.

Dimensional analysis of differential equations

The dimensional analysis techniques developed in the earlier lecture can be applied directly to 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_0,t_0,N_1,t_1,r,K) = 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 property of all autonomous ODE systems, and the reason we will often 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.\]

If \(\hat{N}(\hat{t},\hat{N} _ 0)\) is the solution of this non-dimensionalized logistic equation for initial condition \(\hat{N} _ 0\), then when converted back to dimensional variables and our original calendar timescale, \(N_1 = K \hat{N}( r (t_1-t _ 0), N_0/K)\).

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 )