Gauss–Legendre method

Gauss-Legendre Runge-Kutta (GLRK) methods solve an ordinary differential equation $\backslash dot\{x\}\; =\; f(x)$ with $x(0)\; =\; x\_0$. The distinguishing feature of GLRK is the estimation of $x(h)\; -\; x\_0\; =\; \backslash int\_0^h\; f(\; x(t)\; )\; \backslash ,\; dt$ with Gaussian quadrature.

$$x(h)\; =\; x(0)\; +\; \backslash frac\{h\}\{2\}\; \backslash sum\_\{i=1\}^\backslash ell\; w\_i\; k\_i\; +\; O(h^\{2\backslash ell\}),$$

where $k\_i\; =\; f(\; x(\; h\; c\_i)\; )$ are the sampled velocities, $w\_i$ are the quadrature weights, $c\_i\; =\; \backslash frac\{1\}\{2\}\; h\; (1+r\_i)$ are the abscissas, and $r\_i$ are the roots $P\_\backslash ell(r\_i)\; =\; 0$ of the Legendre polynomial of degree $\backslash ell$. A further approximation is needed, as $k\_i$ is still impossible to evaluate. To maintain truncation error of order $O(h^\{2\backslash ell\})$, we only need $k\_i$ to order $O(h^\{2\backslash ell-1\})$. The Runge-Kutta implicit definition $k\_i\; =\; f\{\backslash left(x\_0\; +\; h\; \backslash sum\_j\; a\_\{ij\}\; k\_j\; \backslash right)\}$ is invoked to accomplish this. This is an implicit constraint that must be solved by a root finding algorithm like Newton's method. The values of the Runge-Kutta parameters $a\_\{ij\}$ can be determined from a Taylor series expansion in $h$.

The Gauss-Legendre methods are implicit, so in general they cannot be applied exactly. Instead one makes an educated guess of $k\_i$, and then uses Newton's method to converge arbitrarily close to the true solution. Below is a Matlab function which implements the Gauss-Legendre method of order four.

% starting point

x = [ 10.5440; 4.1124; 35.8233];

dt = 0.01;

N = 10000;

x_series = [x];

for i = 1:N

x = gauss_step(x, @lorenz_dynamics, dt, 1e-7, 1, 100);

x_series = [x_series x];

end

plot3( x_series(1,:), x_series(2,:), x_series(3,:) );

set(gca,'xtick',[],'ytick',[],'ztick',[]);

title('Lorenz Attractor');

return;

function [td, j] = lorenz_dynamics(state)

% return a time derivative and a Jacobian of that time derivative

x = state(1);

y = state(2);

z = state(3);

sigma = 10;

beta = 8/3;

rho = 28;

td = [sigma*(y-x); x*(rho-z)-y; x*y-beta*z];

j = [-sigma, sigma, 0;

rho-z, -1, -x;

y, x, -beta];

end

function x_next = gauss_step( x, dynamics, dt, threshold, damping, max_iterations )

[d,~] = size(x);

sq3 = sqrt(3);

if damping > 1 || damping <= 0

error('damping should be between 0 and 1.')

end

% Use explicit Euler steps as initial guesses

[k,~] = dynamics(x);

x1_guess = x + (1/2-sq3/6)*dt*k;

x2_guess = x + (1/2+sq3/6)*dt*k;

[k1,~] = dynamics(x1_guess);

[k2,~] = dynamics(x2_guess);

a11 = 1/4;

a12 = 1/4 - sq3/6;

a21 = 1/4 + sq3/6;

a22 = 1/4;

error = @(k1, k2) [k1 - dynamics(x+(a11*k1+a12*k2)*dt); k2 - dynamics(x+(a21*k1+a22*k2)*dt)];

er = error(k1, k2);

iteration=1;

while (norm(er) > threshold && iteration < max_iterations)

fprintf('Newton iteration %d: error is %f.\n', iteration, norm(er) );

iteration = iteration + 1;

[~, j1] = dynamics(x+(a11*k1+a12*k2)*dt);

[~, j2] = dynamics(x+(a21*k1+a22*k2)*dt);

j = [eye(d) - dt*a11*j1, -dt*a12*j1;

-dt*a21*j2, eye(d) - dt*a22*j2];

k_next = [k1;k2] - damping * linsolve(j, er);

k1 = k_next(1:d);

k2 = k_next(d+(1:d));

er = error(k1, k2);

end

if norm(er) > threshold

error('Newton did not converge by %d iterations.', max_iterations);

end

x_next = x + dt / 2 * (k1 + k2);

end

This algorithm is surprisingly cheap. The error in $k\_i$ can fall below $10^\{-12\}$ in as few as 2 Newton steps. The only extra work compared to explicit Runge-Kutta methods is the computation of the Jacobian.

File:Lorenz Attractor from Gauss-Legendre.png

At the cost of adding an additional implicit relation, these methods can be adapted to have time reversal symmetry. In these methods, the averaged position $(x\_f+x\_i)/2$ is used in computing $k\_i$ instead of just the initial position $x\_i$ in standard Runge-Kutta methods. The method of order 2 is just an implicit midpoint method.

: $k\_1\; =\; f\backslash left(\backslash frac\{x\_f+x\_i\}\{2\}\backslash right)$

: $x\_f\; =\; x\_i\; +\; h\; k\_1$

The method of order 4 with 2 stages is as follows.

: $k\_1\; =\; f\backslash left(\; \backslash frac\{x\_f+x\_i\}\{2\}\; -\; \backslash frac\{\backslash sqrt\{3\}\}\{6\}\; h\; k\_2\backslash right)$

: $k\_2\; =\; f\backslash left(\; \backslash frac\{x\_f+x\_i\}\{2\}\; +\; \backslash frac\{\backslash sqrt\{3\}\}\{6\}\; h\; k\_1\backslash right)$

: $x\_f\; =\; x\_i\; +\; \backslash frac\{h\}\{2\}(k\_1\; +\; k\_2)$

The method of order 6 with 3 stages is as follows.

: $k\_1\; =\; f\backslash left(\; \backslash frac\{x\_f\; +\; x\_i\}\{2\}\; -\; \backslash frac\{\backslash sqrt\{15\}\}\{15\}\; h\; k\_2\; -\; \backslash frac\{\backslash sqrt\{15\}\}\{30\}\; h\; k\_3\; \backslash right)$

: $k\_2\; =\; f\backslash left(\; \backslash frac\{x\_f\; +\; x\_i\}\{2\}\; +\; \backslash frac\{\backslash sqrt\{15\}\}\{24\}\; h\; k\_1\; -\; \backslash frac\{\backslash sqrt\{15\}\}\{24\}\; h\; k\_3\; \backslash right)$

: $k\_3\; =\; f\backslash left(\; \backslash frac\{x\_f\; +\; x\_i\}\{2\}\; +\; \backslash frac\{\backslash sqrt\{15\}\}\{30\}\; h\; k\_1\; +\; \backslash frac\{\backslash sqrt\{15\}\}\{15\}\; h\; k\_2\; \backslash right)$

: $x\_f\; =\; x\_i\; +\; \backslash frac\{h\}\{18\}(\; 5\; k\_1\; +\; 8k\_2\; +\; 5k\_3)$

{{reflist}}

• {{Citation | last1=Iserles | first1=Arieh | author1-link=Arieh Iserles | title=A First Course in the Numerical Analysis of Differential Equations | publisher=Cambridge University Press | isbn=978-0-521-55655-2 | year=1996|url=https://books.google.com/books?id=7Zofw3SFTWIC&q=%22Gauss%E2%80%93Legendre%22}}.

{{DEFAULTSORT:Gauss-Legendre method}}

Category:Runge–Kutta methods