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

        function g = testEquationGain(obj, z)
            default_arg('z', sym('z'));
            s = obj.nStages();

            b = sym(obj.b);
            A = sym(obj.a);
            one = sym(ones(s,1));
            I = sym(eye(s));

            g = abs(1 + z*b*inv(I-z*A)*one);
        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