Mercurial > repos > public > sbplib

comparison +time/Rungekutta.m @ 918:679f4ddd982f feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add properties for stage approximations and stage rates in the Rungekutta class.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 26 Nov 2018 16:23:27 -0800
parents 8732d6bd9890
children 384ca2331a12
comparison
equal deleted inserted replaced
917:878652b22157 918:679f4ddd982f
1 classdef Rungekutta < time.Timestepper 1 classdef Rungekutta < time.Timestepper
2 properties 2 properties
3 F % RHS of the ODE 3 F % RHS of the ODE
4 k % Time step 4 dt % Time step
5 t % Time point 5 t % Time point
6 v % Solution vector 6 v % Solution vector
7 n % Time level 7 n % Time level
8 scheme % The scheme used for the time stepping, e.g rk4, rk6 etc. 8 scheme % The scheme used for the time stepping, e.g rk4, rk6 etc.
9 coeffs % Butcher tableau coefficients
10 V % All stage approximations in most recent time step
11 K % All stage rates in most recent time step
9 end 12 end
10 13
11 14
12 methods 15 methods
13 % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with 16 % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with
14 % timestep k and initial conditions v = v0 17 % timestep dt and initial conditions v = v0
15 function obj = Rungekutta(F, k, t0, v0, method) 18 function obj = Rungekutta(F, dt, t0, v0, method)
16 default_arg('method',"rk4"); 19 default_arg('method',"rk4");
17 obj.F = F; 20 obj.F = F;
18 obj.k = k; 21 obj.dt = dt;
19 obj.t = t0; 22 obj.t = t0;
20 obj.v = v0; 23 obj.v = v0;
21 obj.n = 0; 24 obj.n = 0;
25
26 % Extract the coefficients for the specified method
27 % used for the RK updates from the Butcher tableua.
28 [s,a,b,c] = time.rk.butcherTableau(method);
29 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
30
22 % TODO: method "rk4" is also implemented in the butcher tableau, but the rungekutta_4.m implementation 31 % TODO: method "rk4" is also implemented in the butcher tableau, but the rungekutta_4.m implementation
23 % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it. 32 % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it.
24 if (method == "rk4") 33 if (method == "rk4")
25 obj.scheme = @time.rk.rungekutta_4; 34 obj.scheme = @time.rk.rungekutta_4;
26 else 35 else
27 % Extract the coefficients for the specified method 36 obj.scheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, obj.coeffs);
28 % used for the RK updates from the Butcher tableua.
29 [s,a,b,c] = time.rk.butcherTableau(method);
30 coeffs = struct('s',s,'a',a,'b',b,'c',c);
31 obj.scheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs);
32 end 37 end
33 end 38 end
34 39
35 function [v,t] = getV(obj) 40 % v: Current solution
41 % t: Current time
42 % V: All stage approximations in most recent time step
43 % K: All stage rates in most recent time step
44 % T: Time points (corresponding to V and K) in most recent time step
45 function [v,t,V,T,K] = getV(obj)
36 v = obj.v; 46 v = obj.v;
37 t = obj.t; 47 t = obj.t;
48 V = obj.V;
49 K = obj.K;
50 T = obj.t + obj.dt*obj.coeffs.b;
38 end 51 end
39 52
40 function obj = step(obj) 53 function obj = step(obj)
41 obj.v = obj.scheme(obj.v, obj.t, obj.k, obj.F); 54 [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.dt, obj.F);
42 obj.t = obj.t + obj.k; 55 obj.t = obj.t + obj.dt;
43 obj.n = obj.n + 1; 56 obj.n = obj.n + 1;
44 end 57 end
45 end 58 end
46 end 59 end