Mercurial > repos > public > sbplib

changeset 1112:835c8fa456ec feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 10 Apr 2019 22:22:46 +0200
parents 1a265a376b36 (diff) 3230e4cbdbb4 (current diff)
children 47e86b5270ad
files
diffstat 14 files changed, 475 insertions(+), 354 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/ButcherTableau.m	Wed Apr 10 22:22:46 2019 +0200
@@ -0,0 +1,127 @@
+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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/Explicit.m	Wed Apr 10 22:22:46 2019 +0200
@@ -0,0 +1,139 @@
+classdef Explicit < time.Timestepper
+    properties
+        F       % RHS of the ODE
+        dt      % Time step
+        t       % Time point
+        v       % Solution vector
+        n       % Time level
+        bt
+    end
+
+
+    methods
+        % Timesteps v_t = F(t,v), using the specified ButcherTableau
+        % from t = t0 with timestep dt and initial conditions v(0) = v0
+        function obj = Explicit(F, dt, t0, v0, bt)
+            assertType(bt, 'time.rk.ButcherTableau')
+            obj.F = F;
+            obj.dt = dt;
+            obj.t = t0;
+            obj.v = v0;
+            obj.n = 0;
+
+            assert(bt.isExplicit())
+            obj.bt = bt;
+        end
+
+        % v: Current solution
+        % t: Current time
+        function [v,t] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            s = obj.bt.nStages();
+            a = obj.bt.a;
+            b = obj.bt.b;
+            c = obj.bt.c;
+
+            % Compute rates K
+            K = zeros(length(v), s);
+            for i = 1:s
+                V_i = obj.v;
+                for j = 1:i-1
+                    V_i = V_i + dt*a(i,j)*K(:,j);
+                end
+                K(:,i) = F(t+dt*c(i), V_i);
+            end
+
+            % Compute updated solution
+            v_next = v;
+            for i = 1:s
+                v_next = v_next + dt*b(i)*K(:,i);
+            end
+
+            obj.v = v_next;
+            obj.t = obj.t + obj.dt;
+            obj.n = obj.n + 1;
+        end
+
+        % TBD: Method name
+        % TBD: Parameter name
+        %
+        % Takes a regular step but with discreteRates(:,i) added to RHS for stage i.
+        %  v_t = F(t,v) + discreteRates(:, ...)
+        %
+        % Also returns the stage approximations (V) and stage rates (K).
+        function [v,t, V, K] = stepWithDiscreteData(obj, discreteRates)
+            s = obj.bt.nStages();
+            a = obj.bt.a;
+            b = obj.bt.b;
+            c = obj.bt.c;
+
+            % Compute rates K and stage approximations V
+            K = zeros(length(v), s);
+            V = zeros(length(v), s);
+            for i = 1:s
+                V_i = obj.v;
+                for j = 1:i-1
+                    V_i = V_i + dt*a(i,j)*K(:,j);
+                end
+
+                K_i = F(t+dt*c(i), V_i);
+                K_i = K_i + discreteRates(:,i);
+
+                V(:,i) = V_i;
+                K(:,i) = K_i;
+            end
+
+            % Compute updated updated solution
+            v_next = v;
+            for i = 1:s
+                v_next = v_next + dt*b(i)*K(:,i);
+            end
+
+            obj.v = v_next;
+            obj.t = obj.t + obj.dt;
+            obj.n = obj.n + 1;
+        end
+
+        % Returns a vector of time points, including substage points,
+        % in the time interval [t0, tEnd].
+        % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
+        function tvec = timePoints(obj, t0, tEnd)
+            % TBD: Should this be implemented here or somewhere else?
+            N = round( (tEnd-t0)/obj.dt );
+            tvec = zeros(N*obj.s, 1);
+            s = obj.coeffs.s;
+            c = obj.coeffs.c;
+            for i = 1:N
+                ind = (i-1)*s+1 : i*s;
+                tvec(ind) = ((i-1) + c')*obj.dt;
+            end
+        end
+
+        % Returns a vector of quadrature weights corresponding to grid points
+        % in time interval [t0, tEnd], substage points included.
+        % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
+        function weights = quadWeights(obj, t0, tEnd)
+            % TBD: Should this be implemented here or somewhere else?
+            N = round( (tEnd-t0)/obj.dt );
+            b = obj.coeffs.b;
+            weights = repmat(b', N, 1);
+        end
+    end
+
+    methods(Static)
+        % TBD: Function name
+        function ts = methodFromStr(F, dt, t0, v0, methodStr)
+            try
+                bt = time.rk.ButcherTableau.(method);
+            catch
+                error('Runge-Kutta method ''%s'' is not implemented', methodStr)
+            end
+
+            ts = time.rk.Explicit(F, dt, t0, v0, bt);
+        end
+    end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/ExplicitSecondOrder.m	Wed Apr 10 22:22:46 2019 +0200
@@ -0,0 +1,119 @@
+classdef ExplicitSecondOrder < time.Timestepper
+    properties
+        F       % RHS of the ODE
+        dt      % Time step
+        t       % Time point
+        v, vt    % Solution state
+        n       % Time level
+        bt
+    end
+
+
+    methods
+        % Timesteps v_tt = F(t,v,vt), using the specified ButcherTableau
+        % from t = t0 with timestep dt and initial conditions v(0) = v0
+        function obj = ExplicitSecondOrder(F, dt, t0, v0, v0t, bt)
+            assertType(bt, 'time.rk.ButcherTableau')
+            obj.F = F;
+            obj.dt = dt;
+            obj.t = t0;
+            obj.v = v0;
+            obj.vt = v0t;
+            obj.n = 0;
+
+            assert(bt.isExplicit())
+            obj.bt = bt;
+        end
+
+        function [v,t] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+        end
+
+        function [vt,t] = getVt(obj)
+            vt = obj.vt;
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            s = obj.bt.nStages();
+            a = obj.bt.a;
+            b = obj.bt.b;
+            c = obj.bt.c;
+
+            t = obj.t;
+            v = obj.v;
+            vt = obj.vt;
+            dt = obj.dt;
+
+            k1 = obj.F(t,          v,                            v_t);
+            k2 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t,               v_t + 1/2*dt*k1);
+            k3 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t + 1/4*dt^2*k1, v_t + 1/2*dt*k2);
+            k4 = obj.F(t + dt,     v + dt*v_t + 1/2*dt^2*k2,     v_t + dt*k3);
+
+            % Compute rates K
+            K = zeros(length(v), s);
+            for i = 1:s
+            	U_i = obj.v;
+                V_i = obj.vt;
+                for j = 1:i-1
+                    U_i = U_i % + dt*a(i,j)*K(:,j);
+                    V_i = V_i % + dt*a(i,j)*K(:,j);
+                end
+                K(:,i) = F(t+dt*c(i), U_i, V_i);
+            end
+
+            % Compute updated solution
+            v_next = v;
+            vt_next = vt;
+            for i = 1:s
+                v_next  = v_next % + dt*b(i)*K(:,i);
+                vt_next = vt_next % + dt*b(i)*K(:,i);
+            end
+
+            obj.v  = v_next;
+            obj.vt = vt_next;
+            obj.t = obj.t + obj.dt;
+            obj.n = obj.n + 1;
+        end
+
+
+        % Returns a vector of time points, including substage points,
+        % in the time interval [t0, tEnd].
+        % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
+        function tvec = timePoints(obj, t0, tEnd)
+            % TBD: Should this be implemented here or somewhere else?
+            N = round( (tEnd-t0)/obj.dt );
+            tvec = zeros(N*obj.s, 1);
+            s = obj.coeffs.s;
+            c = obj.coeffs.c;
+            for i = 1:N
+                ind = (i-1)*s+1 : i*s;
+                tvec(ind) = ((i-1) + c')*obj.dt;
+            end
+        end
+
+        % Returns a vector of quadrature weights corresponding to grid points
+        % in time interval [t0, tEnd], substage points included.
+        % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
+        function weights = quadWeights(obj, t0, tEnd)
+            % TBD: Should this be implemented here or somewhere else?
+            N = round( (tEnd-t0)/obj.dt );
+            b = obj.coeffs.b;
+            weights = repmat(b', N, 1);
+        end
+    end
+
+    methods(Static)
+        % TBD: Function name
+        function ts = methodFromStr(F, dt, t0, v0, methodStr)
+            try
+                bt = time.rk.ButcherTableau.(method);
+            catch
+                error('Runge-Kutta method ''%s'' is not implemented', methodStr)
+            end
+
+            ts = time.rk.Explicit(F, dt, t0, v0, bt);
+        end
+    end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/get_rk4_time_step.m	Wed Apr 10 22:22:46 2019 +0200
@@ -0,0 +1,18 @@
+% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values.
+function k = get_rk4_time_step(lambda,l_type)
+    default_arg('l_type','complex')
+	assertIsMember(l_type, {'real', 'imag', 'complex'})
+
+    rad = abs(lambda);
+    switch l_type
+    	case 'real'
+	        % kl > -2.7852
+	        k = 2.7852/rad;
+    	case 'imag'
+		    % |kl| < 2.8284
+		    k = 2.8284/rad;
+		case 'complex'
+	        % |kl| < 2.5
+	        k = 2.5/rad;
+    end
+end
--- a/+time/+rk4/get_rk4_time_step.m	Wed Apr 10 22:11:11 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,21 +0,0 @@
-% Calculate the size of the largest time step given the largest evalue for a operator with pure imaginary e.values.
-function k = get_rk4_time_step(lambda,l_type)
-    default_arg('l_type','complex')
-
-    rad = abs(lambda);
-    if strcmp(l_type,'real')
-        % Real eigenvalue
-        % kl > -2.7852
-        k = 2.7852/rad;
-
-    elseif strcmp(l_type,'imag')
-        % Imaginary eigenvalue
-        % |kl| < 2.8284
-        k = 2.8284/rad;
-    elseif strcmp(l_type,'complex')
-        % |kl| < 2.5
-        k = 2.5/rad;
-    else
-        error('l_type must be one of ''real'',''imag'' or ''complex''.')
-    end
-end
\ No newline at end of file
--- a/+time/+rk4/rk4_stability.m	Wed Apr 10 22:11:11 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,58 +0,0 @@
-function rk_stability()
-    ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4));
-    circ  = @(z)(abs(z));
-
-
-    % contour(X,Y,z)
-    ax = [-4 2 -3 3];
-    % hold on
-    fcontour(ruku4,[1,1],[-3, 0.6],[-3.2, 3.2])
-    hold on
-    r = 2.6;
-    fcontour(circ,[r,r],[-3, 0.6],[-3.2, 3.2],'r')
-    hold off
-    % contour(X,Y,z,[1,1],'b')
-    axis(ax)
-    title('4th order Runge-Kutta stability region')
-    xlabel('Re')
-    ylabel('Im')
-    axis equal
-    grid on
-    box on
-    hold off
-    % surf(X,Y,z)
-
-
-    rk4roots()
-end
-
-function fcontour(f,levels,x_lim,y_lim,opt)
-    default_arg('opt','b')
-    x = linspace(x_lim(1),x_lim(2));
-    y = linspace(y_lim(1),y_lim(2));
-    [X,Y] = meshgrid(x,y);
-    mu = X+ 1i*Y;
-
-    z = f(mu);
-
-    contour(X,Y,z,levels,opt)
-
-end
-
-
-function rk4roots()
-    ruku4 = @(z)(abs(1 + z +(1/2)*z.^2 + (1/6)*z.^3 + (1/24)*z.^4));
-    % Roots for real evalues:
-    F = @(x)(abs(ruku4(x))-1);
-    real_x = fzero(F,-3);
-
-    % Roots for imaginary evalues:
-    F = @(x)(abs(ruku4(1i*x))-1);
-    imag_x1 = fzero(F,-3);
-    imag_x2 = fzero(F,3);
-
-
-    fprintf('Real x = %f\n',real_x)
-    fprintf('Imag x = %f\n',imag_x1)
-    fprintf('Imag x = %f\n',imag_x2)
-end
\ No newline at end of file
--- a/+time/+rk4/rungekutta_4.m	Wed Apr 10 22:11:11 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,10 +0,0 @@
-% Takes one time step of size k using the rungekutta method
-% starting from v_0 and where the function F(v,t) gives the
-% time derivatives.
-function v = rungekutta_4(v, t , k, F)
-    k1 = F(v         ,t      );
-    k2 = F(v+0.5*k*k1,t+0.5*k);
-    k3 = F(v+0.5*k*k2,t+0.5*k);
-    k4 = F(v+    k*k3,t+    k);
-    v = v + (1/6)*(k1+2*(k2+k3)+k4)*k;
-end
\ No newline at end of file
--- a/+time/+rk4/rungekutta_6.m	Wed Apr 10 22:11:11 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-% Takes one time step of size k using the rungekutta method
-% starting from v_0 and where the function F(v,t) gives the
-% time derivatives.
-function v = rungekutta_6(v, t , k, F)
-    s = 7
-    k = zeros(length(v),s)
-    a = zeros(7,6);
-    c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1];
-    b = [1/12, 0, 0, 0, 5/12, 5/12, 1/12];
-    a = [
-        0,                           0,                          0,                       0,                     0,                 0;
-        4/7,                         0,                          0,                       0,                     0,                 0;
-        115/112,                     -5/16,                      0,                       0,                     0,                 0;
-        589/630,                     5/18,                       -16/45,                  0,                     0,                 0;
-        229/1200 - 29/6000*sqrt(5),  119/240 - 187/1200*sqrt(5), -14/75 + 34/375*sqrt(5), -3/100*sqrt(5),        0,                 0;
-        71/2400 - 587/12000*sqrt(5), 187/480 - 391/2400*sqrt(5), -38/75 + 26/375*sqrt(5), 27/80 - 3/400*sqrt(5), (1+sqrt(5))/4,     0;
-        -49/480 + 43/160*sqrt(5),    -425/96 + 51/32*sqrt(5),    52/15 - 4/5*sqrt(5),     -27/16 + 3/16*sqrt(5), 5/4 - 3/4*sqrt(5), 5/2 - 1/2*sqrt(5);
-    ]
-
-    for i = 1:s
-        u = v
-        for j = 1: i-1
-            u = u + h*a(i,j) * k(:,j)
-        end
-        k(:,i) = F(t+c(i)*k,u)
-    end
-
-    for i = 1:s
-        v = v + k*b(i)*k(:,i)
-    end
-end
--- a/+time/Rk4SecondOrderNonlin.m	Wed Apr 10 22:11:11 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,77 +0,0 @@
-classdef Rk4SecondOrderNonlin < time.Timestepper
-    properties
-        F
-        k
-        t
-        w
-        m
-
-        D
-        E
-        S
-
-        n
-    end
-
-
-    methods
-        function obj = Rk4SecondOrderNonlin(D, E, S, k, t0, v0, v0t)
-            default_arg('S',0);
-            default_arg('E',0);
-
-            if isnumeric(S)
-                S = @(v,t)S;
-            end
-
-            if isnumeric(E)
-                E = @(v)E;
-            end
-
-            obj.k = k;
-            obj.t = t0;
-            obj.w = [v0; v0t];
-
-            m = length(v0);
-            function wt = F(w,t)
-                v  = w(1:m);
-                vt = w(m+1:end);
-
-                % Def: w = [v; vt]
-                wt(1:m,1) = vt;
-                wt(m+1:2*m,1) = D(v)*v + E(v)*vt + S(v,t);
-
-            end
-
-            obj.F = @F;
-            obj.D = D;
-            obj.E = E;
-            obj.S = S;
-            obj.m = m;
-            obj.n = 0;
-        end
-
-        function [v,t] = getV(obj)
-            v = obj.w(1:end/2);
-            t = obj.t;
-        end
-
-        function [vt,t] = getVt(obj)
-            vt = obj.w(end/2+1:end);
-            t = obj.t;
-        end
-
-        function obj = step(obj)
-            obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
-            obj.t = obj.t + obj.k;
-            obj.n = obj.n + 1;
-        end
-    end
-
-
-    methods (Static)
-        function k = getTimeStep(lambda)
-            k = rk4.get_rk4_time_step(lambda);
-        end
-    end
-
-end
\ No newline at end of file
--- a/+time/Rungekutta4.m	Wed Apr 10 22:11:11 2019 +0200
+++ b/+time/Rungekutta4.m	Wed Apr 10 22:22:46 2019 +0200
@@ -1,36 +1,23 @@
 classdef Rungekutta4 < time.Timestepper
     properties
-        D
-        S
         F
-        k
+        dt
         t
         v
-        m
         n
     end
 
 
     methods
-        function obj = Rungekutta4(D, S, k, t0, v0)
-            obj.D = D;
-            obj.k = k;
+        % Create a time stepper for
+        %   v_t = F(t,v),  v(t0) = v0
+        % with step size dt.
+        function obj = Rungekutta4(F, dt, t0, v0)
+            obj.F = F;
+            obj.dt = dt;
             obj.t = t0;
             obj.v = v0;
-            obj.m = length(v0);
             obj.n = 0;
-
-            if S == 0
-                obj.S = zeros(obj.m,1);
-            else
-                obj.S = S;
-            end
-
-            if S == 0
-                obj.F = @(v,t)(obj.D*v);
-            else
-                obj.F = @(v,t)(obj.D*v + obj.S);
-            end
         end
 
         function [v,t] = getV(obj)
@@ -39,17 +26,24 @@
         end
 
         function obj = step(obj)
-            obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F);
-            obj.t = obj.t + obj.k;
+            t = obj.t;
+            v = obj.v;
+            dt = obj.dt;
+
+            k1 = obj.F(t, v);
+            k2 = obj.F(t + 0.5*dt, v + 0.5*dt*k1);
+            k3 = obj.F(t + 0.5*dt, v + 0.5*dt*k2);
+            k4 = obj.F(t + dt, v + dt*k3);
+
+            obj.v = v + dt*(1/6)*(k1+2*(k2+k3)+k4);
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
 
-
     methods (Static)
-        function k = getTimeStep(lambda)
-            k = rk4.get_rk4_time_step(lambda);
+        function dt = getTimeStep(lambda)
+            dt = rk4.get_rk4_time_step(lambda);
         end
     end
-
 end
\ No newline at end of file
--- a/+time/Rungekutta4SecondOrder.m	Wed Apr 10 22:11:11 2019 +0200
+++ b/+time/Rungekutta4SecondOrder.m	Wed Apr 10 22:22:46 2019 +0200
@@ -1,115 +1,57 @@
 classdef Rungekutta4SecondOrder < time.Timestepper
     properties
         F
-        k
-        t
-        w
-        m
-        D
-        E
-        S
-        M
-        C
-        n
+        dt
+        t, n
+        v, v_t
     end
 
 
     methods
-        % Solves u_tt = Du + Eu_t + S by
-        % Rewriting on first order form:
-        %   w_t = M*w + C(t)
-        % where
-        %   M = [
-        %      0, I;
-        %      D, E;
-        %   ]
-        % and
-        %   C(t) = [
-        %      0;
-        %      S(t)
-        %   ]
-        % D, E, S can either all be constants or all be function handles,
-        % They can also be omitted by setting them equal to the empty matrix.
-        function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
-            obj.D = D;
-            obj.E = E;
-            obj.S = S;
-            obj.m = length(v0);
+        % Create a time stepper for
+        %   v_tt = F(t,v,v_t),  v(t0) = v0, v_t(t0) = v0t
+        % with step size dt, by rewriting on first order form
+        function obj = Rungekutta4SecondOrder(F, dt, t0, v0, v0t)
+            obj.F = F;
+            obj.dt = dt;
+            obj.t = t0;
             obj.n = 0;
 
-
-            if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle')
-                default_arg('D', @(t)sparse(obj.m, obj.m));
-                default_arg('E', @(t)sparse(obj.m, obj.m));
-                default_arg('S', @(t)sparse(obj.m, 1)    );
-
-                if ~isa(D, 'function_handle')
-                    D = @(t)D;
-                end
-                if ~isa(E, 'function_handle')
-                    E = @(t)E;
-                end
-                if ~isa(S, 'function_handle')
-                    S = @(t)S;
-                end
-
-                obj.k = k;
-                obj.t = t0;
-                obj.w = [v0; v0t];
-
-                % Avoid matrix formulation because it is VERY slow
-                obj.F = @(w,t)[
-                    w(obj.m+1:end);
-                    D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t);
-                ];
-            else
-
-                default_arg('D', sparse(obj.m, obj.m));
-                default_arg('E', sparse(obj.m, obj.m));
-                default_arg('S', sparse(obj.m, 1)    );
-
-                I = speye(obj.m);
-                O = sparse(obj.m,obj.m);
-
-                obj.M = [
-                    O, I;
-                    D, E;
-                ];
-                obj.C = [
-                    zeros(obj.m,1);
-                                 S;
-                ];
-
-                obj.k = k;
-                obj.t = t0;
-                obj.w = [v0; v0t];
-
-                obj.F = @(w,t)(obj.M*w + obj.C);
-            end
+            obj.v = v0;
+            obj.v_t = v0t;
         end
 
         function [v,t] = getV(obj)
-            v = obj.w(1:end/2);
+            v = obj.v
             t = obj.t;
         end
 
         function [vt,t] = getVt(obj)
-            vt = obj.w(end/2+1:end);
+            vt = obj.v_t;
             t = obj.t;
         end
 
         function obj = step(obj)
-            obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
-            obj.t = obj.t + obj.k;
+            t = obj.t;
+            v = obj.v;
+            v_t = obj.v_t;
+            dt = obj.dt;
+
+            k1 = obj.F(t,          v,                            v_t);
+            k2 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t,               v_t + 1/2*dt*k1);
+            k3 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t + 1/4*dt^2*k1, v_t + 1/2*dt*k2);
+            k4 = obj.F(t + dt,     v + dt*v_t + 1/2*dt^2*k2,     v_t + dt*k3);
+
+            obj.v = v + dt*v_t + dt^2*(1/6)*(k1 + k2 + k3);
+            obj.v_t = v_t + dt*(1/6)*(k1 + 2*k2 + 2*k3 + k4);
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end
 
-
     methods (Static)
         function k = getTimeStep(lambda)
             k = rk4.get_rk4_time_step(lambda);
         end
     end
-
-end
\ No newline at end of file
+end
--- a/+time/Rungekutta4proper.m	Wed Apr 10 22:11:11 2019 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,42 +0,0 @@
-classdef Rungekutta4proper < time.Timestepper
-    properties
-        F
-        k
-        t
-        v
-        m
-        n
-    end
-
-
-    methods
-        % Timesteps v_t = F(v,t), using RK4 fromt t = t0 with timestep k and initial conditions v = v0
-        function obj = Rungekutta4proper(F, k, t0, v0)
-            obj.F = F;
-            obj.k = k;
-            obj.t = t0;
-            obj.v = v0;
-            obj.m = length(v0);
-            obj.n = 0;
-        end
-
-        function [v,t] = getV(obj)
-            v = obj.v;
-            t = obj.t;
-        end
-
-        function obj = step(obj)
-            obj.v = time.rk4.rungekutta_4(obj.v, obj.t, obj.k, obj.F);
-            obj.t = obj.t + obj.k;
-            obj.n = obj.n + 1;
-        end
-    end
-
-
-    methods (Static)
-        function k = getTimeStep(lambda)
-            k = rk4.get_rk4_time_step(lambda);
-        end
-    end
-
-end
\ No newline at end of file
--- a/+time/Timestepper.m	Wed Apr 10 22:11:11 2019 +0200
+++ b/+time/Timestepper.m	Wed Apr 10 22:22:46 2019 +0200
@@ -1,13 +1,14 @@
 classdef Timestepper < handle
     properties (Abstract)
         t
-        k
+        k % TBD: should this be a method instead?
         n
     end
 
     methods (Abstract)
-         [v,t] = getV(obj)
-         obj = step(obj)
+        % Returns the solution vector v at timestep t.
+        [v,t] = getV(obj)
+        obj = step(obj)
     end
 
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/blargh.m	Wed Apr 10 22:22:46 2019 +0200
@@ -0,0 +1,20 @@
+% Takes nSteps with ts and evaluates the provided function handles
+% for each step and collects the output
+function varargout = blargh(ts, nSteps, varargin)
+    N = length(varargin);
+
+    for j = 1:N
+        Fj = varargin{j}(ts);
+        varargout{j} = zeros(length(Fj),nSteps+1);
+        varargout{j}(:,1) = Fj;
+    end
+
+    for i = 2:nSteps+1
+        ts.step();
+
+        for j = 1:N
+            varargout{j}(:, i) = varargin{j}(ts);
+        end
+    end
+end
+% TODO: Name and discuss this funcition