7. Ordinary Differential Equations#
7.1. Introduction#
Many of the laws of physics are formulated in terms of differential equations.
Part of the power of computational tools is that it is easy to solve almost every differential equation. E.g. one can go beyond the small displacement in the description of oscillations, and explore interesting nonlinear phenomena.
Let’s begin with a recap of some mathematical preliminaries that pertain to differential equations, and then move on, in this chapter, to to discuss algorithms for solving the ordinary differential equations.
7.2. Mathematical Preliminaries#
ORDER: The order of a differential equation refers to the degree of the derivative.
e.g. the most general form for a first-order differential equation is:
\( \frac{ \mathrm{d} y } { \mathrm{d} x } = f(x,y)\).
e.g. even if \(f(x,y)\) is a nasty function of \(x\) and \(y\):
\( \frac{ \mathrm{d} y } { \mathrm{d} x } = x^2 y + x^{17} + y^{13}\),
the equation is still first order.
A general form of a second-orrder differential equation is:
\( \frac{ \mathrm{d}^2 y } { \mathrm{d} x^2 } + \lambda \frac{ \mathrm{d} y } { \mathrm{d} x } = f(x, \frac{ \mathrm{d} y } { \mathrm{d} x }, y)\).
The derivative on the RHS may involved any power of the first derivative as well.
In the above, \(x\) is known as the independent variable and \(y\) is the dependent variable.
ORDINARY VS. PARTIAL
Ordinary diffrential equations contain just one independent variable, e.g. \(x\) in our discussion.
In contrast, partial differential equations contain several independent variables.
E.g. the Schrödinger equation in three dimensions:
\( i \hbar \frac{ \partial \Psi(\vec{x},t) } {\partial t} = -\frac{ \hbar^2 } { 2 m } \left[ \frac{ \partial^2 \Psi} { \partial x^2 } + \frac{ \partial^2 \Psi} { \partial y^2 } + \frac{ \partial^2 \Psi} { \partial z^2 } \right] + V(\vec{x}) \Psi(\vec{x},t)\),
is a partial differential equation with four independent variables \(\vec{x} = (x,y,z)\) and \(t\), and one dependent variable \(\Psi(\vec{x},t)\).
LINEAR VS. NONLINEAR
Part of the strength of computational methods is that we ar eno longer limited to solving linear equations.
A linear equation is one in which only the first power of \(y\) or of \(\mathrm{d}^n y / \mathrm{d}x^n\) appears, e.g.:
\( \frac{ \mathrm{d} y } { \mathrm{d} x } = g(x) y(x)\) is linear,
\( \frac{ \mathrm{d} y } { \mathrm{d} x } = \lambda y(x) - \lambda^2 y^2(x)\) is nonlinear.
WARNING: An important property of linear equations is no longer valid for nonlinear equations: the principle of superposition, that lets us add solutions together to form new ones: e.g. if \(A(x)\) is a solution and \(B(x)\) is a solution to a linear differential equation, then their sum:
\(y(x) = a A(x) + b B(x)\), is also a solution, for arbritrary values of the constants \(a\) and \(b\). This would no longer be true for a nonlinear equation!
Initial and boundary conditions
The general solution of a first-order differential equaiton always contains one arbitrary constant. The general solution of a second-order differential equation contains two, and so on.
For any specific problem, these constants are fixed by the initial conditions. E.g. for a first-order equation for \(y(x)\), the sole initial condition necessary is \(y(x=x_0)\), i.e. the value of \(y\) at some value of \(x\).
For a second-order equation, two initial conditions are required. These may e.g. be the position and velocity at some time when solving a problem involving Newton’s second law of motion:
\( \frac{ \mathrm{d}^2 y } { \mathrm{d} t^2 } = \frac{1}{m} F(y,t)\).
In addition to the initial conditions, it is possible to further restrict the solutions of differential equations, e.g. through boundary conditions that constrain one solution to have fixed values at the boundaries of the solution space.
7.3. Dynamic form for ODEs#
The most general form of an ODE is a set of \(M\) coupled first-order differential equations:
\( \frac{ \mathrm{d} \vec{y} } { \mathrm{d} x } = \vec{f}(x,\vec{y})\),
where \(x\) is the independent variable and \(\vec{y}\) is a vector of \(M\) dependent variables. The vector of functions \(f\) is an \(M\)-component vector as well.
Differential equations of higher order can be written in this first-order coupled form, by introducing auxiliary functions.
Consider, for example, the one-dimensional motion (in the \(z\) direction), of a particle of mass \(m\) under a force field F(z), described by the second-order differential equation (Newton’s second law):
\( m\frac{ \mathrm{d}^2 z } { \mathrm{d} t^2 } = F(z)\).
If we define the momentum as: \(p(t) = m \frac{ \mathrm{d} z } { \mathrm{d} t }\),
then Newton’s law becomes two coupled first-order equations (known as Hamilton’s equations):
\(\frac{ \mathrm{d} p } { \mathrm{d} t } = F(z)\),
\(\frac{ \mathrm{d} z } { \mathrm{d} t } = \frac{p}{m}\),
which are of the desired general form.
It is therefore sufficient to consider in detail only methods for first-order ODEs!
Since the matrix structure of coupled differential equations is of this natural form, we can focus on the case where there is only one independent variable, which can be generalized readily.
Thus, we will begin by consider the methods for solving:
\(\frac{ \mathrm{d} y } { \mathrm{d} x } = f(x,y)\), for a single dependent variable \(y(x)\).
In this chapter, we will discuss methods for solving ODEs, with emphasis on the initial value problem, i.e.:
Find \(y(x)\) given the value of \(y\) at some initial point, say \(y(x=0) = y_0\).
This is the kind of problem that occurs, e.g. when we are given the initial position and momentum of a particle and we wish to find its subsequent motion.
Later on, we will discuss the equally important boundary value and eigenvalue problems.
7.4. ODE Algorithms#
7.4.1. Euler’s Method#
We are interested in the solution of \(\frac{ \mathrm{d} y } { \mathrm{d} x } = f(x,y)\), with the initial condition \(y(x=0) = y_0\).
More specifically, we are usually interested in the value of \(y\) at particular vlaue of \(x\), say \(x=1\).
The general strategy is to divide the interval \([0,1]\) into a large number, \(N\), of equally spaced subintervals of length \(h=1/N\) and then to develop a recursion formula relating \(y_n\) to \(\{y_{n-1}, y_{n-2},...\}\), where \(y_n\) is our approximation to \(y(x_n = nh)\). Such a recursion formula would then allow a step-by-step integration of the differential equation from \(x=0\) to \(x=1\).
One of the simplest algorithms is known as Euler’s method, in which we consider the differential equation at the point \(x_n\), and replace the derivative on the LHS by its forward-difference approximation:
\(\frac{ y_{n+1} - y_n } { h } + \mathcal{O}(h) = f(x_n , y_n)\),
so that, if we solve for \(y_{n+1}\), we obtain a recursion relation expressing \(y_{n+1}\) in terms of the “previous” value, \(y_n\):
\(y_{n+1} = y_n + h f(x_n, y_n) + \mathcal{O}(h^2)\).
The formula has a “local” error (i.e. the error made in taking a single step from \(y_n \rightarrow y_{n+1}\), that is \(\mathcal{O}(h^2)\). The “global” error made in finding \(y(1)\), after taking \(N\) such steps in integrating from \(x=0\) to \(x=1\) is then \(N \mathcal{O}(h^2) \sim \frac{1}{h} \mathcal{O}(h^2) \sim \mathcal{O}(h)\).
Evidently, this error decreases only linearly with decreasing step size, so that if we halve \(h\) (and perform twice as many steps), we halve the error of the final answer.
Let’s apply this method to a simple problem.
7.4.1.1. Example 7.1: An Application of Euler’s Method#
Apply Euler’s method to solve:
\(\frac{ \mathrm{d} y } { \mathrm{d} x } = -xy\) with initial condition \(y(0) = 1\) from \(x=0\) to \(x=3\).
Use various step sizes: \(h=0.500, 0.200, 0.100, 0.050, 0.020, 0.010, 0.005, 0.002, 0.001\) and calculate the error relative to the analytical solution \(y(x) = \exp{-x^2/2}\) at the points \(x=1\) and \(x=3\). Does it scale as you expect with \(h\)?
Plot the \(y(x)\) for a few step sizes, e.g. \(h=0.500, 0.050, 0.001\), as well as the analytical function.
7.4.1.2. Example 7.2: Evaluation of the Integration Algorithm#
A simple and often stringent test of an accurate numerical integration is to use the final value of \(y\) obtained as the initial condition to integrate backward from the final value of \(x\) to the starting point. The extent to which the resulting value of \(y\) differs from the original initial condition is then a measure of the inaccuracy.
Apply this test to Example 7.1.
Although Euler’s method seems to work quite well, it is generally unsatisfactory due to its lower-order accuracy.
The algorithm is not recommended for general use, but it is commonly used to start off more precise algorithms.
7.4.2. Runge-Kutta Methods#
There is a lot of freedom in writing down algorithms for integrating differential equations. There exists a large number of them, each having their own pecularities and advantages.
One very convenient and widely used class of methods are the Runge-Kutta algorithms, which come in varying orders of accuracy.
Here, we derive the second-order version (rk2) to give the spirit of the approach and then state the equations for the fourth-order method.
We begin by writing down the formal integral of our differential equation:
\(\frac{ \mathrm{d} y } { \mathrm{d} x } = f(x,y)\),
as:
\(y(x) = \int f(x,y) \mathrm{d} x\),
which leads to:
\(y_{n+1} = y_n + \int_{x_n}^{x_{n+1}} f(x,y) \mathrm{d} x\).
We then expand \(f(x,y)\) in a Taylor series about the midpoint of the integration interval:
\(f(x,y) \simeq f(x_{n+1/2}, y_{n+1/2}) + (x-x_{n+1/2}) \left.\frac{ \mathrm{d} f } { \mathrm{d} x } \right|_{x_{n+1/2}} + \mathcal{O}(h^2)\).
Substituting into the expression for \(y_{n+1}\), the linear term integrates to zero because it is equally positive and negative in the interval \([x_n, x_{n+1}]\):
\(\int_{x_n}^{x_{n+1}} (x-x_{n+1/2}) \left.\frac{ \mathrm{d} f } { \mathrm{d} x } \right|_{x_{n+1/2}} \mathrm{d}x = \left.\frac{ \mathrm{d} f } { \mathrm{d} x } \right|_{x_{n+1/2}} \int_{x_n}^{x_{n+1}} (x-x_{n+1/2}) \mathrm{d}x = 0\).
Then:
\(y_{n+1} = y_n + h f(x_{n+1/2}, y_{n+1/2}) + \mathcal{O}(h^3)\),
where the error arises from the quadratic term in the Taylor series.
Although it seems as if we need to know the value of \(y_{n+1/2}\), appearing in \(f\) in the right-hand side of this equation for it to be of any use, this is not quite true:
Since the error is already \(\mathcal{O}(h^3)\), an approximation to \(y_{n+1}\) whose error is \(\mathcal{O}(h^2)\) is good enough. This is provided by the simple Euler’s method:
\(y_{n+1/2} = y_n + \frac{h}{2} f(x_n, y_n) + \mathcal{O}(h^2)\),
and thus, if we define \(k_1 = h f(x_n, y_n)\), we can write:
\(y_{n+1} = y_n + h f ( x_n + \frac{h}{2}, y_n + \frac{1}{2} k_1 ) + \mathcal{O}(h^3)\).
This is the second-order Runge-Kutta algorithm. It requires the evaluation of \(f\) twice for each step.
A fourth-order algorithm, which requires \(f\) to be evaluated four times for each integration step and has a local accuracy of \(\mathcal{O}(h^5)\) has been found by experience to give the best balance between accuracy and computational effort. It can be written as follows, with the \(k_i\)s as intermediate variables:
\(k_1 = h f(x_n, y_n)\),
\(k_2 = h f(x_n + \frac{1}{2} h, y_n + \frac{1}{2} k_1)\),
\(k_3 = h f(x_n + \frac{1}{2} h, y_n + \frac{1}{2} k_2)\),
\(k_4 = h f(x_n + h, y_n + k_3)\),
\(y_{n+1} = y_n + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) + \mathcal{O}(h^5)\).
7.4.2.1. Example 7.3: Try out the second- and fourth-order Runge-Kutta methods on the problem defined in Example 7.1.#
Compare the accuracy of the methods.
7.4.3. Application: Nonlinear Oscillators#
With several algorithms for solving ODEs in our possession, we can now examine the behavior of oscillators beyond the harmonic approximation.
Newton’s second law provides the equation of motion for an oscillator, which could be, e.g. a mass attached to a spring moving in one dimension, \(x\):
\(F_k(x) + F_\mathrm{ext}(x,t) = m \frac{ \mathrm{d}^2 x}{\mathrm{d}t^2}\),
where \(F_k(x)\) is the restoring force, e.g. exerted by a spring, and \(F_\mathrm{ext}(x,t)\) is an external (driving force), that may also depend on time.
Let’s examine two models that can describe the departure from linearity, with no external forces being applied (\(F_\mathrm{ext} = 0\)).
For the first model, let’s introduce a potential that is a harmonic oscillator for small displacements \(x\), but also contains a perturbation that introduces a nonlinear term to the force for large values of \(x\):
\(V_k(x) \simeq \frac{1}{2} kx^2 (1 - \frac{2}{3} \alpha x)\).
Taking he derivative, we can derive \(F_k(x)\):
\(F_k(x) = - \frac{ \mathrm{d} V}{\mathrm{d}x} = -kx ( 1 - \alpha x)\),
and the ODE that we need to consider has the form:
\(m \frac{ \mathrm{d}^2 x}{\mathrm{d}t^2} = -kx (1 - \alpha x)\).
For \(\alpha x \ll 1\), we recover simple harmonic motion, but as \(x \rightarrow 1/\alpha\), the anharmonic effects would be large.
As long as \(x < 1/\alpha\), there will be a restoring force, and the motion will be periodic. In general, there will be an asymmetry in the motion to the right and left of the equilibrium position.
As a second model for nonlinear oscillators, we will consider that the spring’s potential function is proportional to some arbitrary even power \(p\) of the displacement from equilibrium:
\(V(x) = \frac{1}{p} kx^p\), (\(p\) even),
such that the force is a restoring force:
\(F_k(x) = - \frac{ \mathrm{d} V}{\mathrm{d}x} = - k x^{p-1}\).
For \(p=2\), we recover the harmonic oscillator. As \(p\) increases, the potential becomes more and more like a square well, where the mass almost moves freely inbetween “collisions” with the wall at \(x\simeq \pm 1\). Regardless of the value of \(p\), the motion will be periodic, but will only be harmonic for \(p=2\).
Newton’s law gives the second-order ODE that we need to solve:
\(m \frac{ \mathrm{d}^2 x}{\mathrm{d}t^2} = - k x^{p-1}\).
Following our discussion on the numerical solution of ODEs, we can reduce the second-order ODE into two coupled first-order equations via an auxiliary variable, in this case the momentum \(p\):
\( \frac{ \mathrm{d}x}{\mathrm{d}t} = \frac{p}{m}\),
\( \frac{ \mathrm{d}p}{\mathrm{d}t} = F_k(x)\).
Let’s begin by examining simple harmonic motion, that we already know a lot about, before we consider the nonlinear cases.
7.4.3.1. Example 7.4: Oscillator warmup: Simple Harmonic Motion#
(a) Generalize the single-variable methods that we have derived to the two-variable case relevant for oscillatory motion.
(b) Integrate the equations of motion for single harmonic motion, i.e. with \(F_k(x) = -kx\), for values of \(k=4\pi^2\) and \(m=1\) and verify that the period \(T\) is what you expect from your Physics I course (recall: \(T=2\pi \sqrt{m/k}\)). You may choose the initial conditions (e.g. \(p=0\) and \(x=x_0\) at \(t=t_0\), where \(x_0\) would be the amplitude).
We are now in a position to investigate nonlinear oscillations. Let’s consider our two models in the next example.
7.4.3.2. Example 7.5: Nonlinear Oscillator Models#
For this example, we will use our generalized rk4 function from Example 7.4.
Consider the restoring forces arising from both potentials:
\(F_k(x) = - kx ( 1 - \alpha x)\),
and:
\(F_k(x) = - k x^{p-1}\).
Consider values of \(p = 6\) and one anhmaronic strength in \( 0 \leq \alpha x \leq 2\).
(a) Check that the solution remains periodic with constant amplitude for a given initial condition, regardless of how nonlinear you make the force. In particular, check that the maximum speed occurs at \(x=0\), and that the velocity \(v=0\) at the maximum \(x\), the latter being a consequence of energy conservation.
(b) Verify that nonharmonic oscillators are nonisosynchronous, i.e. that vibrations with different amplitudes have different periods.
(c) Devise an algorithm to determinet he period \(T\) of the oscillation, by recording times at which the mass passes through the origin. Note that because the motion may be asymmetric, you must record at least three times to deduce the period.
(d) Construct a graph of the deduced period as a function of initial amplitude.