rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
lin_ds_rb_simulation_euler_forward.m
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)
3 %
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
7 
8 % Bernard Haasdonk 2.4.2009
9 
10 if model.verbose >= 9
11  disp('entered reduced simulation ');
12 end;
13 
14 % possibly shrink reduced_data detected by this routine
15 reduced_data_t = lin_ds_reduced_data_subset(model,reduced_data);
16 
17 %model.affine_decomp_mode = 'complete';
18 model.decomp_mode = 2; % coefficients
19 
20 model = model.set_time(model,0);
21 %model.t = 0;
22 
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;
27 else
28  save_time_index = ones(1,model.nt+1);
29  save_time_indices = 1:(model.nt+1);
30 end;
31 
32 x0coeff = model.x0_function_ptr(model,[]);
33 xr = lincomb_sequence(reduced_data_t.x0r,x0coeff);
34 
35 
36 % instead of large dense matrix allocation equivalent:
37 % a' Mtemp b =
38 % a' G b - (a' W) (V' G b) - (a' G V) (W' b) + (a' W) (V' G V) (W' b)
39 %
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:
44 %
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
49 
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;
57 end;
58 
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));
62 
63 if model.error_estimation
64  DeltaX = zeros(1,length(save_time_indices));
65 end;
66 
67 time_indices = zeros(1,length(save_time_indices));
68 time = zeros(1,length(save_time_indices));
69 
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
73 
74 if save_time_index(t_ind+1)
75  Xr(:,t_column) = xr(:);
76  time_indices(t_column) = t_ind;
77  time(t_column) = t;
78  %u = model.u_function_ptr(model);
79 
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;
84 
85  if model.error_estimation
86  DeltaXt = model.state_bound_constant_C1*sqrt(max(e0sqr,0));
87  DeltaX(t_column) = DeltaXt;
88  end;
89  t_column = t_column + 1;
90 end;
91 
92 model.dt = model.T/model.nt;
93 
94 operators_required = 1;
95 if ~isfield(model,'data_const_in_time')
96  model.data_const_in_time = 0;
97 end,
98 
99 % loop over time
100 for t_ind = 1:model.nt
101 % keyboard;
102  if model.verbose >= 10
103  disp(['entered time-loop step ',num2str(t_ind)]);
104  end;
105  if model.verbose == 9
106  fprintf('.');
107  end;
108 
109  t = (t_ind-1)*model.dt;
110  model = model.set_time(model,t);
111 % model.t = t;
112 
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);
118 
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!!
127  end;
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);
134 % model.t = t;
135  if model.data_const_in_time
136  operators_required = 0;
137  end;
138  end;
139 
140  u = model.u_function_ptr(model);
141  ddt_xr = Ar * xr + Br * u;
142  xr_old = xr;
143  xr = xr + model.dt * ddt_xr;
144 
145  if model.error_estimation
146 
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;
152 
153  % error estimator at next time
154  DeltaXt = DeltaXt + ...
155  model.state_bound_constant_C1 * model.dt ...
156  * sqrt(max(Rnormsqr,0));
157  end;
158 
159  % store result
160  if save_time_index(t_ind+1)
161  Xr(:,t_column) = xr;
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);
168 % model.t = t;
169  Y(:,t_column) = Cr * xr + D * u;
170  time(t_column) = t;
171 
172  if model.error_estimation
173  DeltaX(t_column) = DeltaXt;
174  end;
175 
176  t_column = t_column + 1;
177  end;
178 
179 end;
180 
181 if model.verbose == 9
182  fprintf('\n');
183 end;
184 
185 rb_sim_data = [];
186 rb_sim_data.Xr = Xr;
187 rb_sim_data.Y = Y;
188 if model.error_estimation
189  rb_sim_data.DeltaX = DeltaX;
190  rb_sim_data.DeltaY = DeltaX * model.output_bound_constant_C2;
191 end;
192 rb_sim_data.time = time;
193 rb_sim_data.time_indices = time_indices;
194 
195 %if model.debug
196 % disp('short before leavin lin_ds_rb_simulation');
197 %end;
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
Definition: verbose.m:17