4 %
function, which performs a reduced basis online simulation
for
5 % the parameter vector mu, which is assumed to be set in the model
6 %
class, i.e. a previous model = model.set_mu(model,[...]) is assumed
9 % allowed dependency of data: Nmax, N, M, mu
10 % not allowed dependency of data: H
11 % allowed dependency of computation: Nmax, N, M, mu
12 % not allowed dependency of computation: H
13 % Unknown at
this stage: ---
15 % Required fields of model as required by the numerical operators
16 % mu_names : the cell array of names of mu-components and
17 %
for each of these stringt, a corresponding
18 % field in model is expected.
19 % T : end-time of simulation
20 % nt : number of timesteps to compute
22 % L_I_inv_norm_bound: upper bound on implicit
operator inverse norm
23 % constant in time, a single scalar
24 % L_E_norm_bound: upper bounds on
explicit operator
25 % constant in time, a single scalar
27 % optional fields of model:
28 % data_const_in_time :
if this flag is set, then only operators
29 %
for first time instance are computed
30 % name_output_functional:
if this field is existent, then an
31 % output estimation is performed and error.estimtations
32 % starting_tim_step: in t-partition
case this is the starting time step
33 %
for the actual t-partitioni time step
34 % stopping_time_step: in t-partition
this is the stopping time step
for
35 % the actual t-partition
37 % optional fileds of reduced_data:
38 % Delta0: initial error
40 % generated fields of simulation_data:
41 % a: time sequence of solution coefficients, columns are
42 %
'a(:,k)' = `a^(k-1)`
43 % Delta:time sequence of `L^2`-posteriori error estimates
44 %
'Delta(k)' = `\Delta^{k-1}_N`
45 % or energy-norm-posterior error estimates
46 %
'Delta_energy(k)' = `\bar \Delta^{k-1}_N`
47 % depending on the field
"error_norm" in model.
49 %
if model.name_output_functional is set, then additionally, a
50 % sequence of output estimates
's(U(:,))' and error bound Delta_s is returned
52 % see the ***_gen_reduced_data() routine for specifications of the
53 % fields of reduced_data
55 % Bernard Haasdonk 15.5.2007
56 % Markus Dihlmann Feb 2011
58 % extended to primal/dual formulaton by Dominik Garmatter:
60 % There are now 3 controlfields (model.want_dual, model.want_impproved_output and model.use_scm).
61 % These enable the following cases:
64 % - model.want_improved_output = 0 and model.want_dual = 1 or 0:
65 % => lin_evol_rb_simulation_data_primal_dual(model, reduced_data) with
66 % reduced_data including the primal or dual reduced_data. Simply performs a
67 % dual or primal rb_simulation. In the dual case a functional has to be
68 % specified via model.name_output_functional. In the dual case there is
69 % never an output. In the primal case one can choose to cpmoute the reduced
70 % output (and estimator) by specifying a functional via
71 % model.name_output_functional.
73 %- model.want_improved_output = 1 and model.want_dual = 0:
75 % reduced_data including the primal and to some extend dual reduced_data!
80 % In all the cases it is possible to choose model.error_norm to be 'l2'
81 % (standard) or 'energy'. In the 'energy'-case a lower bound to the
82 % stability constant (coercivity or infsup) is required. Either it is set
83 % in model.constant_LB or one uses the SCM (model.use_scm = 1). In case you
84 % want to use SCM you have to set model.error_norm = 'energy' and
85 % model.use_scm = 1 before generating reduced_data. Also see
scm_offline
86 % for other required quantities.
87 % In the above case (model.want_improved_output = 1 and model.want_dual = 0)
88 % the energy error estimator is always used.
90 % special Note: in the european_option_pricing model a functional named
91 % 'theta' includes a partial timederivative and therefore produces a
92 % slightly different dual problem. Therefore the 'theta'-case is sometimes
93 % treated differently. Don't care about this case if you don't use the
94 % european_option_pricing model.
96 % Dominik Garmatter 03.08 2012
98 if ~isfield(model, 'want_dual')
101 if ~isfield(model, 'want_improved_output')
102 model.want_improved_output = 0;
104 if ~isfield(model, 'use_scm')
108 % if you want the improved output you need the dual rb_coefficients and the
109 % dual Delta for the estimator, so they are exracted here and the reduced_data
110 % have to include them if you want the improved output!
111 if model.want_improved_output && model.want_dual == 0
112 dual_data.a = reduced_data.dual_sim_data;
113 dual_data.Delta = reduced_data.dual_Delta;
114 if isequal(size(dual_data.Delta), [1,1]) == 0 % give a warning if the dual estimator isn't measured in the energy norm
115 warning('you might have chosen the wrong error estimator for your dual rb_simulation')
121 disp(['performing simulation for mu = ', mat2str(get_mu(model))]);
124 if ~isequal(model.error_norm,'l2') && ~isequal(model.error_norm,'energy')
125 error('error_norm unknown!!');
128 if isfield(model,'starting_time_step')
129 t_ind_start = model.starting_time_step;
134 if isfield(model, 'stopping_time_step')
135 t_ind_stop = model.stopping_time_step;
137 t_ind_stop = model.nt;
140 % number of reduced Basis vectors
141 nRB = length(reduced_data.a0{1});
143 a = zeros(nRB,t_ind_stop-t_ind_start+1); % matrix of RB-coefficients
144 Delta = zeros(1,t_ind_stop-t_ind_start+1); % vector of error estimators
145 model.dt = model.T/model.nt;
146 % initialise the correction term (only in primal
case with improved output)
147 if model.want_improved_output && model.want_dual == 0
151 % initial data projection and linear combination
152 model.decomp_mode = 2;
153 sa0 = rb_init_values(model,[]); % rb_init_values have to include the differentiation between the primal and the dual
case (see
for example
rb_init_values_primal_dual)
154 a0 = lincomb_sequence(reduced_data.a0,sa0);
156 a(:,1) = a0; % in the dual
case this will be modified in the first timestep
159 if isfield(reduced_data,
'Delta0')
160 Delta(1) = reduced_data.Delta0;
163 LL_I_is_speye = -1; % unknown state here
165 for t = (t_ind_start+1):t_ind_stop % loop over time steps: computation of a(t+1)==a^t
166 % set the time - relevant if your operators are dependant on time.
167 if model.want_dual == 0
168 model.t = (t-1)*model.dt; % forward time in primal case
170 model.t = model.T - (t-1)*model.dt; % the backwards time in dual case
174 disp(['entered time-loop step ',num2str(t)]);
177 % linear combination of all quantities - only once, if data is const in time
178 if (t==t_ind_start+1) || (~model.data_const_in_time)
179 if model.want_dual %dual
180 if strcmp(model.name_output_functional, 'theta') == 1 % in theta-case the M-1-timestep is a little different and thus the residualnorm in that timestep
181 [sLL_I, sLL_E, ~, sK_II, sK_IE, sK_EE, sm_I, sm_E, sm, sK_IdId, sK_IdE, ~, ~, ~] = ...
182 rb_operators(model,[]);
183 % additional quantities for the dual resiual in the theta-case
184 m_I = lincomb_sequence(reduced_data.m_I,sm_I);
185 m_E = lincomb_sequence(reduced_data.m_E,sm_E);
186 m = lincomb_sequence(reduced_data.m,sm);
187 else % non-theta dual case
188 [sLL_I, sLL_E, ~, sK_II, sK_IE, sK_EE, ~, ~, ~, sK_IdId, sK_IdE, ~, ~, ~] = ...
189 rb_operators(model,[]);
191 LL_I = lincomb_sequence(reduced_data.LL_I, sLL_I);
192 LL_E = lincomb_sequence(reduced_data.LL_E, sLL_E);
193 % for the dual residual
194 K_II = lincomb_sequence(reduced_data.K_II, sK_II);
195 K_IE = lincomb_sequence(reduced_data.K_IE, sK_IE);
196 K_EE = lincomb_sequence(reduced_data.K_EE, sK_EE);
197 K_IdId = lincomb_sequence(reduced_data.K_IdId, sK_IdId);
198 K_IdE = lincomb_sequence(reduced_data.K_IdE, sK_IdE);
201 [sLL_I, sLL_E, sbb, sK_II, sK_IE, sK_EE, sm_I, sm_E, sm, ~, ~, ~, ~, ~] = ...
202 rb_operators(model,[]);
203 LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
204 LL_E = lincomb_sequence(reduced_data.LL_E,sLL_E);
205 bb = lincomb_sequence(reduced_data.bb,sbb);
206 K_II = lincomb_sequence(reduced_data.K_II,sK_II);
207 K_IE = lincomb_sequence(reduced_data.K_IE,sK_IE);
208 K_EE = lincomb_sequence(reduced_data.K_EE,sK_EE);
209 m_I = lincomb_sequence(reduced_data.m_I,sm_I);
210 m_E = lincomb_sequence(reduced_data.m_E,sm_E);
211 m = lincomb_sequence(reduced_data.m,sm);
212 if model.want_improved_output % imrproved output includes the quantities for the correction term
213 [sLL_I, sLL_E, sbb, sK_II, sK_IE, sK_EE, sm_I, sm_E, sm, ~, ~, sLL_I_correct, sLL_E_correct, sbb_correct] = ...
214 rb_operators(model,[]);
215 LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
216 LL_E = lincomb_sequence(reduced_data.LL_E,sLL_E);
217 bb = lincomb_sequence(reduced_data.bb,sbb);
218 K_II = lincomb_sequence(reduced_data.K_II,sK_II);
219 K_IE = lincomb_sequence(reduced_data.K_IE,sK_IE);
220 K_EE = lincomb_sequence(reduced_data.K_EE,sK_EE);
221 m_I = lincomb_sequence(reduced_data.m_I,sm_I);
222 m_E = lincomb_sequence(reduced_data.m_E,sm_E);
223 m = lincomb_sequence(reduced_data.m,sm);
224 LL_I_correct = lincomb_sequence(reduced_data.LL_I_correct, sLL_I_correct);
225 LL_E_correct = lincomb_sequence(reduced_data.LL_E_correct, sLL_E_correct);
226 bb_correct = lincomb_sequence(reduced_data.bb_correct, sbb_correct);
230 % fix for the dual endtimecondition: up to know it is only u^{du,M} = P_{X_N}(v_l) but we need u^{du,M} = -L_I^{du,-1}P_{X_N}(v_l)
231 if (t==t_ind_start+1) && model.want_dual
232 a(:,1) = - LL_I \ a0; % a0 was exactly P_{X_N}(v_l)
235 if model.data_const_in_time
236 if model.want_dual % dual
237 % save the original L_E
for the last timestep of the dual problem
238 % (it is a(:,t_ind_stop) = L_E*a(:,t_ind_stop+1), i.e. L_I =
239 % Id in the last timestep of the dual problem)
242 if strcmp(model.name_output_functional,
'theta') == 1 % prepare the inhomogeniety of the theta-
case
245 LL_I = speye(size(LL_I));
251 LL_I = speye(size(LL_I));
254 elseif max(max(LL_I- speye(size(LL_I))))<1e-12
259 if model.want_dual % solve LL_I * a(:,t_adj) = LL_E * a(:,t_adj+1), i.e. the dual problem
260 if (t==t_ind_start+1) && strcmp(model.name_output_functional,
'theta') == 1 % slightly different dual problem in
case of the theta-functional
261 rhs = LL_E * a(:,t-t_ind_start) + a0;
262 elseif (t==t_ind_stop) % the last timestep of the dual problem
263 if model.data_const_in_time % use the saved L_E here (remember LL_I_is_speye = 1 anyway in
this case)
264 rhs = LL_E_save * a(:,t-t_ind_start);
265 else % otherwise use the current L_E and set L_I = Id.
266 rhs = LL_E * a(:,t-t_ind_start);
267 LL_I = speye(size(LL_I));
270 else % the normal timesteps of the dual problem
271 rhs = LL_E * a(:,t-t_ind_start);
274 if LL_I_is_speye == 1
275 a(:, t+1-t_ind_start) = rhs;
276 elseif LL_I_is_speye == -1 && isequal(LL_I, speye(size(LL_I)))
277 a(:, t+1-t_ind_start) = rhs;
279 [a(:, t+1-t_ind_start), flag] = bicgstab(LL_I,rhs,[],1000);
281 disp(['error in system solution, solver return flag = ', ...
286 % the norm of the dual residual
287 if strcmp(model.name_output_functional, 'theta') == 1 && (t==t_ind_start+1) % M-1-timestep in theta-case differs!
289 a(:,t+1-t_ind_start)' * K_II * a(:,t+1-t_ind_start) ...
290 -2* a(:,t+1-t_ind_start)' * K_IE * a(:,t-t_ind_start) ...
291 + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start) ...
293 - 2* a(:,t+1-t_ind_start)' * m_I ...
294 + 2* a(:,t-t_ind_start)' * m_E;
295 elseif (t==t_ind_stop)
297 a(:,t+1-t_ind_start)' * K_IdId * a(:,t+1-t_ind_start) ...
298 -2* a(:,t+1-t_ind_start)' * K_IdE * a(:,t-t_ind_start) ...
299 + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start);
302 a(:,t+1-t_ind_start)' * K_II * a(:,t+1-t_ind_start) ...
303 -2* a(:,t+1-t_ind_start)' * K_IE * a(:,t-t_ind_start) ...
304 + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start);
307 else % solve LL_I* a(:,t+1) = LL_E * a(:,t) + bb, i.e. the primal problem
308 rhs = LL_E * a(:,t-t_ind_start) + bb;
309 if (LL_I_is_speye == 1)
310 a(:,t+1-t_ind_start) = rhs;
311 elseif (LL_I_is_speye == -1) && isequal(LL_I, speye(size(LL_I)))
312 a(:,t+1-t_ind_start) = rhs;
314 [a(:,t+1-t_ind_start), flag] = bicgstab(LL_I,rhs,[],1000);
316 disp(['error in system solution, solver return flag = ', ...
322 a(:,t+1-t_ind_start)' * K_II * a(:,t+1-t_ind_start) ...
323 -2* a(:,t+1-t_ind_start)' * K_IE * a(:,t-t_ind_start) ...
324 + a(:,t-t_ind_start)' * K_EE * a(:,t-t_ind_start) ...
326 - 2* a(:,t+1-t_ind_start)' * m_I ...
327 + 2* a(:,t-t_ind_start)' * m_E;
328 if model.want_improved_output && model.want_dual == 0 % sum up the correction term for the improved output
329 scalarproduct = (- a(:,t+1-t_ind_start)' * LL_I_correct * dual_data.a(:,t+1-t_ind_start)...
330 + a(:,t-t_ind_start)' * LL_E_correct * dual_data.a(:,t+1-t_ind_start)...
331 + bb_correct * dual_data.a(:,t+1-t_ind_start));
332 correction_term = correction_term + scalarproduct;
336 % ensure real estimators (could be complex due to numerics)
338 res_norm = sqrt(res_norm_sqr);
342 % either use the estimator measured in l2 norm or energy norm. In the
343 % energy case just square the residuals - the rest will be done later.
344 if isequal(model.error_norm, 'l2')
345 Delta(t+1-t_ind_start) = model.L_I_inv_norm_bound * (res_norm + ...
346 model.L_E_norm_bound * Delta(t-t_ind_start));
347 elseif isequal(model.error_norm, 'energy') || (model.want_improved_output == 1 && model.want_dual == 0)
348 Delta(t+1-t_ind_start) = res_norm^2;
352 % Rest of energy error estimator
353 if isequal(model.error_norm, 'energy') || (model.want_improved_output == 1 && model.want_dual == 0)
354 % constant is a lower bound to alpha or beta (depending on your model)
355 % and can either be an empiric constant set in model.constant_LB or it
360 constant = model.constant_LB;
362 Delta = 1/constant*sqrt(sum(Delta)); % sum up the Deltas and take the sqrt
365 % return simulation result
366 simulation_data = [];
367 simulation_data.a = a;
368 simulation_data.Delta = Delta;
371 % in the primal case you can choose to have a reduced output
372 if isfield(model,'name_output_functional') && model.want_dual == 0
373 % ignore this 'theta'-case if you do not use the european option
374 % pricing model. See also the special note at the beginning of this function.
375 if strcmp(model.name_output_functional, 'theta') == 1
376 % the partial time-derivative in the theta-functional is approximated by \frac(P_N^k+1 - P_N^k}{model.deltaT}
377 % so we need a Matrix with the difference of the coefficients to form
378 % u_N^k+1 - u_N^k = RB * Matrix.
379 % This Matrix has the form:
381 % a_1^1-a_1^0 a_1^2-a_1^1 ... a_1^M-a_1^M-1
382 % a_2^1-a_2^0 ... a_2^M-a_2^M-1
386 % a_N^1-a_N^0 a_N^2-a_N^1 ... a_N^M-a_N^M-1
388 %
for t^k und k = 0, ... ,M and N being the number of basis vectors.
390 dP_dt = zeros(reduced_data.N , model.nt);
391 for i = 1 : size(dP_dt,2)
392 dP_dt(:,i) = a(:,i+1) - a(:,i);
394 if model.want_improved_output
395 % use both improved estimators (dual and primal) and multiply with
396 % the stability constant, since both estimators already include a
398 s = reduced_data.s_RB * dP_dt(:,end) - correction_term; % in the improved output the functional is only evaluated at the last timestep
399 Delta_s = Delta*dual_data.Delta*constant;
401 s = reduced_data.s_RB * dP_dt;
402 if isequal(model.error_norm,'l2')
403 % normal sequence of l2-measured estimators
404 Delta_dt = zeros(1, model.nt);
405 for i = 1 : size(dP_dt,2)
406 Delta_dt(i) = simulation_data.Delta(1,i+1) + simulation_data.Delta(1,i); % this builds the standard theta error estimator, which was a sum of 2 l2-error estimators
408 Delta_s = reduced_data.s_l2norm .* Delta_dt; % continuousity constant * RB-Koeffs
409 elseif isequal(model.error_norm,'energy')
410 % estimator measured in energy norm
411 Delta_s = 2*reduced_data.s_l2norm * Delta; % 2*continuousity constant of the functional in theta case
414 simulation_data.s = s(:);
415 simulation_data.Delta_s = Delta_s(:);
417 if model.want_improved_output
418 % use both improved estimators (dual and primal) and multiply with
419 % the stability constant, since both estimators already include a
421 s = reduced_data.s_RB * a(:,end) - correction_term; % in the improved output the functional is only evaluated at the last timestep
422 Delta_s = Delta.*dual_data.Delta*constant;
424 s = reduced_data.s_RB * a;
425 if isequal(model.error_norm,'l2')
426 % normal sequence of l2-measured estimators
427 Delta_s = reduced_data.s_l2norm .* Delta; % continuousity constant * RB-Koeffs
428 elseif isequal(model.error_norm,'energy')
429 % estimator measured in energy norm
430 Delta_s = reduced_data.s_l2norm * Delta; % remember: Delta is different now
433 simulation_data.s = s(:);
434 simulation_data.Delta_s = Delta_s(:);
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
function reduced_data = lin_evol_gen_reduced_data_primal_dual(model, detailed_data)
reduced_data = lin_evol_gen_reduced_data_primal_dual(model, detailed_data)
function reduced_data = lin_evol_primal_dual_gen_reduced_data(model, detailed_data)
reduced_data = lin_evol_primal_dual_gen_reduced_data(model, detailed_data)
function [ constant_LB , constant_UB ] = scm_lower_bound(model, reduced_data)
[constant_LB, constant_UB] = scm_lower_bound(model, reduced_data)
function a0 = rb_init_values_primal_dual(model, detailed_data)
a0 = rb_init_values_primal_dual(model, detailed_data)
function simulation_data = lin_evol_primal_dual_rb_simulation(model, reduced_data)
simulation_data = lin_evol_primal_dual_rb_simulation(model, reduced_data)
function scm_offline_data = scm_offline(model, detailed_data, M_train, D_train)
scm_offline_data = scm_offline(model, detailed_data, M_train, D_train)
function simulation_data = lin_evol_rb_simulation_primal_dual(model, reduced_data)
simulation_data = lin_evol_rb_simulation_primal_dual(model, reduced_data)
function detailed_data = lin_evol_primal_dual_gen_detailed_data(model, model_data)
detailed_data = lin_evol_primal_dual_gen_detailed_data(model, model_data)