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