Gauss–Legendre Runge–Kutta methods

Remarks

The method is the collocation method at the nodes $c_i$: it finds the polynomial $u$ of degree $s$ with $u(t_n) = y_n$ and $u'(t_n + c_i h) = f(t_n + c_i h, u(t_n + c_i h))$ for each $i$, and sets $y_{n+1} = u(t_n + h)$. The weights $b_j$ and nodes $c_j$ form the Gauss–Legendre quadrature rule on $[0, 1]$, which is exact up to degree $2s - 1$; collocation theory turns this into order $2s$. Example: $s = 1$ gives $c_1 = \frac12$, $a_{11} = \frac12$, $b_1 = 1$, the implicit midpoint method $y_{n+1} = y_n + h f\big(t_n + \frac h2, \frac{y_n + y_{n+1}}{2}\big)$, with $R(z) = \frac{1 + z/2}{1 - z/2}$; clause 3 reads $\frac12 + \frac12 = 1$. Example: $s = 2$ gives $c = \frac12 \mp \frac{\sqrt3}{6}$, $b = (\frac12, \frac12)$, $a_{11} = a_{22} = \frac14$, $a_{12} = \frac14 - \frac{\sqrt3}{6}$, $a_{21} = \frac14 + \frac{\sqrt3}{6}$, a method of order $4$. The price of these properties is that each step needs the solution of a nonlinear system of dimension $sd$ for a system of dimension $d$.

Going deeper

For a student

The family links three topics: Gaussian quadrature, which fixes the nodes; collocation, which gives the order; and geometric integration, through the symplecticity condition. For Hamiltonian systems they are standard high-order symplectic methods, used for long-time simulations in celestial mechanics and molecular dynamics. Their stability function satisfies $\abs{R(iy)} = 1$ on the imaginary axis and $\abs{R(z)} \to 1$ as $z \to -\infty$, so unlike backward Euler they do not damp very stiff components.