view +time/Rungekutta4SecondOrder.m @ 985:ad6de007e7f6 feature/timesteppers

Convert Rungekutta4SecondOrder to take F(t,v,v_t) instead of matrices
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 08 Jan 2019 12:34:29 +0100
parents 0585a2ee7ee7
children a99f00896b8e
line wrap: on
line source

classdef Rungekutta4SecondOrder < time.Timestepper
    properties
        F
        dt
        t
        n

        G % RHS for rewritten system
        w
    end


    methods
        % 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;

            m = length(v0);
            obj.w = [v0; v0t];
            obj.G = @(t, w)[
                w(m+1:end);
                F(t,w(1:m), w(m+1:end));
            ];
        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)
            % TODO: Avoid extra storage
            w = obj.w;
            dt = obj.dt;

            k1 = obj.G(t, w);
            k2 = obj.G(t + 0.5*dt, w + 0.5*dt*k1);
            k3 = obj.G(t + 0.5*dt, w + 0.5*dt*k2);
            k4 = obj.G(t + dt, w + dt*k3);

            obj.w = w + 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);
        end
    end

end