Mercurial > repos > public > sbplib

comparison +time/Rungekutta4SecondOrder.m @ 991:a99f00896b8e feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Rungekutta4SecondOrder native second order
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Jan 2019 10:17:00 +0100
parents ad6de007e7f6
children f2988a63c3aa
comparison
equal deleted inserted replaced
990:1066bb31bc95 991:a99f00896b8e
1 classdef Rungekutta4SecondOrder < time.Timestepper 1 classdef Rungekutta4SecondOrder < time.Timestepper
2 properties 2 properties
3 F 3 F
4 dt 4 dt
5 t 5 t, n
6 n 6 v, v_t
7
8 G % RHS for rewritten system
9 w
10 end 7 end
11 8
12 9
13 methods 10 methods
14 % Create a time stepper for 11 % Create a time stepper for
18 obj.F = F; 15 obj.F = F;
19 obj.dt = dt; 16 obj.dt = dt;
20 obj.t = t0; 17 obj.t = t0;
21 obj.n = 0; 18 obj.n = 0;
22 19
23 m = length(v0); 20 obj.v = v0;
24 obj.w = [v0; v0t]; 21 obj.v_t = v0t;
25 obj.G = @(t, w)[
26 w(m+1:end);
27 F(t,w(1:m), w(m+1:end));
28 ];
29 end 22 end
30 23
31 function [v,t] = getV(obj) 24 function [v,t] = getV(obj)
32 v = obj.w(1:end/2); 25 v = obj.v
33 t = obj.t; 26 t = obj.t;
34 end 27 end
35 28
36 function [vt,t] = getVt(obj) 29 function [vt,t] = getVt(obj)
37 vt = obj.w(end/2+1:end); 30 vt = obj.v_t;
38 t = obj.t; 31 t = obj.t;
39 end 32 end
40 33
41 function obj = step(obj) 34 function obj = step(obj)
42 % TODO: Avoid extra storage 35 t = obj.t;
43 w = obj.w; 36 v = obj.v;
37 v_t = obj.v_t;
44 dt = obj.dt; 38 dt = obj.dt;
45 39
46 k1 = obj.G(t, w); 40 k1 = obj.F(t, v, v_t);
47 k2 = obj.G(t + 0.5*dt, w + 0.5*dt*k1); 41 k2 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t, v_t + 1/2*dt*k1);
48 k3 = obj.G(t + 0.5*dt, w + 0.5*dt*k2); 42 k3 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t + 1/4*dt^2*k1, v_t + 1/2*dt*k2);
49 k4 = obj.G(t + dt, w + dt*k3); 43 k4 = obj.F(t + dt, v + dt*v_t + 1/2*dt^2*k2, v_t + dt*k3);
50 44
51 obj.w = w + dt*(1/6)*(k1+2*(k2+k3)+k4); 45 obj.v = v + dt*v_t + dt^2*(1/6)*(k1 + k2 + k3);
46 obj.v_t = v_t + dt*(1/6)*(k1 + 2*k2 + 2*k3 + k4);
52 obj.t = obj.t + obj.dt; 47 obj.t = obj.t + obj.dt;
53 obj.n = obj.n + 1; 48 obj.n = obj.n + 1;
54 end 49 end
55 end 50 end
56
57 51
58 methods (Static) 52 methods (Static)
59 function k = getTimeStep(lambda) 53 function k = getTimeStep(lambda)
60 k = rk4.get_rk4_time_step(lambda); 54 k = rk4.get_rk4_time_step(lambda);
61 end 55 end
62 end 56 end
63
64 end 57 end