# Differential equations

( previous, home, next )

Differential equations are difficult, confusing, weak, and unrealistic. They often mislead students as to the nature of systems. Mathematicians have had difficulty defining a derivative and there is a reason. Derivatives do not exist except in a mathematician’s imagination. Nowhere in nature does nature take a derivative. Nature only integrates, that is, accumulates in stocks. Casting behavior in terms of differential equations leaves many students with an ambiguous or even reversed sense of the direction of causality. I have had MIT students argue that water flows out of the faucet because the level of water in the glass is rising; that seems natural to them if the flow has been defined as the derivative of the water level in the glass. -- J. W. Forrester, 2009 Early weather map by including wind directions by Elias Loomis, 1842. First map, drawn by Edmond Halley, attempting to describe wind patterns across the earths oceans, 1686

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

The classical idea of 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 as functions of its independent variable 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 transcendental functions like $$\tan()$$ or $$\sinh()$$. For example, the general closed-form solution of the differential equation $\begin{gather} \frac{d^2y}{dt^2} = y \end{gather}$ is given by $\begin{gather} y(t) = A \sinh(t) + B \cosh(t). \end{gather}$

A closed-form general solution allows us to quickly calculate all sorts of useful relationships that we want to know.

For many differential equations, we have not found closed-form solutions. At first, this seems a liability -- it prevents us from using our ordinary differential equation to predict the future. On the other hand, it means differential equations have much greater expressive power than our standard library of functions, and expressiveness is an 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 old 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 \operatorname{Ai}(x) + C_B \operatorname{Bi}(x)$$, where $$\operatorname{Ai}(x)$$ and $$\operatorname{Bi}(x)$$ are linearly independent Airy functions. We know a great deal about Airy functions, including the locations of their zeros, series approximations, integral representations, derivatives, and asymptotic behaviors -- many of the same things we know about common trigonometric functions. So, if you work with them enough, Airy's odd beasts become familiar friends. The same can be said of many other important special functions including Bessel functions, Jacobi elliptic 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 # 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.

### 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$$. 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.

# 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 )