Intro for those who've never had DEs

Q: What is a DE?

A: A DE is a differential equation, tying a function to its derivative(s). An ODE is an ordinary DE: it links function f(x) to its ordinary derivative(s), not any partial derivative(s).

Classic Example of an ODE: Population Growth

Assume that the increase in a population over a given short (instantaneous would be nice!) interval ${\displaystyle \left.h\right.}$is just proportional to the population itself: that is,

${\displaystyle \left.P(t+h)=P(t)+P(t)\cdot h\cdot r\right.}$

Then

${\displaystyle \left.P'(t)\approx {\frac {P(t+h)-P(t)}{h}}={\frac {P(t)\cdot h\cdot r}{h}}=rP(t)\right.}$.

In the limit as we let ${\displaystyle \left.h\rightarrow 0\right.}$, we have

${\displaystyle \left.P'(t)=rP(t)\right.}$.

This is a classic and very important DE, which we now want to solve. When I asked my students how to solve this, they said "Separation of variables." I said "Don't be silly! Use the general method for solving DEs."

They said that they didn't know what that was, so I told them:

The general method for solving DEs is to stare at them until a solution comes to you.


This is much faster than separation of variables: you simply ask yourself "Do I know a function which is essentially its own derivative, up to a constant?" To which I hope you reply to yourself "Of course! An exponential function."

${\displaystyle \left.P(t)=\alpha e^{rt}\right.}$.

We've found an infinite number of solutions. So now the question is, how do we choose between them? Well, we can't in the absence of more information. But if we have the initial value of the population, ${\displaystyle \left.P(0)\right.}$ at time ${\displaystyle \left.t=0\right.}$, then we can give the exact and unique solution,

${\displaystyle \left.P(t)=P(0)e^{rt}\right.}$.

We've solved the initial value problem.

Classic Example: the Simple Pendulum

You might visit this website, for a nice animation of this problem, along with several derivations. Mine is similar.

We'll derive this using a little vector calculus and a little bit of Newtonian physics. This image will do to give you the idea, but I actually think of it a little differently (thinking of the tension on the pendulum arm):

If we take as our origin the point on the ceiling to which the pendulum is attached, then the vector position of the pendulum bob is

${\displaystyle \left.{\underline {p}}(t)=L\langle \sin(\theta (t)),-\cos(\theta (t))\rangle \right.}$

In the following, I'll suppress the dependence of ${\displaystyle \left.\theta \right.}$ on time, replacing ${\displaystyle \left.\theta (t)\right.}$ with ${\displaystyle \left.\theta \right.}$ (for clarity).

${\displaystyle \left.{\underline {v}}(t)=L{\frac {d\theta }{dt}}\langle \cos(\theta ),\sin(\theta )\rangle \right.}$

${\displaystyle \left.{\underline {a}}(t)=L\left[{\frac {d^{2}\theta }{dt^{2}}}\langle \cos(\theta ),\sin(\theta )\rangle +\left({\frac {d\theta }{dt}}\right)^{2}\langle -\sin(\theta ),\cos(\theta )\rangle \right]\right.}$ or ${\displaystyle \left.{\underline {a}}(t)=L\langle {\frac {d^{2}\theta }{dt^{2}}}\cos(\theta )-\left({\frac {d\theta }{dt}}\right)^{2}\sin(\theta ),{\frac {d^{2}\theta }{dt^{2}}}\sin(\theta )+\left({\frac {d\theta }{dt}}\right)^{2}\cos(\theta )\rangle \right.}$

Now we turn to the physics: Newton's second law is ${\displaystyle \left.{\underline {F}}=m{\underline {a}}\right.}$.

The forces on the pendulum bob are gravity, operating in the negative y direction (of size mg, where m is the mass and g the acceleration due to gravity), and the tension (of size T) on the bob from the pendulum arm. The unbalanced forces give rise to accelerations. If we take the bob at rest, and move it to the right, then we can write the following for the force:

${\displaystyle \left.{\underline {F}}=\langle -T\sin(\theta ),T\cos(\theta )-mg\rangle \right.}$

where ${\displaystyle \left.\theta \right.}$ is the displacement in angle. We equate the force and mass times acceleration and simplify to get ${\displaystyle \left.\langle -T\sin(\theta ),T\cos(\theta )-mg\rangle =mL\langle {\frac {d^{2}\theta }{dt^{2}}}\cos(\theta )-\left({\frac {d\theta }{dt}}\right)^{2}\sin(\theta ),{\frac {d^{2}\theta }{dt^{2}}}\sin(\theta )+\left({\frac {d\theta }{dt}}\right)^{2}\cos(\theta )\rangle \right.}$

Equating components, we get

${\displaystyle \left.{\frac {d^{2}\theta }{dt^{2}}}\cos(\theta )-\left({\frac {d\theta }{dt}}\right)^{2}\sin(\theta )=-{\frac {T}{mL}}\sin(\theta )\right.}$

and

${\displaystyle \left.{\frac {d^{2}\theta }{dt^{2}}}\sin(\theta )+\left({\frac {d\theta }{dt}}\right)^{2}\cos(\theta )={\frac {T}{mL}}\cos(\theta )-{\frac {g}{L}}\right.}$

Our objective is to eliminate the term ${\displaystyle \left.\left({\frac {d\theta }{dt}}\right)^{2}\right.}$ from these two equations. We can do this by multiplying each by a trig function:

${\displaystyle \left.{\frac {d^{2}\theta }{dt^{2}}}\cos ^{2}(\theta )-\left({\frac {d\theta }{dt}}\right)^{2}\sin(\theta )\cos(\theta )=-{\frac {T}{mL}}\sin(\theta )\cos(\theta )\right.}$

and

${\displaystyle \left.{\frac {d^{2}\theta }{dt^{2}}}\sin ^{2}(\theta )+\left({\frac {d\theta }{dt}}\right)^{2}\sin(\theta )\cos(\theta )={\frac {T}{mL}}\sin(\theta )\cos(\theta )-{\frac {g}{L}}\sin(\theta )\right.}$

Adding these two equations leads to magical results: we get elimination of the tension, components, too, and the result is (almost) our lovely ODE:

${\displaystyle \left.{\frac {d^{2}\theta }{dt^{2}}}\left(\sin ^{2}(\theta )+\cos ^{2}(\theta )\right)=-{\frac {g}{L}}\sin(\theta )\right.}$

or, invoking Pythagoras,

${\displaystyle \left.{\frac {d^{2}\theta }{dt^{2}}}=-{\frac {g}{L}}\sin(\theta )\right.}$

Strategies for Solution

The pendulum ODE is interesting, because our method of solution (staring) doesn't immediately give rise to a solution (at least I didn't stare one down). There are two approaches you can take (at least), when staring doesn't work (and so you give up on the exact solution):

1. Replace the ODE with a related one that you like better, or
2. Solve the one you've got numerically.

In either event, you're not actually solving the given ODE -- you're just hoping that your approximate method is giving you an estimate that is close enough for your purposes.

Linearize the ODE

As an example, using our pendulum problem, we could assume that ${\displaystyle \left.\theta \right.}$ is small, and hence that (using Taylor series) ${\displaystyle \left.\sin(\theta )\approx \theta \right.}$ (with error ${\displaystyle \left.O(\theta ^{3})\right.}$). Then our ODE

${\displaystyle \left.{\frac {d^{2}\theta }{dt^{2}}}=-{\frac {g}{L}}\sin(\theta )\right.}$

is approximately given by

${\displaystyle \left.{\frac {d^{2}\theta }{dt^{2}}}=-{\frac {g}{L}}\theta \right.}$

which I can stare down: solutions are sines and cosines. In fact, the general solution is

${\displaystyle \left.\theta (t)=\alpha \sin(\omega t)+\beta \cos(\omega t)\right.}$,

where ${\displaystyle \left.\omega ={\sqrt {\frac {g}{L}}}\right.}$.

This is a linear combination of the solutions we know, and works because the equation is linear. One way to know that you have a linear equation is if the linear combination of known solutions is another solution!

We may hope will be a good approximation so long as ${\displaystyle \left.\theta \right.}$ is small (which means that you don't pull the bob too far from equilibrium -- its rest position).

The problem becomes an IVP, and we find a particular solution, when we know two things about our problem at ${\displaystyle \left.t=0\right.}$ (usually the position and velocity). Two constraints allow us to find the two unknown coefficients ${\displaystyle \left.\alpha \right.}$ and ${\displaystyle \left.\beta \right.}$ (the general rule being that you need one constraint for each unknown).

Use Numerics

Alternatively we can use numerics. But in the process of approximating derivatives, etc., we're also not solving the given problem. On the other hand, we don't have to be nearly as clever as we were in using the Taylor series.

One additional problem in the pendulum problem is that we don't have a first-order problem (it's not about derivatives -- it's about second derivatives!). In this case, we can be clever and use a trick to turn the second-order ODE into a first-order system:

Let

${\displaystyle \left.\gamma (t)={\frac {d\theta }{dt}}\right.}$

Then

${\displaystyle \left.{\frac {d\gamma }{dt}}={\frac {d^{2}\theta }{dt^{2}}}\right.}$

Hence we have a system

${\displaystyle \left.{\begin{bmatrix}{\frac {d\theta }{dt}}\\{\frac {d\gamma }{dt}}\end{bmatrix}}={\frac {d}{dt}}{\begin{bmatrix}\theta \\\gamma \end{bmatrix}}={\begin{bmatrix}\gamma \\-{\frac {g}{L}}\sin(\theta )\end{bmatrix}}\right.}$

This type of non-linear first-order ODE will be the typical subject of our numerical attacks.

Elementary Theory of Initial-Value Problems

The Big Picture

First of all we need to say a little about what differential equations and initial-value problems are, and the conditions under which they have solutions, and when those solutions are unique (it's nice to know that what you're hunting for is out there, and that it's alone).

Then we need to realize that, because we're solving these problems numerically, we're not going to be solving the given initial-value problem, but rather one close to that given (we call that a "perturbed" problem -- it's been jiggled, but hopefully not too badly).

So here is a brief overview of conditions under which solutions (or approximate solutions) of a perturbed problem will give reasonable information about the real problem.

Definitions

• differential equation: an equation linking a function ${\displaystyle y}$ and its derivatives and independent variables.
• initial-value problem: with time as the independent variable, a differential equation as well as enough initial conditions to uniquely determine the solution.
• well-posed initial-value problem: the initial-value problem

${\displaystyle \left(1\right)\qquad \qquad \qquad \qquad \qquad {\frac {dy}{dt}}=f(t,y),\quad a\leq t\leq b,\quad y(a)=\alpha }$

is well-posed if

• A unique solution ${\displaystyle y(t)}$ exists, and
• For any ${\displaystyle \left.\epsilon >0\right.}$ there exists a positive constant ${\displaystyle \left.k(\epsilon )\right.}$ such that whenever ${\displaystyle \left.|\epsilon _{0}|<\epsilon \right.}$ and ${\displaystyle \left.\delta (t)\right.}$ is continuous with ${\displaystyle |\delta (t)|<\epsilon }$ on ${\displaystyle [a,b]}$, a unique solution ${\displaystyle z(t)}$ to

${\displaystyle \left(2\right)\qquad \qquad \qquad \qquad \qquad {\frac {dz}{dt}}=f(t,z)+\delta (t),\quad a\leq t\leq b,\quad y(a)=\alpha +\epsilon _{0}}$

exists with
${\displaystyle \left.|z(t)-y(t)|

The problem (2) is called a perturbed problem associated with the given problem.

• Lipschitz condition: A function ${\displaystyle \left.f(t,y)\right.}$ satisfies a Lipschitz condition in ${\displaystyle \left.y\right.}$ on a set ${\displaystyle \left.D\subset \mathbb {R} ^{2}\right.}$ if a constant ${\displaystyle \left.L>0\right.}$ exists with
${\displaystyle \left.|f(t,y_{1})-f(t,y_{2})|\leq L|y_{1}-y_{2}|\right.}$
• A set ${\displaystyle \left.D\subset \mathbb {R} ^{2}\right.}$ is said to be convex if whenever ${\displaystyle \left.(t_{1},y_{1})\right.}$ and ${\displaystyle \left.(t_{2},y_{2})\right.}$ belong to ${\displaystyle \left.D\right.}$ and ${\displaystyle \left.\lambda \right.}$ is in ${\displaystyle \left.[0,1]\right.}$, the point ${\displaystyle \left.((1-\lambda )t_{1}+\lambda t_{2},(1-\lambda )y_{1}+\lambda y_{2})\right.}$ also belongs to ${\displaystyle \left.D\right.}$.

Theorems/Formulas

Theorem 5.4: Suppose that ${\displaystyle \left.D=\{(t,y)|a\leq t\leq b,-\infty and that ${\displaystyle \left.f(t,y)\right.}$ is continuous on ${\displaystyle \left.D\right.}$. If ${\displaystyle \left.f\right.}$ satisfies a Lipschitz condition on ${\displaystyle \left.D\right.}$ in the variable ${\displaystyle \left.y\right.}$, then the initial value problem has a unique solution ${\displaystyle \left.y(t)\right.}$ for ${\displaystyle \left.a\leq t\leq b\right.}$.

Theorem 5.3: If ${\displaystyle \left.f(t,y)\right.}$ is defined on a convex set ${\displaystyle \left.D\subset \mathbb {R} ^{2}\right.}$. If a constant ${\displaystyle \left.L>0\right.}$ exists with

${\displaystyle \left|{\frac {\partial f}{\partial y}}(t,y)\right|\leq L\quad \forall (t,y)\in D,}$

then ${\displaystyle \left.f\right.}$ satisfies a Lipschitz condition in ${\displaystyle \left.y\right.}$ on ${\displaystyle \left.D\right.}$ with Lipschitz constant ${\displaystyle \left.L\right.}$.

Theorem 5.6: Suppose that ${\displaystyle \left.D=\{(t,y)|a\leq t\leq b,-\infty and that ${\displaystyle \left.f(t,y)\right.}$ is continuous on ${\displaystyle \left.D\right.}$. If ${\displaystyle \left.f\right.}$ satisfies a Lipschitz condition in ${\displaystyle \left.y\right.}$ on ${\displaystyle \left.D\right.}$, then the initial value problem (1) is well-posed.

Properties/Tricks/Hints/Etc.

So the up-shot is that our initial-value problem satisfies a Lipschitz condition, we're in good shape: we're going to be able to get an approximation numerically.