Code covered by the BSD License  

Highlights from
First-order stiff ordinary differential equation solver

from First-order stiff ordinary differential equation solver by Ali Dinler
Runs 20 implicit and semi-implicit methods for first-order initial value stiff ODEs.

[RKa,RKb,RKc,RKns,RKord,RKtype,RKname]=implicit_setmethod(name)
 function [RKa,RKb,RKc,RKns,RKord,RKtype,RKname]=implicit_setmethod(name)
% SETMETHOD - Returns Butcher tableau coefficients for Runge-Kutta methods
%
% CALL [RKa,RKb,RKc,RKns,RKord,RKtype,RKname]=implicit_setmethod(name);
%
% INPUT:
%        name   - abbreviation for the method to be used
%  
%                 'GL2' 'IM2'        - Gauss-Legendre/implicit midpoint, order 2
%                 'GL4'              - Gauss-Legendre, order 4
%                 'GL6'              - Gauss-Legendre, order 6
%                 'radIA1'           - Radau IA method of order 1
%				  'radIA3'	         - Radau IA method of order 3
%				  'radIA5'	         - Radau IA method of order 5
%                 'radIIA1' 'BE'     - Radau IIA method of order 1
%				  'radIIA3'	         - Radau IIA method of order 3
%				  'radIIA5'	         - Radau IIA method of order 5
%				  'lobIIIA2' 'TrapM' - Lobatto IIIA method of order 2/Trapezoidal method
%				  'lobIIIA4'	     - Lobatto IIIA method of order 4
%				  'lobIIIA6'	     - Lobatto IIIA method of order 6
%				  'lobIIIB2'	     - Lobatto IIIB method of order 2
%				  'lobIIIB4'	     - Lobatto IIIB method of order 4
%				  'lobIIIB6'	     - Lobatto IIIB method of order 6
%				  'lobIIIC2'	     - Lobatto IIIC method of order 2
%				  'lobIIIC4'	     - Lobatto IIIC method of order 4
%				  'lobIIIC6'	     - Lobatto IIIC method of order 6
%				  'SDIRK3'			 - Singly diagolally implicit RK method of order 3
%				  'SDIRK4'			 - Singly diagolally implicit RK method of order 4					
%
% OUTPUT:
%        RKa    - matrix A from the method's Butcher tableau
%        RKb    - vector b from the method's Butcher tableau
%        RKc    - vector c from the method's Butcher tableau or the
%                 abscissae values in case of Magnus methods
%        RKns   - number of stages in the method
%        RKord  - order of the method
%        RKtype - type of the method ('explicit', 'implicit', 'SDIRK')
%        RKname - name of the method (the same as input "name")
% 
%	Written by Ali Dinler 2006



if (strcmp(name,'GL2')|strcmp(name,'IM2')),
  % Gauss-Legendre/implicit midpoint of order 2
  RKc = [1/2];
  RKa = [1/2];
  RKb = [1];
  RKord = 2;
  RKns = 1;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'GL4'),
  % Gauss-Legendre method of order 4
  RKc = [(3-sqrt(3))/6 (3+sqrt(3))/6];
  RKa = [       1/4       (3-2*sqrt(3))/12;
         (3+2*sqrt(3))/12        1/4  ];
  RKb = [1/2 1/2];
  RKord = 4;
  RKns = 2;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'GL6'),
  % Gauss-Legendre method of order 6
  RKc = [(5-sqrt(15))/10 1/2 (5+sqrt(15))/10];
  RKa = [       5/36         (10-3*sqrt(15))/45 (25-6*sqrt(15))/180;
         (10+3*sqrt(15))/72       2/9           (10-3*sqrt(15))/72;
         (25+6*sqrt(15))/180 (10+3*sqrt(15))/45       5/36  ];
  RKb = [5/18 4/9 5/18];
  RKord = 6;
  RKns = 3;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'radIA1'),
  % Radau IA method of order 1
  RKc = [0];
  RKa = [1];
  RKb = [1];
  RKord = 1;
  RKns = 1;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'radIA3'),
  % Radau IA method of order 3
  RKc = [0,2/3];
  RKa = [1/4,-1/4;1/4,5/12];
  RKb = [1/4,3/4];
  RKord = 3;
  RKns = 2;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'radIA5'),
  % Radau IA method of order 5
  RKc = [0,(6-sqrt(6)/10),(6+sqrt(6)/10)];
  RKa = [1/9,(-1-sqrt(6))/18,(-1+sqrt(6))/18;
         1/9,(88+7*sqrt(6))/360,(88-43*sqrt(6))/360;
         1/9,(88+43*sqrt(6))/360,(88-7*sqrt(6))/360];
  RKb = [1/9,(16+sqrt(6))/36,(16-sqrt(6))/36];
  RKord = 5;
  RKns = 3;
  RKtype = 'implicit';
  RKname = name;
  RKname = name;
elseif (strcmp(name,'radIIA1')|strcmp(name,'BE')),
  % Radau IIA method of order 1/Bakward Euler
  RKc = [1];
  RKa = [1];
  RKb = [1];
  RKord = 1;
  RKns = 1;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'radIIA3'),
  % Radau IIA method of order 3
  RKc = [1/3,1];
  RKa = [5/12,-1/12;3/4,1/4];
  RKb = [3/4,1/4];
  RKord = 3;
  RKns = 2;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'radIIA5'),
  % Radau IIA method of order 5
  RKc = [(4-sqrt(6))/10,(4+sqrt(6))/10,1];
  RKa = [(88-7*sqrt(6))/360,(296-169*sqrt(6))/1800,(-2+3*sqrt(6))/225;
         (296+169*sqrt(6))/1800,(88+7*sqrt(6))/360,(-2-3*sqrt(6))/225;
         (16-sqrt(6))/36,(16+sqrt(6))/36,1/9];
  RKb = [(16-sqrt(6))/36,(16+sqrt(6))/36,1/9];
  RKbhat = '?';
  RKord = 5;
  RKns = 3;
  RKtype = 'implicit';
  RKname = name;
elseif (strcmp(name,'lobIIIA2')|strcmp(name,'TrapM')),
  % Lobatto IIIA method of order 2/Trapezoidal method
  RKc = [0,1];
  RKa = [0,0;1/2,1/2];
  RKb = [1/2,1/2];
  RKord = 2;
  RKns = 2;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'lobIIIA4'),
  % Lobatto IIIA method of order 4
  RKc = [0,1/2,1];
  RKa = [0,0,0;
         5/24,1/3,-1/24;
         1/6,2/3,1/6];
  RKb = [1/6,2/3,1/6];
  RKord = 4;
  RKns = 3;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'lobIIIA6'),
  % Lobatto IIIA method of order 6
  RKc = [0,(5-sqrt(5))/10,(5+sqrt(5))/10,1];
  RKa = [0,0,0,0;(11+sqrt(5))/120,(25-sqrt(5))/120,(25-13*sqrt(5))/120,(-1+sqrt(5))/120;
        (11-sqrt(5))/120,(25+13*sqrt(5))/120,(25+sqrt(5))/120,(-1-sqrt(5))/120;
        1/12,5/12,5/12,1/12];
  RKb = [1/12,5/12,5/12,1/12];
  RKord = 6;
  RKns = 4;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'lobIIIB2'),
  % Lobatto IIIB method of order 2
  RKc = [0,1];
  RKa = [1/2,0;1/2,0];
  RKb = [1/2,1/2];
  RKord = 2;
  RKns = 2;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'lobIIIB4'),
  % Lobatto IIIB method of order 4
  RKc = [0,1/2,1];
  RKa = [1/6,-1/6,0;1/6,1/3,0;1/6,5/6,0];
  RKb = [1/6,2/3,1/6];
  RKord = 4;
  RKns = 3;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'lobIIIB6'),
  % Lobatto IIIB method of order 6
  RKc = [0,(5-sqrt(5))/10,(5+sqrt(5))/10,1];
  RKa = [1/12,(-1-sqrt(5))/24,(-1+sqrt(5))/24,0;
         1/12,(25+sqrt(5))/120,(25-13*sqrt(5))/120,0;
         1/12,(25+13*sqrt(5))/120,(25-sqrt(5))/120,0;
         1/12,(11-sqrt(5))/24,(11+sqrt(5))/24,0];
  RKb = [1/12,5/12,5/12,1/12];
  RKord = 6;
  RKns = 4;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'lobIIIC2'),
  % Lobatto IIIC method of order 2
  RKc = [0,1];
  RKa = [1/2,-1/2;1/2,1/2];
  RKb = [1/2,1/2];
  RKord = 2;
  RKns = 2;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'lobIIIC4'),
  % Lobatto IIIC method of order 4
  RKc = [0,1/2,1];
  RKa = [1/6,-1/3,1/6;1/6,5/12,-1/12;1/6,2/3,1/6];
  RKb = [1/6,2/3,1/6];
  RKord = 4;
  RKns = 3;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'lobIIIC6'),
  % Lobatto IIIC method of order 6
  RKc = [0,(5-sqrt(5))/10,(5+sqrt(5))/10,1];
  RKa = [1/12,-sqrt(5)/12,sqrt(5)/12,-1/12;
         1/12,1/4,(10-7*sqrt(5))/60,sqrt(5)/60;
         1/12,(10+7*sqrt(5))/60,1/4,-sqrt(5)/60;
         1/12,5/12,5/12,1/12];
  RKb = [1/12,5/12,5/12,1/12];
  RKord = 6;
  RKns = 4;
  RKtype = 'implicit';
  RKname = name;
elseif strcmp(name,'SDIRK3'),
  % Singly Diagonally Implicit RK method of order 3
  gamma=0.4358665215;
  RKc = [gamma,0.7179332608,1];
  RKa = [gamma,0,0;
         0.2820667392,gamma,0;
         1.208496649,-0.644363171,gamma];
  RKb = [1.208496649,-0.644363171,gamma];
  RKord = 3;
  RKns = 3;
  RKtype = 'diagonally implicit';
  RKname = name;
elseif strcmp(name,'SDIRK4'),
  % Singly Diagonally Implicit RK method of order 4
  gamma=1/4;
  RKc = [gamma,3/4,11/20,1/2,1];
  RKa = [gamma,0,0,0,0;
         1/2,gamma,0,0,0;
         17/50,-1/25,gamma,0,0;
         371/1360,-137/2720,15/544,gamma,0;
         25/24,-49/48,125/16,-85/12,gamma];
  RKb = [25/24,-49/48,125/16,-85/12,gamma];
  RKord = 4;
  RKns = 5;
  RKtype = 'diagonally implicit';
  RKname = name;
else,
  error(['**The method ' name ' is unknown.']);
end;
return;

Contact us at files@mathworks.com