Mercurial > repos > public > sbplib

comparison +time/+rk/General.m @ 995:10c5eda235b7 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Full use of butcher tableau in time.rk.General. Inline rungekutta step methods
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Jan 2019 22:57:13 +0100
parents 44e7e497c3b7
children
comparison
equal deleted inserted replaced
994:2f89959fb9f0 995:10c5eda235b7
4 dt % 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 9 bt
10 V % All stage approximations in most recent time step 10 V % All stage approximations in most recent time step
11 K % All stage rates in most recent time step 11 K % All stage rates in most recent time step
12 end 12 end
13 13
14 14
15 methods 15 methods
16 % Timesteps v_t = F(v,t), using the specified ButcherTableau 16 % Timesteps v_t = F(t,v), using the specified ButcherTableau
17 % from t = t0 with timestep dt and initial conditions v(0) = v0 17 % from t = t0 with timestep dt and initial conditions v(0) = v0
18 function obj = General(F, dt, t0, v0, bt, discreteData) 18 function obj = General(F, dt, t0, v0, bt)
19 assertType(bt, 'time.rk.ButcherTableau') 19 assertType(bt, 'time.rk.ButcherTableau')
20 default_arg('discreteData', []);
21 obj.F = F; 20 obj.F = F;
22 obj.dt = dt; 21 obj.dt = dt;
23 obj.t = t0; 22 obj.t = t0;
24 obj.v = v0; 23 obj.v = v0;
25 obj.n = 0; 24 obj.n = 0;
26 25
27 assert(bt.isExplicit()) 26 assert(bt.isExplicit())
28 obj.coeffs = struct('s',bt.nStages,'a',bt.a(:,1:end-1),'b',bt.b,'c',bt.c); 27 obj.bt = bt;
29
30 if isempty(discreteData)
31 obj.scheme = @(v,t,n) time.rk.rungekutta(v, t, dt, F, obj.coeffs);
32 else
33 obj.scheme = @(v,t,n) time.rk.rungekuttaDiscreteData(v, t, dt, F, obj.coeffs, discreteData, n);
34 end
35 end 28 end
36 29
37 % v: Current solution 30 % v: Current solution
38 % t: Current time 31 % t: Current time
39 % V: All stage approximations in most recent time step 32 % V: All stage approximations in most recent time step
40 % K: All stage rates in most recent time step 33 % K: All stage rates in most recent time step
41 % T: Time points (corresponding to V and K) in most recent time step 34 % T: Time points (corresponding to V and K) in most recent time step
42 function [v,t,V,T,K] = getV(obj) 35 function [v,t] = getV(obj)
43 v = obj.v; 36 v = obj.v;
44 t = obj.t; 37 t = obj.t;
45 V = obj.V;
46 K = obj.K;
47 T = obj.t + obj.dt*obj.coeffs.b;
48 end 38 end
49 39
50 function obj = step(obj) 40 function obj = step(obj)
51 [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.n); 41 s = obj.bt.nStages();
42 a = obj.bt.a;
43 b = obj.bt.b;
44 c = obj.bt.c;
45
46 % Compute rates K
47 K = zeros(length(v), s);
48 for i = 1:s
49 V_i = obj.v;
50 for j = 1:i-1
51 V_i = V_i + dt*a(i,j)*K(:,j);
52 end
53 K(:,i) = F(t+dt*c(i), V_i);
54 end
55
56 % Compute updated solution
57 v_next = v;
58 for i = 1:s
59 v_next = v_next + dt*b(i)*K(:,i);
60 end
61
62 obj.v = v_next;
63 obj.t = obj.t + obj.dt;
64 obj.n = obj.n + 1;
65 end
66
67 % TBD: Method name
68 % TBD: Parameter name
69 %
70 % Takes a regular step but with discreteRates(:,i) added to RHS for stage i.
71 % v_t = F(t,v) + discreteRates(:, ...)
72 %
73 % Also returns the stage approximations (V) and stage rates (K).
74 function [v,t, V, K] = stepWithDiscreteData(obj, discreteRates)
75 s = obj.bt.nStages();
76 a = obj.bt.a;
77 b = obj.bt.b;
78 c = obj.bt.c;
79
80 % Compute rates K and stage approximations V
81 K = zeros(length(v), s);
82 V = zeros(length(v), s);
83 for i = 1:s
84 V_i = obj.v;
85 for j = 1:i-1
86 V_i = V_i + dt*a(i,j)*K(:,j);
87 end
88
89 K_i = F(t+dt*c(i), V_i);
90 K_i = K_i + discreteRates(:,i);
91
92 V(:,i) = V_i;
93 K(:,i) = K_i;
94 end
95
96 % Compute updated updated solution
97 v_next = v;
98 for i = 1:s
99 v_next = v_next + dt*b(i)*K(:,i);
100 end
101
102 obj.v = v_next;
52 obj.t = obj.t + obj.dt; 103 obj.t = obj.t + obj.dt;
53 obj.n = obj.n + 1; 104 obj.n = obj.n + 1;
54 end 105 end
55 106
56 % Returns a vector of time points, including substage points, 107 % Returns a vector of time points, including substage points,
57 % in the time interval [t0, tEnd]. 108 % in the time interval [t0, tEnd].
58 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already. 109 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
59 function tvec = timePoints(obj, t0, tEnd) 110 function tvec = timePoints(obj, t0, tEnd)
111 % TBD: Should this be implemented here or somewhere else?
60 N = round( (tEnd-t0)/obj.dt ); 112 N = round( (tEnd-t0)/obj.dt );
61 tvec = zeros(N*obj.s, 1); 113 tvec = zeros(N*obj.s, 1);
62 s = obj.coeffs.s; 114 s = obj.coeffs.s;
63 c = obj.coeffs.c; 115 c = obj.coeffs.c;
64 for i = 1:N 116 for i = 1:N
69 121
70 % Returns a vector of quadrature weights corresponding to grid points 122 % Returns a vector of quadrature weights corresponding to grid points
71 % in time interval [t0, tEnd], substage points included. 123 % in time interval [t0, tEnd], substage points included.
72 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already. 124 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
73 function weights = quadWeights(obj, t0, tEnd) 125 function weights = quadWeights(obj, t0, tEnd)
126 % TBD: Should this be implemented here or somewhere else?
74 N = round( (tEnd-t0)/obj.dt ); 127 N = round( (tEnd-t0)/obj.dt );
75 b = obj.coeffs.b; 128 b = obj.coeffs.b;
76 weights = repmat(b', N, 1); 129 weights = repmat(b', N, 1);
77 end 130 end
78 end 131 end