Mercurial > repos > public > sbplib

comparison +time/+rk/Explicit.m @ 996:3b903011b1a9 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename time.rk.General to time.rk.Explicit and fix some errors
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Jan 2019 23:01:17 +0100
parents +time/+rk/General.m@10c5eda235b7
children d4fe089b2c4a
comparison
equal deleted inserted replaced
995:10c5eda235b7 996:3b903011b1a9
1 classdef Explicit < time.Timestepper
2 properties
3 F % RHS of the ODE
4 dt % Time step
5 t % Time point
6 v % Solution vector
7 n % Time level
8 scheme % The scheme used for the time stepping, e.g rk4, rk6 etc.
9 bt
10 V % All stage approximations in most recent time step
11 K % All stage rates in most recent time step
12 end
13
14
15 methods
16 % Timesteps v_t = F(t,v), using the specified ButcherTableau
17 % from t = t0 with timestep dt and initial conditions v(0) = v0
18 function obj = Explicit(F, dt, t0, v0, bt)
19 assertType(bt, 'time.rk.ButcherTableau')
20 obj.F = F;
21 obj.dt = dt;
22 obj.t = t0;
23 obj.v = v0;
24 obj.n = 0;
25
26 assert(bt.isExplicit())
27 obj.bt = bt;
28 end
29
30 % v: Current solution
31 % t: Current time
32 % V: All stage approximations in most recent time step
33 % K: All stage rates in most recent time step
34 % T: Time points (corresponding to V and K) in most recent time step
35 function [v,t] = getV(obj)
36 v = obj.v;
37 t = obj.t;
38 end
39
40 function obj = step(obj)
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;
103 obj.t = obj.t + obj.dt;
104 obj.n = obj.n + 1;
105 end
106
107 % Returns a vector of time points, including substage points,
108 % in the time interval [t0, tEnd].
109 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
110 function tvec = timePoints(obj, t0, tEnd)
111 % TBD: Should this be implemented here or somewhere else?
112 N = round( (tEnd-t0)/obj.dt );
113 tvec = zeros(N*obj.s, 1);
114 s = obj.coeffs.s;
115 c = obj.coeffs.c;
116 for i = 1:N
117 ind = (i-1)*s+1 : i*s;
118 tvec(ind) = ((i-1) + c')*obj.dt;
119 end
120 end
121
122 % Returns a vector of quadrature weights corresponding to grid points
123 % in time interval [t0, tEnd], substage points included.
124 % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
125 function weights = quadWeights(obj, t0, tEnd)
126 % TBD: Should this be implemented here or somewhere else?
127 N = round( (tEnd-t0)/obj.dt );
128 b = obj.coeffs.b;
129 weights = repmat(b', N, 1);
130 end
131 end
132
133 methods(Static)
134 % TBD: Function name
135 function ts = methodFromStr(F, dt, t0, v0, methodStr)
136 try
137 bt = time.rk.ButcherTableau.(method);
138 catch
139 error('Runge-Kutta method ''%s'' is not implemented', methodStr)
140 end
141
142 ts = time.rk.Explicit(F, dt, t0, v0, bt);
143 end
144 end
145 end