Figure 8.3

% Illustration of GL4 collocation approximation

% Butcher tableau for GL4.
tab = struct('a', [1/4, 1/4 - sqrt(3)/6; 1/4 + sqrt(3)/6, 1/4], ...
             'b', [1/2, 1/2], ...
             'c', [1/2 - sqrt(3)/6, 1/2 + sqrt(3)/6]);

% Run one step of GL4.
h = 2*pi()/5;
x0 = [1; 0];
[x1, ~, ~, k] = butcherintegration(x0, @stiffode, h, tab.a, tab.b);

% Get exact solution.
t = linspace(0, h, 1001);
x = [cos(t); -sin(t)];
xdot = [-sin(t); -cos(t)];

% Get collocation polynomial approximation.
k1 = k(:,1);
k2 = k(:,2);
c1 = tab.c(1);
c2 = tab.c(2);

rho1 = h*(k1*c2 - k2*c1)/(c1 - c2);
rho2 = 0.5*h*(k1 - k2)/(c1 - c2);

xc = NaN(2, length(t));
xdotc = NaN(2, length(t));
xnode = NaN(2, length(tab.c));
for i = 1:2
    xc(i,:) = x0(i) - rho1(i)*(t/h) + rho2(i)*(t/h).^2;
    xdotc(i,:) = -rho1(i)/h + 2*rho2(i)/h*(t/h);
    xnode(i,:) = x0(i) - rho1(i)*tab.c + rho2(i)*tab.c.^2;
end

% Get tangent lines for derivative.
dt = 0.03; % Width for tangents.
tt = [tab.c(1) + [-dt, dt], NaN(), tab.c(2) + [-dt, dt]];
xt = [xnode(:,1) - h*k1*dt, xnode(:,1) + h*k1*dt, NaN(2, 1), ...
      xnode(:,2) - h*k2*dt, xnode(:,2) + h*k2*dt];

% Make a plot.
figure();
t = t/h;
for i = 1:2
    subplot(2, 2, i);
    hold('on');
    plot(t, x(i,:), '-k', t, xc(i,:), '--k');
    plot(tt, xt(i,:), '-k', 'linewidth', 5);
    ylabel(sprintf('x_%d', i), 'rotation', 0);
    xlabel('\tau/h');

    subplot(2, 2, 2 + i);
    hold('on');
    plot(t, xdot(i,:), '-k', t, xdotc(i,:), '--k');
    plot(tab.c, [k1(i), k2(i)], 'ok', 'markerfacecolor', 'k');
    ylabel(sprintf('dx_%d/dt', i), 'rotation', 0);
    xlabel('\tau/h');
end

% Save data.
data = struct('t', t, 'x', x, 'xdot', xdot, 'xc', xc, 'xdotc', xdotc, ...
              'tt', tt, 'xt', xt, 'h', h, 'k', k, 'c', tab.c, 'xn', xnode);
save('-v7', 'gl4illustration.mat', '-struct', 'data');

function [xdot, jac] = stiffode(x)
% [xdot, jac] = stiffode(x)
%
% Example stiff ODE. First output is dx/dt; second is Jacobian of dx/dt.
A = [0 1; -1 0];
lambda = 500;
rho = (x'*x - 1);
xdot = A*x - lambda*rho*x;
jac = A - eye(2)*lambda*rho - 2*lambda*(x*x');
end%function

function [xplus, converged, iter, k] = butcherintegration(x, func, h, a, b, Niter)
% [xplus, converged, iter, k] = butcherintegration(x, f, h, a, b, [Niter=100])
%
% Performs numerical integration of f starting from x for h time units.
%
% f must be a function handle that returns dx/dt. If a defines an implicit
% method, then f must also return the Jacobian of dx/dt. In either case, f must
% be time-invariant.
%
% a and b should be the coefficients from the Butcher tableau for the
% integration scheme. If the a matrix defines an explicit scheme (i.e., zero on
% and above the diagonal, then the explicit formulas will be used. Otherwise,
% the implicit equations are solved using Newton's method.
%
% Niter is the maximum number of iterations to perform. xguess is a guess for
% x at the next timestep.
narginchk(5, 6);
if nargin() < 6
    Niter = 100;
end

% Check arguments.
x = x(:);
Nx = length(x);
b = h*b(:);
Ns = length(b);
if ~isequal(size(a), [Ns, Ns])
    error('a must be a square matrix with a column for each element of b!');
end
a = h*a;

% Decide if tableau is explicit or not.
if nnz(triu(a)) == 0
    k = butcher_explicit(x, func, a);
    converged = true();
    iter = 0;
else
    [k, converged, iter] = butcher_implicit(x, func, a, Niter);
end

% Advance state.
xplus = x + k*b;

end%function

% *****************************************************************************

function k = butcher_explicit(x, func, a)
% Performs an explicit integration step.
Nx = length(x);
Ns = size(a, 1);
k = zeros(Nx, Ns);
for i = 1:Ns
    k(:,i) = func(x + k*a(i,:)');
end
end%function

% *****************************************************************************

function [k, converged, iter] = butcher_implicit(x, func, a, Niter)
% Solves implicit integration using Newton's Method. Timestep is assumed 1.
narginchk(4, 4);
Nx = length(x);
Ns = size(a, 1);;
k = zeros(Nx, Ns);
converged = false();
iter = 0;
while ~converged && iter < Niter
    iter = iter + 1;

    dX = k*a';
    if Ns == 1
        [f, J] = func(x + dX);
        J = a*J;
    else
        f = zeros(Nx*Ns, 1);
        J = zeros(Nx*Ns, Nx*Ns);
        for i = 1:Ns
            I = ((i - 1)*Nx + 1):(i*Nx);
            [f(I), Ji] = func(x + dX(:,i));
            J(I,:) = kron(a(i,:), Ji);
        end
    end

    f = f - k(:);
    J = J - eye(Nx*Ns);
    dk = -J\f;
    k = k + reshape(dk, [Nx, Ns]);

    converged = norm(dk) < 1e-8;
end

end%function

Figure 8.3:

Polynomial approximation \tilde {x}_1(t) and true trajectory x_1(t) of the first state and its derivative.

Code for Figure 8.3

main.m

stiffode.m

butcherintegration.m