Mercurial > repos > public > sbplib

changeset 990:1066bb31bc95 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
--- /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
--- 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