Mercurial > repos > public > sbplib

comparison +time/Rungekutta.m @ 931:384ca2331a12 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Rungekutta class allow for discrete data.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 03 Dec 2018 16:26:44 -0800
parents 679f4ddd982f
children 3860dad28239
comparison
equal deleted inserted replaced
930:34d882bffae4 931:384ca2331a12
13 13
14 14
15 methods 15 methods
16 % 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
17 % timestep dt and initial conditions v = v0 17 % timestep dt and initial conditions v = v0
18 function obj = Rungekutta(F, dt, t0, v0, method) 18 function obj = Rungekutta(F, dt, t0, v0, method, discreteData)
19 default_arg('method',"rk4"); 19 default_arg('method',"rk4");
20 default_arg('discreteData', []);
20 obj.F = F; 21 obj.F = F;
21 obj.dt = dt; 22 obj.dt = dt;
22 obj.t = t0; 23 obj.t = t0;
23 obj.v = v0; 24 obj.v = v0;
24 obj.n = 0; 25 obj.n = 0;
26 % Extract the coefficients for the specified method 27 % Extract the coefficients for the specified method
27 % used for the RK updates from the Butcher tableua. 28 % used for the RK updates from the Butcher tableua.
28 [s,a,b,c] = time.rk.butcherTableau(method); 29 [s,a,b,c] = time.rk.butcherTableau(method);
29 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c); 30 obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
30 31
31 % TODO: method "rk4" is also implemented in the butcher tableau, but the rungekutta_4.m implementation 32 if isempty(discreteData)
32 % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it. 33 % TODO: method "rk4" is also implemented in the butcher tableau, but the rungekutta_4.m implementation
33 if (method == "rk4") 34 % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it.
34 obj.scheme = @time.rk.rungekutta_4; 35 if (method == "rk4")
36 obj.scheme = @(v,t,dt,F,n) time.rk.rungekutta_4(v ,t, dt, F);
37 else
38 obj.scheme = @(v,t,dt,F,n) time.rk.rungekutta(v, t, dt, F, obj.coeffs);
39 end
35 else 40 else
36 obj.scheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, obj.coeffs); 41 obj.scheme = @(v,t,dt,F,n) time.rk.rungekuttaDiscreteData(v, t, dt, F, obj.coeffs, discreteData, n);
37 end 42 end
38 end 43 end
39 44
40 % v: Current solution 45 % v: Current solution
41 % t: Current time 46 % t: Current time
49 K = obj.K; 54 K = obj.K;
50 T = obj.t + obj.dt*obj.coeffs.b; 55 T = obj.t + obj.dt*obj.coeffs.b;
51 end 56 end
52 57
53 function obj = step(obj) 58 function obj = step(obj)
54 [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.dt, obj.F); 59 [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.dt, obj.F, obj.n);
55 obj.t = obj.t + obj.dt; 60 obj.t = obj.t + obj.dt;
56 obj.n = obj.n + 1; 61 obj.n = obj.n + 1;
57 end 62 end
58 end 63 end
59 end 64 end