Mercurial > repos > public > sbplib
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