1 function rb_sim_data = lin_ds_rb_simulation_euler_forward(model,reduced_data)
2 %
function rb_sim_data = lin_ds_rb_simulation_euler_forward(model,reduced_data)
4 % simulation of reduced linear system with
explicit Euler forward
5 % discretization. This method is deprecated, as same result can be
6 % obtained with theta-scheme, theta=0 in lin_ds_rb_simulation
8 % Bernard Haasdonk 2.4.2009
11 disp(
'entered reduced simulation ');
14 % possibly shrink reduced_data detected by
this routine
15 reduced_data_t = lin_ds_reduced_data_subset(model,reduced_data);
17 %model.affine_decomp_mode =
'complete';
18 model.decomp_mode = 2; % coefficients
20 model = model.set_time(model,0);
23 if isfield(model,
'save_time_indices')
24 save_time_index = zeros(1,model.nt+1);
25 save_time_index(model.save_time_indices(:)+1) = 1;
26 save_time_indices = find(save_time_index==1)-1;
28 save_time_index = ones(1,model.nt+1);
29 save_time_indices = 1:(model.nt+1);
32 x0coeff = model.x0_function_ptr(model,[]);
33 xr = lincomb_sequence(reduced_data_t.x0r,x0coeff);
36 % instead of large dense matrix allocation equivalent:
38 % a' G b - (a' W) (V' G b) - (a' G V) (W' b) + (a' W) (V' G V) (W' b)
40 % these components can not be summed up here,
41 % as we must be able to extract online data for other reduced dimensions!!!
42 % hence no longer computation of m here
43 % but decomposition as follows:
45 % m00: 2d cell array of values x0{q}
'G x0{q'}
46 % VtGx0: 1d cell array of of vectors (V
' G x0{q})
47 % Wtx0: 1d cell array of of vectors (W' x0{q})
48 % VtGV : matrix V
' G V
50 if model.error_estimation
51 %e0sqr = lincomb_sequence2(reduced_data_t.m,x0coeff,x0coeff);
52 x0tGx0 = lincomb_sequence2(reduced_data_t.m00,x0coeff,x0coeff);
53 Wtx0 = lincomb_sequence(reduced_data_t.Wtx0,x0coeff);
54 VtGx0 = lincomb_sequence(reduced_data_t.VtGx0,x0coeff);
55 VtGV = reduced_data_t.VtGV;
56 e0sqr = x0tGx0 - (Wtx0)'*VtGx0 - (VtGx0)
' * Wtx0 + (Wtx0)' * VtGV * Wtx0;
59 D = model.D_function_ptr(model,[]);
60 Xr = zeros(length(xr(:)),length(save_time_indices));
61 Y = zeros(size(reduced_data_t.Cr{1},1),size(Xr,2));
63 if model.error_estimation
64 DeltaX = zeros(1,length(save_time_indices));
67 time_indices = zeros(1,length(save_time_indices));
68 time = zeros(1,length(save_time_indices));
70 t_column = 1; % next column index to be filled in output
71 t_ind = 0; % current t_ind between 0 and nt
72 t = 0; % absolute time between 0 and T
74 if save_time_index(t_ind+1)
75 Xr(:,t_column) = xr(:);
76 time_indices(t_column) = t_ind;
78 %u = model.u_function_ptr(model);
80 Ccoeff = model.C_function_ptr(model,[]);
81 Cr = lincomb_sequence(reduced_data_t.Cr,Ccoeff);
82 u = model.u_function_ptr(model);
83 Y(:,t_column) = Cr * xr + D * u;
85 if model.error_estimation
86 DeltaXt = model.state_bound_constant_C1*sqrt(max(e0sqr,0));
87 DeltaX(t_column) = DeltaXt;
89 t_column = t_column + 1;
92 model.dt = model.T/model.nt;
94 operators_required = 1;
95 if ~isfield(model,'data_const_in_time')
96 model.data_const_in_time = 0;
100 for t_ind = 1:model.nt
103 disp(['entered time-loop step ',num2str(t_ind)]);
109 t = (t_ind-1)*model.dt;
110 model = model.set_time(model,t);
113 if operators_required
114 Acoeff = model.A_function_ptr(model,[]);
115 Ar = lincomb_sequence(reduced_data_t.Ar,Acoeff);
116 Bcoeff = model.B_function_ptr(model,[]);
117 Br = lincomb_sequence(reduced_data_t.Br,Bcoeff);
119 if model.error_estimation
120 M1 = lincomb_sequence2(reduced_data_t.M1,Acoeff,Acoeff);
121 M2 = lincomb_sequence2(reduced_data_t.M2,Bcoeff,Bcoeff);
122 M3 = reduced_data_t.M3;
123 M4 = lincomb_sequence2(reduced_data_t.M4,Bcoeff,Acoeff);
124 M5 = lincomb_sequence(reduced_data_t.M5,Acoeff);
125 M6 = lincomb_sequence(reduced_data_t.M6,Bcoeff);
126 % output matrices required at next time!!
128 model = model.set_time(model,t+model.dt);
129 % model.t = t+model.dt;
130 Ccoeff = model.C_function_ptr(model,[]);
131 Cr = lincomb_sequence(reduced_data_t.Cr,Ccoeff);
132 D = model.D_function_ptr(model,[]);
133 model = model.set_time(model,t);
135 if model.data_const_in_time
136 operators_required = 0;
140 u = model.u_function_ptr(model);
141 ddt_xr = Ar * xr + Br * u;
143 xr = xr + model.dt * ddt_xr;
145 if model.error_estimation
147 % residual at previous timestep => "rectangular" integration
148 Rnormsqr = xr_old' * M1 * xr_old + u' * M2 * u + ...
149 ddt_xr' * M3 * ddt_xr + ...
150 2*u' * M4 * xr_old - 2*ddt_xr' * M5 * xr_old - ...
151 2 * ddt_xr' * M6 * u;
153 % error estimator at next time
154 DeltaXt = DeltaXt + ...
155 model.state_bound_constant_C1 * model.dt ...
156 * sqrt(max(Rnormsqr,0));
160 if save_time_index(t_ind+1)
162 time_indices(t_column) = t_ind;
163 % for output u required at next timestep!
164 model = model.set_time(model,t+model.dt);
165 % model.t = t+model.dt;
166 u = model.u_function_ptr(model);
167 model = model.set_time(model,t);
169 Y(:,t_column) = Cr * xr + D * u;
172 if model.error_estimation
173 DeltaX(t_column) = DeltaXt;
176 t_column = t_column + 1;
188 if model.error_estimation
189 rb_sim_data.DeltaX = DeltaX;
190 rb_sim_data.DeltaY = DeltaX * model.output_bound_constant_C2;
192 rb_sim_data.time = time;
193 rb_sim_data.time_indices = time_indices;
196 % disp('
short before leavin lin_ds_rb_simulation');
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...