Figure 8.4(b)

Code for Figure 8.4(b)

main.m

% Example numerical integration of a stiff ODE.

% Butcher tableaus (only a and b) for the three methods.
tabs = struct();
tabs.ieu = struct('a', 1, 'b', 1);
tabs.mid = struct('a', 0.5, 'b', 1);
tabs.gl4 = struct('a', [1/4, 1/4 - sqrt(3)/6; 1/4 + sqrt(3)/6, 1/4], ...
                  'b', [1/2, 1/2]);
tabnames = fieldnames(tabs);

% Run each method.
f = @stiffode;
x0 = [1; 0];
Neval = 10;
T = 2*pi();
x = struct();
for n = 1:length(tabnames)
    n = tabnames{n};
    Nstep = Neval/length(tabs.(n).b); % Use Neval function evaluations.
    h = T/Nstep;
    x.(n) = zeros(2, Nstep + 1);
    x.(n)(:,1) = x0;
    for i = 1:Nstep
        x.(n)(:,i + 1) = butcherintegration(x.(n)(:,i), f, h, ...
                                            tabs.(n).a, tabs.(n).b);
    end
end

% Add exact solution.
theta = linspace(0, 2*pi(), 257);
x.exact = [cos(theta); -sin(theta)];
tabnames = [{'exact'}; tabnames];

% Make a plot.
errorplot2d(x, T, 'SouthEast');

stiffode.m

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

butcherintegration.m

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

errorplot2d.m

function errorplot2d(x, T, location)
% errorplot2d(x, T, [location])
%
% Makes a plot of integration error over the given time interval T.
%
% x must be a struct with fields to plot. All step sizes are assumed to be
% uniform. If x contains an 'exact' field, it is handled specially.
narginchk(2, 3);
if nargin() < 3
    location = {};
else
    location = {'location', location};
end

% Create figure and axes.
figure();

ax1 = subplot(2, 2, 1);
hold(ax1, 'on');
ylabel(ax1, 'x_1', 'rotation', 0);
xlabel(ax1, 't');

ax2 = subplot(2, 2, 3);
hold(ax2, 'on');
ylabel(ax2, 'x_2', 'rotation', 0);
xlabel(ax2, 't');

axp = subplot(2, 2, [2, 4]);
hold(axp, 'on');
ylabel(axp, 'x_2', 'rotation', 0);
xlabel(axp, 'x_1');

% Plot exact solution.
if isfield(x, 'exact')
    t = linspace(0, T, size(x.exact, 2));
    plot(ax1, t, x.exact(1,:), '-k');
    plot(ax2, t, x.exact(2,:), '-k');
    plot(axp, x.exact(1,:), x.exact(2,:), '-k');
    exact = {'exact'};
else
    exact = {};
end

% Plot other solutions.
fields = setdiff(fieldnames(x), {'exact'});
markers = {'o', 's', 'd', '^', 'v', '>', '<', 'p', 'h'};
markers = markers(mod(0:(length(fields) - 1), length(markers)) + 1);
colors = hsv(length(fields));
for i = 1:length(fields)
    f = fields{i};
    t = linspace(0, T, size(x.(f), 2));
    style = {markers{i}, 'color', colors(i,:)};
    plot(ax1, t, x.(f)(1,:), style{:});
    plot(ax2, t, x.(f)(2,:), style{:});
    plot(axp, x.(f)(1,:), x.(f)(2,:), style{:});
end

% Add legend.
legend(axp, exact{:}, fields{:}, location{:});

end%function

Figure 8.4(b):

Performance of implicit integration methods on a stiff ODE. Simulation result for M=10 points.

Code for Figure 8.4(b)

main.m

stiffode.m

butcherintegration.m

errorplot2d.m