Numerical Integration: 2nd and 4th order Runge-Kutta schemes

Second-order Runge-Kutta scheme (RK2)


There are several variations of the second order Runge-Kutta schemes of numerical integration.

They are predictor-correction methods that it involves 2 stages. A step that predicts the value at the next time step (using the explicit Euler's method as predictor) and another step that applies a correction to the estimated value.

Put simply, it improves upon the Euler method by adding a midpoint in the step which increases the accuracy by one order.

The local error at each step is order \(O(h^3)\) giving a global error of order \(O(h^2)\).

Heun method


The Heun's method is expressed as follows:
\[y_{n+1} = y_n + \frac{h}{2}[f(t_n, y_n) + f(t_{n+1},y_{n+1})]\]
where:
  • \(y_n\) is the approximation of solution at time \(t_n\)
  • \(y_{n+1}\) is the approximation of solution at time \(t_{n+1}\)
  • \(h\) is the step size
  • \(f(t,y)\) is the function that defines the derivative of \(y\) with repsect to \(t\), ie \(\dot y = f(t, y)\)

The two steps are as follows:
  1. In the first step, calculate the estimate of slope of solution at beginning of time step. Then calculate a provisional value of the solution at the end of the time step, by using explicit Euler.
  2. Calculate the estimate of slope of solution at provisional value of the solution. Then use this estimate to calculate a new improved estimate of solution at end of time step.
To make the above single line formula easier to understand, the method may be expressed as follows:
\[y_{n+1}=h\frac{1}{2}\left[k_1+k_2  \right]\]
where:
\begin{aligned} k_{1}&=\ f(t_{n},y_{n}),\\ k_{2}&=\ f(t_{n+1},hf(t_{n+1}, y_{n+1})) = f(t_{n+1},hk_1) \end{aligned}

Written in pseudo computer code, this becomes:
\[\begin{aligned}p_1&=x_n\\v_1&=v_n\\a_1&=a(p_1,v_1)\\p_2&=p_1+v_1\cdot\Delta t\\v_2&=v_1+a_1\cdot \Delta t\\a_2&=a(p_2,v_2)\\x_{n+1}&=x_n+\frac{v_1+v_2}{2}\cdot \Delta t\\v_{n+1}&=v_n+\frac{a_1+a_2}{2}\cdot \Delta t \end{aligned} \]

Comparison to basic Euler


In essence, where as the basic forward Euler method takes the slope of the solution curve at the beginning of the timestep \((k_1)\) and extrapolates the solution at the end of the time step \(t_n\).

The Heun method estimates a better slope by taking the average of the slope at the beginning of the time step and the end of the timestep. However we do not know the slope at \(t_{n+1}\). So we estimate it using basic Euler.


Midpoint method


Explicit Midpoint method


The explicit midpoint method, sometimes known as the modified Euler method is:
\[y_{n+1}=y_n + hf(t_n+h/2, y_n+hf(t_n,y_n))\]
or,
\[y_{n+1}=y_n+hk_2\]
where:
\[\begin{aligned} k_{1}&=f(t_{n},y_{n}), \\  k_{2}&=f(t_{n+\frac{1}{2}},y_{n+\frac{1}{2}}) = f(t_{n+\frac{1}{2}},hk_1)     \end{aligned}\]

Implicit Midpoint method

\[y_{n+1}=y_n + hf(t_{n+\frac{h}{2}}, \frac{1}{2}(y_n+y_{n+1}))\]



Fourth-order Runge Kutta method (RK4)



The "classic Runge-Kutta method" or simply the "Runge-Kutta method" refers to the following.

The formula is given by

\[\begin{aligned} y_{n+1}&=y_{n}+{\frac {h}{6}}\left(k_{1}+2k_{2}+2k_{3}+k_{4}\right),\\ t_{n+1}&=t_{n}+h\\ \end{aligned} \]
for n = 0, 1, 2, 3, ..., using:
\[{\displaystyle {\begin{aligned}k_{1}&=\ f(t_{n},y_{n}),\\k_{2}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{1}}{2}}\right),\\k_{3}&=\ f\!\left(t_{n}+{\frac {h}{2}},y_{n}+h{\frac {k_{2}}{2}}\right),\\k_{4}&=\ f\!\left(t_{n}+h,y_{n}+hk_{3}\right)\end{aligned}}}\]


Here \(y_{n+1}\) is the RK4 approximation of \(y(t_{n+1})\), and the next value \((y_{n+1})\) is determined by the present value \((y_{n})\) plus the weighted average of four increments, where each increment is the product of the size of the interval, h, and an estimated slope specified by function f on the right-hand side of the differential equation.




Useful References


Numerical Integration: Verlet integration


Position Verlet


 To a close approximation (actually performing a Taylor Series Expansion about  x(t±Δt) ), that might look like this:

\[x(t+\Delta t) = x(t) + v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 + \frac{1}{6}b(t) \Delta t^3 + \mathcal{O}(\Delta t^4)\]


This means that if we need to find the next  x, we need the current  x ,  v ,  a , etc. However, because few people calculate the jerk term, our error is typically  \(\mathcal{O}(\Delta t^3)\) . That said, we can calculate  x  with less knowledge and higher accuracy if we play a trick! Let's say we want to calculate  x  of the previous timestep. Again, to a close approximation, that might look like this:

\[x(t-\Delta t) = x(t) - v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 - \frac{1}{6}b(t) \Delta t^3 + \mathcal{O}(\Delta t^4)\]


Now, we have two equations to solve for two different timesteps in x, one of which we already have. If we add the two equations together and solve for x(t+Δt), we find that

\[x(t+ \Delta t) = 2x(t) - x(t-\Delta t) + a(t)\Delta t^2 + \mathcal{O}(\Delta t^4)\]
So, this means we can find our next  x simply by knowing our current  x , the  x  before that, and the acceleration! No velocity necessary! In addition, this drops the error to \(\mathcal{O}(\Delta t^4)\) , which is great! Here is what it looks like in code:


  
  XXXXXXXXXXXXXXXXX





Now, obviously this poses a problem; what if we want to calculate a term that requires velocity, like the kinetic energy, \(\frac{1}{2}mv^2\). In this case, we certainly cannot get rid of the velocity! Well, we can find the velocity to  O(Δt2)  accuracy by using the Stormer-Verlet method. We have the equations for  x(t+Δt)  and  x(t−Δt) above, so let's start there. If we subtract the latter from the former, we get the following:
\[x(t+\Delta t) - x(t - \Delta t) = 2v(t)\Delta t + \frac{1}{3}b(t)\Delta t^3\]
When we solve for  v(t), we get
\[\begin{align} v(t) &= \frac{x(t+\Delta t) - x(t-\Delta t)}{2\Delta t} + \frac{b(t) \Delta t^3}{3 \Delta t} \\  v(t) &= \frac{x(t+\Delta t) - x(t-\Delta t)}{2\Delta t} + \mathcal{O}(\Delta t^2). \end{align}\]

Note that the 2 in the denominator makes sense because we are going over 2 timesteps. It's essentially solving  \(v=\frac{Δx}{Δt}\). In addition, we can calculate the velocity of the next timestep like so
\[v(t+\Delta t) = \frac{x(t+\Delta t) - x(t)}{\Delta t} + \mathcal{O}(\Delta t)\]


Velocity Verlet

Now, let's say we actually need the velocity to calculate out next timestep. Well, in this case, we simply cannot use the above approximation and instead need to use the Velocity Verlet algorithm.

In some ways, this algorithm is even simpler than above. We can calculate everything like

\[\begin{align} x(t+\Delta t) &=x(t) + v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 \\ a(t+\Delta t) &= f(x(t+\Delta t)) \\ v(t+\Delta t) &= v(t) + \frac{1}{2}(a(t) + a(t+\Delta t))\Delta t \end{align}\]
which is literally the kinematic equation above, solving for  x ,  v , and  a  every timestep. You can also split up the equations like so
\[\begin{align} v(t+\frac{1}{2}\Delta t) &= v(t) + \frac{1}{2}a(t)\Delta t \\ x(t+\Delta t) &=x(t) + v(t+\frac{1}{2}\Delta t)\Delta t \\ a(t+\Delta t) &= f(x(t+\Delta t)) \\ v(t+\Delta t) &= v(t+\frac{1}{2}\Delta t) + \frac{1}{2}a(t+\Delta t)\Delta t \end{align}\]

Here is the velocity Verlet method in code:



  XXXXXXXXXXXXXXXXXX


Even though this method is more widely used than the simple Verlet method mentioned above, it unfortunately has an error term of \(\mathcal{O}(\Delta t^2)\), which is two orders of magnitude worse. 


Verlet Leapfrog Algorithm

Takes the form

\[\begin{aligned}x_{i+1}&=x_{i}+v_{i}\Delta t+{\frac {1}{2}}a_{i}\Delta t^{2},\\v_{i+1}&=v_{i}+{\frac {1}{2}}(a_{i}+a_{i+1})\Delta t.\end{aligned}\]


Useful References

https://en.wikipedia.org/wiki/Verlet_integration

https://julesjacobs.com/2019/03/15/leapfrog-verlet.html

https://pikuma.com/blog/verlet-integration-2d-cloth-physics-simulation


https://www.algorithm-archive.org/contents/verlet_integration/verlet_integration.html

Numerical Integration: 2nd and 4th order Runge-Kutta schemes

Second-order Runge-Kutta scheme (RK2) There are several variations of the second order Runge-Kutta schemes of numerical integration. They ar...