Mercurial > repos > public > sbplib
changeset 990:1066bb31bc95 feature/timesteppers
Create class for butcher tableau
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Wed, 09 Jan 2019 09:09:15 +0100 |
| parents | e41c93d7ab08 |
| children | a99f00896b8e |
| files | +time/+rk/ButcherTableau.m +time/+rk/butcherTableauFromStr.m |
| diffstat | 2 files changed, 127 insertions(+), 42 deletions(-) [+] |
line wrap: on
line diff
diff -r e41c93d7ab08 -r 1066bb31bc95 +time/+rk/ButcherTableau.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rk/ButcherTableau.m Wed Jan 09 09:09:15 2019 +0100 @@ -0,0 +1,115 @@ +classdef ButcherTableau + properties + a,b,c + end + + methods + % A ButcherTableau describes a specific rungekutta method where + % y(n+1) = y(n) + dt*b(i)*k(i) + % k(i) = F(t + c(i)*dt, Y(i)) + % Y(i) = y(i) + dt*a(i,j)*k(j) + % where repeating indecies imply summation + function obj = ButcherTableau(a,b,c) + s = length(c); + assertSize(a, [s,s]); + assertLength(b, s); + + obj.a = a; + obj.b = b; + obj.c = c; + end + + function s = nStages(obj) + s = length(obj.c); + end + + function b = isExplicit(obj) + b = all(all(triu(obj.a)==0)); + end + + % TBD: Add functions for checking accuracy, stability? + end + + methods(Static) + % TVD (Total Variational Diminishing) + function bt = tvd_3() + a = [ + 0, 0, 0; + 1, 0, 0; + 1/4, 1/4, 0; + ]; + b = [1/6, 1/6, 2/3]; + c = [0 1 1/2]; + + bt = time.rk.ButcherTableau(a,b,c); + end + + % Standard RK4 + function bt = rk4() + a = [ + 0, 0, 0, 0; + 1/2, 0, 0, 0; + 0, 1/2, 0, 0; + 0, 0, 1, 0; + ]; + + b = [1/6 1/3 1/3 1/6]; + c = [0, 1/2, 1/2, 1]; + + bt = time.rk.ButcherTableau(a,b,c); + end + + % 3/8 RK4 (Kuttas method). Lower truncation error, more flops. + % Irreducible, unlike standard rk4. + function bt = rk4_3_8() + a = [ + 0, 0, 0, 0; + 1/3, 0, 0, 0; + -1/3, 1, 0, 0; + 1, -1, 1, 0; + ]; + + b = [1/8 3/8 3/8 1/8]; + c = [0, 1/3, 2/3, 1]; + + bt = time.rk.ButcherTableau(a,b,c); + end + + % Runge-Kutta 6 from Alshina07 + function bt = rk6() + a = zeros(7,7); + + a(2,1) = 4/7; + + a(3,1) = 115/112; + a(3,2) = -5/16; + + a(4,1) = 589/630; + a(4,2) = 5/18; + a(4,3) = -16/45; + + a(5,1) = 229/1200 - 29/6000*sqrt(5); + a(5,2) = 119/240 - 187/1200*sqrt(5); + a(5,3) = -14/75 + 34/375*sqrt(5); + a(5,4) = -3/100*sqrt(5); + + a(6,1) = 71/2400 - 587/12000*sqrt(5); + a(6,2) = 187/480 - 391/2400*sqrt(5); + a(6,3) = -38/75 + 26/375*sqrt(5); + a(6,4) = 27/80 - 3/400*sqrt(5); + a(6,5) = (1+sqrt(5))/4 + + a(7,1) = -49/480 + 43/160*sqrt(5); + a(7,2) = -425/96 + 51/32*sqrt(5); + a(7,3) = 52/15 - 4/5*sqrt(5); + a(7,4) = -27/16 + 3/16*sqrt(5); + a(7,5) = 5/4 - 3/4*sqrt(5); + a(7,6) = 5/2 - 1/2*sqrt(5); + + b = [1/12 0 0 0 5/12 5/12 1/12]; + c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1]; + + bt = time.rk.ButcherTableau(a,b,c); + end + end +end
diff -r e41c93d7ab08 -r 1066bb31bc95 +time/+rk/butcherTableauFromStr.m --- a/+time/+rk/butcherTableauFromStr.m Wed Jan 09 08:56:42 2019 +0100 +++ b/+time/+rk/butcherTableauFromStr.m Wed Jan 09 09:09:15 2019 +0100 @@ -6,46 +6,16 @@ % @return b - weights for summing stages % @return c - time step coefficents for intermediate stages function [s,a,b,c] = butcherTableauFromStr(method) -switch method - case "tvd-3" - % TVD (Total Variational Diminishing) - s = 3; - a = zeros(s,s-1); - a(2,1) = 1; - a(3,1) = 1/4; a(3,2) = 1/4; - b = [1/6, 1/6, 2/3]; - c = [0 1 1/2]; - case "rk4" - % Standard RK4 - s = 4; - a = zeros(s,s-1); - a(2,1) = 1/2; - a(3,1) = 0; a(3,2) = 1/2; - a(4,1) = 0; a(4,2) = 0; a(4,3) = 1; - b = [1/6 1/3 1/3 1/6]; - c = [0, 1/2, 1/2, 1]; - case "rk4-3/8" - % 3/8 RK4 (Kuttas method). Lower truncation error, more flops. - % Irreducible, unlike standard rk4. - s = 4; - a = zeros(s,s-1); - a(2,1) = 1/3; - a(3,1) = -1/3; a(3,2) = 1; - a(4,1) = 1; a(4,2) = -1; a(4,3) = 1; - b = [1/8 3/8 3/8 1/8]; - c = [0, 1/3, 2/3, 1]; - case "rk6" - % Runge-Kutta 6 from Alshina07 - s = 7; - a = zeros(s,s-1); - a(2,1) = 4/7; - a(3,1) = 115/112; a(3,2) = -5/16; - a(4,1) = 589/630; a(4,2) = 5/18; a(4,3) = -16/45; - a(5,1) = 229/1200 - 29/6000*sqrt(5); a(5,2) = 119/240 - 187/1200*sqrt(5); a(5,3) = -14/75 + 34/375*sqrt(5); a(5,4) = -3/100*sqrt(5); - a(6,1) = 71/2400 - 587/12000*sqrt(5); a(6,2) = 187/480 - 391/2400*sqrt(5); a(6,3) = -38/75 + 26/375*sqrt(5); a(6,4) = 27/80 - 3/400*sqrt(5); a(6,5) = (1+sqrt(5))/4; - a(7,1) = -49/480 + 43/160*sqrt(5); a(7,2) = -425/96 + 51/32*sqrt(5); a(7,3) = 52/15 - 4/5*sqrt(5); a(7,4) = -27/16 + 3/16*sqrt(5); a(7,5) = 5/4 - 3/4*sqrt(5); a(7,6) = 5/2 - 1/2*sqrt(5); - b = [1/12 0 0 0 5/12 5/12 1/12]; - c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1]; - otherwise + try + bt = time.rk.ButcherTableau.(method); + catch error('Runge-Kutta method %s is not implemented', method) -end \ No newline at end of file + end + + assert(bt.isExplicit()); + + s = bt.nStages(); + a = bt.a(:,1:end-1); + b = bt.b; + c = bt.c; +end