function temple5044_rk_order_conditions(A,b,c)
%TEMPLE5044_RK_ORDER_CONDITIONS
%   Verifies how well order conditions (up to order 5)
%   are satisfied by the Runge-Kutta scheme provided
%   via the Butcher tableau (A,b,c).
%   If only A and b are given, c=A*e is assumed, where
%   e=(1,...,1)^T.
%
% (C) 2019/02/06 by Benjamin Seibold
%                   http://www.math.temple.edu/~seibold/

r = size(A,2); % number of stages
e = ones(r,1); % vector of all ones
if nargin<3, c = A*e; end % choose c according to solving t'=1 exactly
b = b(:); % make vectors b and c column vectors
if r~=size(A,1)||r~=length(b)||r~=length(c) % check input dimensions
    error('Input has inconsistent dimensions.')
end

% Calculate error in order conditions (up to order 5)
tau1 = abs(b'*e-1);
tau2 = [abs(b'*c-1/2),norm(c-A*e,1)];
tau3 = abs([b'*A*c-1/6,b'*c.^2-1/3]);
tau4 = abs([b'*A^2*c-1/24,b'*A*c.^2-1/12,b'*((A*c).*c)-1/8,b'*c.^3-1/4]);
tau5 = abs([b'*A^3*c-1/120,...
    b'*A*c.^3-1/20,...
    b'*A^2*c.^2-1/60,...
    b'*A*((A*c).*c)-1/14,...
    b'*((A*c.^2).*c)-1/15,...
    b'*((A^2*c).*c)-1/30,...
    b'*((A*c).*c.^2)-1/10,...
    b'*(A*c).^2-1/20,...
    b'*c.^4-1/5]);

% Evaluate whether conditions are satisfied
tol = 1e-12; % tolerance to call condition satisfied
word = {'no ','yes'}; % verbal formulation of boolean
order_error = [norm(tau1,1),norm(tau2,1),...
    norm(tau3,1),norm(tau4,1),norm(tau5,1)];
order_error_linear = [tau1(1),tau2(1),tau3(1),tau4(1),tau5(1)];
order_word = word((order_error<tol)+1);
order_word_linear = word((order_error_linear<tol)+1);

% Display Butcher tableau
number_format = '%10.4g'; % format of numbers displayed
fprintf('Runge-Kutta method (Butcher tableau):\n')
fmt = [number_format,' |',repmat([number_format],1,r-1),number_format,'\n'];
fprintf(fmt,[c,A]');
for j = 0:r, fprintf('----------'), end
fprintf('---\n')
fmt = ['           |',repmat(number_format,1,r-1),number_format,'\n'];
fprintf(fmt,b');

% Display which order conditions are satisfied
fprintf('\nConsistency (order 1) conditions: %s [%0.0e]',...
    order_word{1},order_error(1))
fprintf('   (linear: %s [%0.0e])\n',...
    order_word_linear{1},order_error_linear(1))
for j = 2:5
    fprintf('Order %d conditions              : %s [%0.0e]',...
        j,order_word{j},order_error(j))
    fprintf('   (linear: %s [%0.0e])\n',...
        order_word_linear{j},order_error_linear(j))
end
fprintf('\n')
