7.1. Background and Some Examples#

The Basic ODE Initial Value Problem#

We consider the problem of solving (approximately) the ordinary differential equation

(7.1)#\[\frac{du}{dt} = f(t, u(t)), a \leq t \leq b\]

with the initial condition

(7.2)#\[u(a) = u_0\]

I will follow the common custom of referring to the independent variable as “time”.

For now, \(u(t)\) is real-valued, but little will change when we later let it be vector-valued (and/or complex-valued).

Notation for the solution of an initial value problem#

Sometimes, we need to be more careful and explicit in describing the function that solves the above initial value problem; then the input parameters \(a\) and \(u_0 = u(a)\) will be included of the function’s formula:

\[u(t) = u(t; a, u_0)\]

(It is standard mathematical convention to separate parameters like \(a\) and \(u_0\) from variables like \(t\) by putting the former after a semicolon.

Examples#

Throughout this chapter, numerical methods will be demonstrated and compared based on the following fairly simple example equations.

Some single first order equations#

Example 7.1 (Integration)

If the derivative depends only on the independent variable \(t\), so that

(7.3)#\[\frac{du}{dt} = f(t), a \leq t \leq b\]

the solution is given by integration:

\[u(t) = u_0 + \int_a^t f(s)\, ds.\]

In particular, with \(u_0 = 0\) the value at \(b\) is

\[u(t) = \int_a^b f(t)\, dt,\]

and this gives us a back-door way to use numerical methods for solving ODEs to evaluate definite integrals.

Example 7.2 (The simplest “genuine” ODE)

The simplest case with \(u\) present in \(f\) is \(f(t, u) = f(u) = u\). But it does not hurt to add a constant, so:

(7.4)#\[ \frac{du}{dt} = k u,\; k \text{ a constant.} \]

The solution is

\[ u(t) = u_0 e^{k(t-a)} \]

We will see that this simple example contains the essence of many ideas that are relevant far more generally.

Example 7.3 (A nonlinear equation, with solutions that blow-up in a finite time)

In the previous examples, \(f(t, u)\) is linear in \(u\) (consider \(t\) as fixed); nonlinearities can lead to more difficult behavior. The equation

(7.5)#\[\frac{du}{dt} = u^2, \, u(a) = u_0\]

This can be solved by separation of variables — or for now you can just verify the solution

\[ u(t) = \frac{1}{T - t}, \, T = a + 1/u_0. \]

Note that if \(u_0 > 0\), the only exists for \(t < T\). (The solution is also valid for \(T > 0\), but that part has no connection to the initial data at \(t=a\).)

Example 7.3 warns us that the IVP might not be well-posed when we set the interval \([a, b]\) in advance: all we can guarantee in general is that a solution exists up to some time \(b\), \(b > a\).

Example 7.4 (A “stiff” equation with disparate time scales)

One common problem in practical situations is differential equations where some phenomena happen on a very fast time scale, but only ever at very small amplitudes, so they have very little relevance to the overall solution. One example is decriptions of some chemical reactions, where some reaction products (like free radicals) are producd in tiny quantities and break down very rapidly, so they change on a very fast time scale but are scarcely relevant to the overall solution.

This disparity of time-sales is called stiffness, from the analogy of a mechanical system in which some components are very stiff and so vibrate at very high frequencies, but typically only at very small amplitudes, or very quicky damped away, so that they can often be safely described by assuming that those stiff parts are completely rigid — do not move at all.

One equation that illustrates this feature is

(7.6)#\[\frac{du}{dt} = -\sin t -k(u - \cos t)\]

where \(k\) is large and positive.

This has the family of solutions is

\[u(t) = \cos t + c e^{-k(t-a)}\]

with \(c = u_0 - \cos(a)\) for the initial value problem \(u(a) = u_0\).

These all get close to \(\cos t\) quickly and then stay nearby, but with a rapid and rapidly decaying “transient” \(c e^{-k t}\).

Many of the most basic and widely use numerical methods (including Euler’s Method thet we meet soon) need to use very small time steps to handle that fast transient, even when it is very small because \(u_0 \approx 1\).

On the other hand there are methods that “supress” these transients, allowing use of larger time steps while still getting an accurate description of the main, slower, phenomena. The simplest of these is the Backward Euler Method that we will see in a later section.

Some systems and higher order equations#

Example 7.5 (Motion of a (Damped) Mass-Spring System in One Dimension)

A simple mathematical model of a damped mass-spring system is

\[\begin{split} \begin{split} M \frac{d^2 u}{d t^2} &= -K u - D \frac{d u}{d t} \\ & \text{with initial conditions} \\ u(a) &= u_0 \\ \left. \frac{dy}{dt} \right|_{t=a} &= v_0 \end{split} \end{split}\]

where \(K\) is the spring constant and \(D\) is the coefficient of friction, or drag.

The first order system form can be left in terms of \(y\) and \(y'\) as

\[\begin{split} \frac{d}{dt} \left[\begin{array}{c} u \\ u' \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -K & -D \end{array}\right] \left[\begin{array}{cc} u \\ u' \end{array}\right] \end{split}\]

Exact solutions#

For testing of numerical methods in this and subsequent sections, here are the exact solutions.

They depend on whether

  • \(D < D_0 := 2 \sqrt{K M}\): underdamped,

  • \(D > D_0\) : overdamped, or

  • \(D = D_0\) : critically damped.

We will mostly explore the first two more “generic” cases.

For the underdamped case, the general solution is

\[ u(t) = e^{-(D/(2M)) (t-a)} [ A \cos(\omega(t-a)) + B \sin (\omega(t-a))], \quad \omega = \frac{\sqrt{4KM - D^2}}{2M} \]

For the above initial conditions, \(A = u_0\) and \(B = (v_0 + u_0 D/(2M)/\omega\).

An important special case of this is the undamped system \(M \frac{d^2 u}{d t^2} = -K u\) for which the solutions become

\[ u(t) = A \cos(\omega(t-a)) + B \sin (\omega(t-a)), \quad \omega =\sqrt{K/M} \]

and it can be verified that the “energy”

\[ E(t) = \frac{M}{2}(y'(t))^2 + \frac{K}{2}(u(t))^2 = \frac{1}{2}(K u_1^2 + {M}u_2^2) \]

is conserved: \(dE/dt = 0\). Conserved quantities can provide a useful check of the acccuracy of numerical method, so we will look at this below.

For the overdamped case, the general solution is

\[ u(t) = A e^{\lambda_+ (t-a)} + B e^{\lambda_- (t-a)}, \quad \lambda_\pm = \frac{-D \pm \Delta}{2M},\; \Delta = \sqrt{D^2 - 4KM} \]

For the above initial conditions, \(A = M(v_0 - \lambda_- u_0)/\Delta\) and \(B = u_0 - A\).

Remark 7.1 (Stiffness)

Fixing \(M\) and scaling \(K = D \to \infty\), \(\Delta = D \sqrt{1 - 4M/D} \approx D - 2M\) so

\[ \lambda_- \approx -\frac{D}{M} + 1 \to -\infty, \quad \lambda_+ \approx -1. \]

Thus the time scales of the two exponential decays become hugely different, with the fast term \(B e^{\lambda_- (t-a)}\) becoming negligible compared to the slower decaying \(A e^{\lambda_+ (t-a)}\).

This is a simple example of stiffness, and influences the choice of a good numerical method for solving such equations.

The variable can be rescaled to the case \(K = M = 1\), so that will be done from now on, but of course you can easily experiment with other parameter values by editing copies of the Jupyter notebooks.

Example 7.6 (The Freely Rotating Pendulum)

Both the above equations are constant coefficient linear, which is convenient for the sake of having exact solution to compare with, but one famous nonlinear example is worth exporing too.

A pendulum with mass \(m\) concentrated at a distance \(L\) from the axis of rotation and that can rotate freely in a vertical plane about that axis and with possible friction proportional to \(D\), can be modeled in terms of its angular position \(\theta\) and angular velocity \(\omega = \theta'\) by

\[ M L \theta'' = -M g\sin\theta - D L \theta', \quad \theta(0) = \theta_0, \theta'(0) = \omega_0 \]

or in system form

\[\begin{split} \frac{d}{dt} \left[\begin{array}{c} \theta \\ \omega \end{array}\right] = \left[\begin{array}{c}\omega \\ -\frac{g}{L}\sin\theta -\frac{D}{M} \omega \end{array}\right] \end{split}\]

These notes will mostly look at the frictionelss case \(D=0\), which has conserved energy

\[ E(\theta, \omega) = \frac{ML}{2} \omega^2 - M g\cos\theta \]

For this, the solution fall into three qualitatively different cases depending on whether the energy is less than, equal to, or greater than the “critical energy” \(Mg\), which is the energy of the unstable stationary solutions \(\theta(t) = \pi (\mod 2\pi)\), \(\omega(t) = 0\): “balancing at the top”:

  • For \(E < Mg\), a solution can never reach the top, so the pendulum rocks back and forth, reach maximum height at \(\theta = \pm \arccos(-E/(Mg))\)

  • For \(E > Mg\), solutions have angular speed \(|\omega| \geq \sqrt{E -Mg} > 0\) so it never drops to zero, and so the direction of rotation can never reverse: solutions rotate in one direction for ever.

  • For \(E = Mg\), one special type of solution is those up-side down stationary ones. Any other solution always has \(|\omega| = \sqrt{E - Mg\cos\theta} > 0\), and so never stops or reverses direction but instead approaches the above stationary point asymptotically both as \(t \to \infty\) and \(t \to \infty\). To visualize concretely, the solution starting at the bottom with \(\theta(0) = 0\), \(\omega(0) = \sqrt{2g/L}\) has \(\theta(t) \to \pm \pi\) and \(\omega(t) \to 0\) as \(t \to \pm\infty\).

Remark 7.2 (Separatrices)

This last kind of special solution is known as a separatrix due to separating the other two qualitatively different sorts of solution. They are also known as heteroclinic orbits, for “asymptotically” starting and ending at different stationary solutions in each time direction — or homoclinic if you consider the angle as a “mod \(2\pi\)” value describing a position, so that \(\theta = \pm\pi\) are the same location and the solutions start and end at the same stationary point.

Example 7.7 (A “Fast-Slow” Equation)

The equation

\[ u'' + (K+1) u' + K u = 0, \quad u(0) = u_0, u'(0) = v_0 \]

which has first order system form

\[\begin{split} \frac{d}{dt} \left[\begin{array}{c} u \\ u' \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -K & -(K+1) \end{array}\right] \left[\begin{array}{c} u \\ u' \end{array}\right] \end{split}\]

This has the general solution

\[ u(t) = A e^{-t} + B e^{-Kt} \]

so for large \(K\), it has two very disparate time scales, with only the slower scale of much significance after an initial transient.

This is a convenient “toy” example for testing two refinements to algorithms:

  • Variable time step sizes, so that they can be short during the initial transient and longer later, when only the \(e^{-t}\) behavior is significant.

  • Implicit methods that can effectively suppress the fast but extremely small \(e^{-kt}\) terms while hanling the larger, slower terms accurately.