Ordinary Differential Equations (ODEs)

Basic Concepts

ODE depends only on one variable and has the structure: \begin{equation} \frac{dy}{dx}=f(x,y) \text{ or } \frac{d^2 y}{dx^2} = f(x,y,\frac{dy}{dx}) \end{equation} A second order ODE can be transformed into a coupled set of first order of ODEs: \begin{equation} p = \frac{dy}{dx};\frac{dp}{dx}=f(x,y,p) \end{equation}

Numerical Methods

1. Euler Algorithm

Approximate derivative by forward difference: \begin{equation} y_{n+1}=y_n+hf(x_n,y_n)+O(h^2) \end{equation} Euler method has mainly educational value. It is not very accurate, and encounters stability problems.

2. Verlet Algorithm

Designed for Newton’s Second Law \(\frac{d^2x}{dt^2}=a(x,t)\). Use the central difference to eliminate the odd terms, keeping the time reversal symmetry. \begin{equation} x_{n+1}=2x_n-x_{n-1}+h^2a_n \end{equation}

3. Leapfrog Algorithm

Leapfrog Algorithm implies the evaluation of position and velocity at different time steps. The algorithm is given by the following steps:

First calculate the velocity half step later by current acceleration. \begin{equation} x^{(1)}(t+h/2)=x^{(1)}(t-h/2)+hx^{(2)}(t)+O(h^2) \end{equation}
Calculate the position at next step. \begin{equation} x(t+h)=x(t)+hx^{(1)}(t+h/2)+O(h^3) \end{equation}
Calculate current acceleration. \begin{equation} x^{(1)}(t+h)=(x^{(1)}(t+h/2)+\frac{h}{2}x^{(2)}(t+h))+O(h^2) \end{equation}

4. Predictor-Corrector Mtehods

In addition to Euler’s method, we can compute the slope \(k_2\) at \(t_{i+1}\) with the result of Euler’s method \(y_{i+1}\) (prediction). Then take the average of two slopes and compute again with: \begin{equation} y_{i+1} \approx y(t_i)+\frac{h}{2} (k_1+k_2) \end{equation}

5. Runge-Kutta Methods

Runge-Kutta methods use a trial step at the midpoint of an interval to cancel out lower-order error terms.

Second-Order Formula

\begin{equation} k_1 = hf(t_i,y_i); k_2 = hf(t_{i+1/2},y_i+k_1/2) \end{equation} \begin{equation} y_{i+1} \approx y_i+k_2+O(h^3) \end{equation}

Fourth-Order Formula

We compute first \(k_1=hf(t_i,y_i)\), which is nothing but the slope at \(t_i\) (Euler’s method).
Then we compute the slope at the midpoint \(t_{i}+h/2\) using Euler’s prediction, as in the second-order RK: \begin{equation} k_2 = hf(t_i+h/2,y_i+k_i/2). \end{equation}
The improved slope at the midpoint is used to further improve the slope of \(y_{i+1/2}\) by computing \begin{equation} k_3 = hf(t_i+h/2,y_i+k_2/2). \end{equation}
With \(k_3\) we can predice \(y_{i+1}=y_i+k_3\) and compute the slope at \(t_i+h\): \begin{equation} k_4 = hf(t_i+h,y_i+k_3). \end{equation}
The final algorithm becomes then \begin{equation} y_{i+1}=y_i+\frac{1}{6}(k_1+2k_2+2k_3+k_4). \end{equation}

Summery

In the case of energy-conserving systems (no damping or external driving forces):

Use the Verlet/leapfrog algorithm
In the case of non-energy-conserving systems (including damping or external driving forces)
Use the RK algorithm

Second-Order ODE Problems

Boundary-Value Problems (BVP)

In this case of problems, we can not use step algorithms directly because only know \(u(0)\) or \(u'(0)\). \begin{equation} \frac{d^2u}{dx^2}=f(x,u,u’), x\in [x_{lb},x_{rb}] \end{equation}

Boundary Conditions

Dirichlet: \(u(x_{lb}),u(x_{rb})\) provided.
Neumann: \(u'(x_{lb}),u'(x_{rb})\) provided.
Periodic: \(u(x_{lb})=u(x_{rb})\) and \(u'(x_{lb})=u'(x_{rb})\).
Mixed: combination of Dirichlet and Neumann.

Eigenvalue Problems

Solution that satisfies BC only exists for specific \(\lambda\)(eigenvalue). For example (1-D TISE): \begin{equation} \frac{d^2\psi}{dx^2}=-\frac{2m}{\hbar^2}[E-V(x)]\psi(x) \end{equation}

Methods

Shooting Method

Shooting method transfers BVP into IVP+root search problem:

Assume unknown initial parameter. (e.g. \(u'(0)=\alpha\))
Compute trial solution \(u_\alpha(1)\) with step forward algorithm (RK4).
Define target function: \(F(\alpha)=u_\alpha(1)-u(1)\).
Search root \(\alpha\) s.t. \(F(\alpha)=0\).

Superposition Shooting Method

If linear equation \(u''+d(x)u'+q(x)u=s(x)\): We need only two trial solutions \(u_{\alpha_0}\) and \(u_{\alpha_1}\). The correct solution of the equation is \begin{equation} u(x)=au_{\alpha_0}(x)+bu_{\alpha_1}(x), \end{equation} where \(a\) and \(b\) are determined from \(u(0)=u_0\) and \(u(1)=u_1\): \begin{equation} a=\frac{u_{\alpha_1}(1)-u_1}{u_{\alpha_1}(1)-u_{\alpha_0}(1)};b=\frac{u_1-u_{\alpha_0}(1)}{u_{\alpha_1}(1)-u_{\alpha_0}(1)} \end{equation}

Numerov’s Algorithm

Consider the equation \begin{equation} \frac{d^2u}{dx^2}+q(x)u(x)=s(x), \end{equation} which does not have a first-order derivative term.

From three-point difference we have the formula \begin{equation} (1+\frac{h^2}{12}q_{n+1})u_{n+1}-2(1+\frac{5h^2}{12}q_n)u_n+(1+\frac{h^2}{12}q_{n-1})u_{n-1}=\frac{h^2}{12}(s_{n+1}+10s_n+s_{n-1})+O(h^6). \end{equation}