3 %
function demonstrating RB approach
for linear dynamical systems
4 % according to MATHMOD 2009 paper
6 % Method also demonstrates
new call of all reduced schemes:
8 % model_data = gen_model_data(model);
10 % sim_data = detailed_simulation(model,model_data);
11 % produces some structure with suitable fields, e.g DOF vector
12 % output sequence, reference to output files, simulation
13 % information, etc. Last parameter is model-
struct, possibly
14 % precedent parameters, e.g. grid, etc.
15 % so not only DOF-vector is the
return argument, but a structure
16 % detailed_data = gen_detailed_data(model,model_data);
17 % produce high-dimensional data, e.g. reduced basis, etc.
19 % reduced_data = gen_reduced_data(model,detailed_data);
20 % produce low-dimensional data, e.g. gram matrices between
21 % reduced basis vectors, subgrid extraction, etc.
22 % mu is not known (so is the previous offline-data or
23 % combination of offline and online data)
24 % rb_sim_data = rb_simulation(model,reduced_data)
25 % perform reduced simulation, i.e. mu known
26 % produces some structure with suitable fields, e.g DOF vector,
29 % Further routines
"extending" datasets, e.g.:
30 % rb_sim_data = rb_reconstruction(rb_sim_data,model)
31 % reduced_data = extract_reduced_data_subset(reduced_data,model)
34 % Bernard Haasdonk 2.4.2009
36 %step = 1 % simulation and visualization of dynamical system
37 %step = 2 % reduced simulation and visualization of approximation
38 step =3 % detailed and reduced simulation and visualization of
42 model = lin_ds_model_default;
44 model.data_const_in_time = 1;
47 model.affinely_decomposed = 1;
48 model.A_function_ptr = @A_func;
49 % ratio of first to second ellipses axes
50 %model.A_axes_ratio = 2;
51 model.A_axes_ratio = 2;
52 model.B_function_ptr = @B_func;
53 model.C_function_ptr = @C_func;
54 model.D_function_ptr = @D_func;
55 model.x0_function_ptr = @x0_func;
56 % shift of x0 in third coordinate
58 model.u_function_ptr = @(model) 0.1; % define u constant
59 %model.u_function_ptr = @(model) sin(4*model.t);
60 %model.u_function_ptr = @(model) 0; % define u constant 0
61 %model.u_function_ptr = @(model) model.t; % define u = 1 * t
62 model.G_matrix_function_ptr = @(model,model_data) eye(3,3);
63 % the following lead to both larger C1 and C2!!!
64 %model.G_matrix_function_ptr = @(model) diag([2,0.5,1]);
65 %model.G_matrix_function_ptr = @(model) diag([0.5,2,1]);
67 % dummy parametrization:, parameter is not used...
68 model.mu_names = {
'A_axes_ratio',
'x0_shift'};
69 %model.mu_ranges = {[0.5,2],[0,1]};
70 model.mu_ranges = {[0.5,0.5],[1,1]};
71 model.affinely_decomposed = 1;
73 % set
function pointers to main model methods
74 %model.gen_model_data = @lin_ds_gen_model_data;
75 %model.detailed_simulation = @lin_ds_detailed_simulation;
76 %model.gen_detailed_data = @lin_ds_gen_detailed_data;
77 %model.gen_reduced_data = @lin_ds_gen_reduced_data;
78 %model.rb_simulation = @lin_ds_rb_simulation;
79 %model.plot_sim_data = @lin_ds_plot_sim_data;
80 %model.rb_reconstruction = @lin_ds_rb_reconstruction;
81 % determin model data, i.e. matrix components, matrix G
83 model_data = gen_model_data(model); % matrix components
85 % be sure to determine correct constants once,
if anything in the
86 % problem setting has changed!!
88 % the following only needs to be performed,
if the model are changed.
90 disp(
'Estimating continuity bound constants')
91 model.estimate_lin_ds_nmu = 10;
92 model.estimate_lin_ds_nX = 100;
93 model.estimate_lin_ds_nt = 100;
95 [C1,C2] = lin_ds_estimate_bound_constants(model,model_data)
96 model.state_bound_constant_C1 = C1;
97 model.output_bound_constant_C2 = C2;
99 model.state_bound_constant_C1 = 2;
100 model.output_bound_constant_C2 = sqrt(3);
102 model.error_estimation = 1;
103 model_data.RB = eye(3,2);
104 model.RB_generation_mode = 'from_detailed_data';
107 case 1 % simulation and visualization
108 % later: sim_data = detailed_simulation(model);
110 sim_data = detailed_simulation(model,model_data);
114 case 2 % reduced simulation
115 %define reduced basis of first two coordinates (i.e. error only
120 % produce high-dimensional data, e.g. reduced basis, etc.
122 detailed_data = gen_detailed_data(model,model_data);
124 % produce low-dimensional data, e.g. gram matrices between
125 % reduced basis vectors, subgrid extraction, etc.
126 % mu is not known (should comprise previous offline and online steps)
128 reduced_data = gen_reduced_data(model,detailed_data);
129 model.N = reduced_data.N;
131 % rb_sim_data = rb_simulation(reduced_data,model)
132 % perform reduced simulation, i.e. mu known
133 % produces some structure with suitable fields, e.g DOF vector
135 rb_sim_data = rb_simulation(model,reduced_data);
137 rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
139 model.plot_title = 'reduced trajectory';
140 %model.plot_indices = 1:3; => default
144 case 3 % detailed and reduced simulation
146 %model.u_function_ptr = @(model) 0; % define u constant 0
147 % => leads to error 0 :-) if x0_shift = 0
149 %model.u_function_ptr = @(model) 0.1; % define u constant
150 %=> error linearly increasing
151 model.u_function_ptr = @(model) 0.3*sin(2*model.t)+0.002*model.t;
154 sim_data = detailed_simulation(model,model_data);
155 detailed_data = gen_detailed_data(model,model_data);
156 model.enable_error_estimator=1;
157 reduced_data = gen_reduced_data(model,detailed_data);
158 model.N = reduced_data.N;
159 rb_sim_data = rb_simulation(model,reduced_data);
160 rb_sim_data = rb_reconstruction(model,detailed_data, ...
163 figure, plot3(sim_data.X(1,:),...
169 plot3(rb_sim_data.X(1,:),...
170 rb_sim_data.X(2,:),...
171 rb_sim_data.X(3,:),...
174 title('exact and reduced trajectories');
175 legend({
'exact',
'reduced'});
178 error.X = sim_data.X-rb_sim_data.X;
179 error.Y = sim_data.Y-rb_sim_data.Y;
181 err = sqrt(sum((detailed_data.G * error.X).*error.X,1));
182 figure, plot([rb_sim_data.DeltaX(:),err(:)]),
183 Ylim =
get(gca,
'Ylim');
184 set(gca,
'Ylim',[0,max(Ylim)]);
185 title(
'rb state error and estimator');
186 legend({
'estimator',
'error'});
188 figure, plot([rb_sim_data.DeltaY(:),abs(error.Y(:))]);
189 title(
'rb output error and estimator');
190 Ylim =
get(gca,
'Ylim');
191 set(gca,
'Ylim',[0,max(Ylim)]);
192 legend({
'estimator',
'error'});
197 error(
'step unknown');
200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
201 % functions defining the dynamical system and input/initial data
202 function A = A_func(model,model_data);
203 A = eval_affine_decomp(@A_components,@A_coefficients,...
206 % A_parameter: ratio of first to second axis of
207 function Acomp = A_components(model,model_data)
208 %Acomp = {[0 -0.25 0; 4 0 0 ; 0 0 0 ]};
209 %Acomp = {[0 -1 0; 1 0 0 ; 0 0 0 ]};
210 %Acomp = {[-0.01 0 0; 0 -0.9 0 ; 0 0 -0.001 ]};
211 Acomp = {[0 1 0; 0 0 0 ; 0 0 0 ],...
212 [0 0 0; 1 0 0 ; 0 0 0 ] ...
215 function Acoeff = A_coefficients(model)
216 Acoeff = [-model.A_axes_ratio, 1/model.A_axes_ratio];
218 function B = B_func(model,model_data)
219 B = eval_affine_decomp(@B_components,@B_coefficients,model,model_data);
221 function Bcomp = B_components(model,model_data)
224 function Bcoeff = B_coefficients(model)
227 function C = C_func(model,model_data)
228 C = eval_affine_decomp(@C_components,@C_coefficients,model,model_data);
230 function Ccomp = C_components(model,model_data)
232 %Ccomp = {[2,0,1],[0,2,1]};
234 function Ccoeff = C_coefficients(model)
236 %Ccoeff = [model.C_weight,1-model.C_weight];
238 function D = D_func(model,model_data)
239 D = 0; % no decomposition required for scheme
241 function x0 = x0_func(model,model_data)
242 x0 = eval_affine_decomp(@x0_components,@x0_coefficients,model,model_data);
244 function x0comp = x0_components(model,model_data)
245 x0comp = {[1;0;0],[1;0;1]};
247 function x0coeff = x0_coefficients(model)
248 x0coeff = [1-model.x0_shift,model.x0_shift];
function p = plot_sim_data(model, model_data, sim_data, plot_params)
function performing the plot of the simulation results as specified in model.
function detailed_data = rb_basis_generation(model, detailed_data)
reduced basis construction with different methods
function demo_lin_ds()
function demonstrating RB approach for linear dynamical systems according to MATHMOD 2009 paper ...