| [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;
|
|