Numerical Methods

When analytical solutions are not possible, we can use numerical methods to approximate the solution of an ODE. This is common because the ODEs are derived from physical laws. Given any differential equation, for example,

$$\stackrel{...}{x}+\ddot{x}\ddot{x}+\sin (\ddot{x})=1$$

we can convert to a first-order system (involving at most the first derivative). Let velocity $v=\dot{x}$ and acceleration $a=\ddot{x}$. Then we have

$$\{\begin{array}{l}\dot{x}=v\\ \dot{v}=a\\ \dot{a}=1-av-\sin (a)\end{array}$$

Let $y=[x,v,a{]}^{\top }$. The ODE becomes $\dot{y}=\mathbf{f}(y)$.

Forward Euler Method

We can discretize time into time-frames.

$$y^{(n)}=y(n\Delta t)$$

Then, we calculate the next state $y^{(n+1)}$ by using the current state $y^{(n)}$ and the derivative ${\dot{y}}^{(n)}$. The general idea is that

$$\frac{y^{(n+1)}-y^{(n)}}{\Delta t}\approx \mathbf{f}(y^{(n)})$$

such that

$$y^{(n+1)}\approx y^{(n)}+\Delta t\cdot \mathbf{f}(y^{(n)})$$

Recall that $\mathbf{f}$ is a vector-valued function that describes the change in the state $y$ at time $t$.

Advantages

• There is an explicit formula formula for the next state.
• Computationally fast and easy to implement. Limitations
• Not very accurate, unless $\Delta t$ is tiny. The error is proportional to $\Delta t^2$ from the Taylor expansion of $y(t+\Delta t)$.
• Can be energy increasing (violates conservation of energy).

Backward Euler Method

The backward Euler method is an implicit method that calculates the next state $y^{(n+1)}$ using the derivative at the next state rather than the current state. The formula is given by

$$\frac{y^{(n+1)}-y^{(n)}}{\Delta t}\approx \mathbf{f}(y^{(n+1)})$$

Advantages

• Energy decreasing (dissipating). This “looks” physical, i.e. it is more realistic as we expect energy to dissipate over time.
• Can take large time steps without instability.
• Can incorporate collision.

Limitations

• Not very accurate, unless $\Delta t$ is tiny.
• We have to solve implicitly for $y^{(n+1)}$ at each step, which can be computationally expensive.

Runge-Kutta Method (RK4)

This method is accurate, stable, and explicit.

$$\begin{array}{rl}{k}_1& =\mathbf{f}(y^{(n)})\\ k_2& =\mathbf{f}\left(y^{(n)}+\frac{\Delta t}{2}{k}_1\right)\\ k_3& =\mathbf{f}\left(y^{(n)}+\frac{\Delta t}{2}{k}_2\right)\\ k_4& =\mathbf{f}\left(y^{(n)}+\Delta t{k}_3\right)\\ y^{(n+1)}& =y^{(n)}+\frac{\Delta t}{6}(k_1+2k_2+2k_3+k_4)\end{array}$$

The idea is that we take a weighted average of the slopes at different points within the time step to get a better estimate of the next state. In most cases, the RK4 method works very well. There are other special algorithms (non-RK4) that are deigned to preserve energy or momentum.

• Variational integrator
• Symplectic integrator
• Lie group integrator

Example 1

Consider the pendulum equation $m\ddot{\theta }=-mg\sin \theta$. Given $({\theta }_i,{\dot{\theta }}_i)$ at time $t_i$, we can use the RK4 method to compute $({\theta }_{i+1},{\dot{\theta }}_{i+1})$ at time $t_{i+1}=t_i+\Delta t$.

Recall from example that

$$f(\theta ,v)=\dot{y}=\left[\begin{array}{c}v\\ -\sin \theta \end{array}\right]$$

(we omit $g$ for simplicity). Then, we can compute the intermediate slopes as follows:

$$k_1=f({\theta }_i,v_i)=\left[\begin{array}{c}{v}_i\\ -\sin {\theta }_i\end{array}\right]$$

To prepare for $k_2$, we need to use $k_1$ to project halfway forward $\frac{\Delta t}{2}$. This gives

$$\left[\begin{array}{c}{\theta }_{i+1/2}^{\ast }\\ v_{i+1/2}^{\ast }\end{array}\right]=\left[\begin{array}{c}{\theta }_i\\ v_i\end{array}\right]+\frac{\Delta t}{2}\left[\begin{array}{c}{v}_i\\ -\sin {\theta }_i\end{array}\right]$$

such that

$$\begin{array}{rl}{\theta }_{i+1/2}^{\ast }& ={\theta }_i+\frac{\Delta t}{2}{v}_i\\ v_{i+1/2}^{\ast }& =v_i-\frac{\Delta t}{2}\sin {\theta }_i\end{array}$$

Then, we can compute $k_2$ using the projected state:

$$k_2=f({\theta }_{i+1/2}^{\ast },v_{i+1/2}^{\ast })=\left[\begin{array}{c}{v}_{i+1/2}^{\ast }\\ -\sin {\theta }_{i+1/2}^{\ast }\end{array}\right]$$

Next, we can compute $k_3$ by projecting forward again using $k_2$:

$$\left[\begin{array}{c}{\theta }_{i+1/2}^{\ast \ast }\\ v_{i+1/2}^{\ast \ast }\end{array}\right]=\left[\begin{array}{c}{\theta }_i\\ v_i\end{array}\right]+\frac{\Delta t}{2}\left[\begin{array}{c}{v}_{i+1/2}^{\ast }\\ -\sin {\theta }_{i+1/2}^{\ast }\end{array}\right]$$

such that

$$\begin{array}{rl}{\theta }_{i+1/2}^{\ast \ast }& ={\theta }_i+\frac{\Delta t}{2}{v}_{i+1/2}^{\ast }\\ v_{i+1/2}^{\ast \ast }& =v_i-\frac{\Delta t}{2}\sin {\theta }_{i+1/2}^{\ast }\end{array}$$

Then we can compute $k_3$ using the projected state:

$$k_3=f({\theta }_{i+1/2}^{\ast \ast },v_{i+1/2}^{\ast \ast })=\left[\begin{array}{c}{v}_{i+1/2}^{\ast \ast }\\ -\sin {\theta }_{i+1/2}^{\ast \ast }\end{array}\right]$$

Finally, we can compute $k_4$ by projecting forward using $k_3$:

$$\left[\begin{array}{c}{\theta }_{i+1}^{\ast \ast \ast }\\ v_{i+1}^{\ast \ast \ast }\end{array}\right]=\left[\begin{array}{c}{\theta }_i\\ v_i\end{array}\right]+\frac{\Delta t}{2}\left[\begin{array}{c}{v}_{i+1/2}^{\ast \ast }\\ -\sin {\theta }_{i+1/2}^{\ast \ast }\end{array}\right]$$

such that

$$\begin{array}{rl}{\theta }_{i+1}^{\ast \ast \ast }& ={\theta }_i+\Delta t{v}_{i+1/2}^{\ast \ast }\\ v_{i+1}^{\ast \ast \ast }& =v_i-\Delta t\sin {\theta }_{i+1/2}^{\ast \ast }\end{array}$$

Then, we can compute $k_4$ using the projected state:

$$k_4=f({\theta }_{i+1}^{\ast \ast \ast },v_{i+1}^{\ast \ast \ast })=\left[\begin{array}{c}{v}_{i+1}^{\ast \ast \ast }\\ -\sin {\theta }_{i+1}^{\ast \ast \ast }\end{array}\right]$$

Plugging everything into the RK4 formula, we get

$$\begin{array}{rl}\left[\begin{array}{c}{\theta }_{i+1}\\ v_{i+1}\end{array}\right]& =\left[\begin{array}{c}{\theta }_i\\ v_i\end{array}\right]+\frac{\Delta t}{6}\left(k_1+2k_2+2k_3+k_4\right)\\ & =\left[\begin{array}{c}{\theta }_i\\ v_i\end{array}\right]+\frac{\Delta t}{6}\left(\left[\begin{array}{c}{v}_i\\ -\sin {\theta }_i\end{array}\right]+2\left[\begin{array}{c}{v}_{i+1/2}^{\ast }\\ -\sin {\theta }_{i+1/2}^{\ast }\end{array}\right]+2\left[\begin{array}{c}{v}_{i+1/2}^{\ast \ast }\\ -\sin {\theta }_{i+1/2}^{\ast \ast }\end{array}\right]+\left[\begin{array}{c}{v}_{i+1}^{\ast \ast \ast }\\ -\sin {\theta }_{i+1}^{\ast \ast \ast }\end{array}\right]\right)\end{array}$$

The choice of $\Delta t$ is crucial for the accuracy and stability of the numerical solution. Graphically, if $\Delta t=0.1$, the numerical solution follows the true solution closely. However, if $\Delta t=0.9$, we preserve no information about the true solution and the numerical solution diverges significantly from the true solution.

Example 1.2 (2nd Order Discretization)

We can also use a second-order discretization method to solve the pendulum equation, or central difference approximation for a second derivative (in our case, acceleration). We consider

$$\begin{array}{rl}{v}_{\text{forward}}& =\frac{{\theta }_{i+1}-{\theta }_i}{\Delta t}\\ v_{\text{backward}}& =\frac{{\theta }_i-{\theta }_{i-1}}{\Delta t}\end{array}$$

and take the average to find the acceleration:

$$\begin{array}{rl}{a}_i& =\frac{v_{\text{forward}}-v_{\text{backward}}}{\Delta t}\\ & =\frac{{\theta }_{i+1}-2{\theta }_i+{\theta }_{i-1}}{\Delta t^2}\end{array}$$

But recall $\ddot{\theta }=a$, so we can rearrange to get

$$\frac{{\theta }_{i+1}-2{\theta }_i+{\theta }_{i-1}}{\Delta t^2}=-\sin {\theta }_i$$