rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
cpde2016.m
Go to the documentation of this file.
1 function cpde2016(input)
2 %CPDE2016 This function reproduces all the results from the CPDE paper
3 %
4 % Andreas Schmidt, 2016
5 
6 % Function for building all the results in the paper
7 if ~iscell(input)
8  input = {input};
9 end
10 
11 % Reset the random number generator to get the same results
12 rng default;
13 
14 n0 = 20;
15 n = n0^2;
16 N = 2*n0^2;
17 
18 % We consider a (damped) wave equation:
19 diffA = my_2d_matrix(n0, '0', '0', '0', '0.1');
20 A = @(mu) [zeros(n) eye(n); diffA -mu(1)*eye(n)];
21 B0 = [zeros(n,1); my_2d_vector(n0, '(0.3 <= x) & (x <= 0.7) & (0.4 <= y) & (y <= 0.6)')];
22 B = @(mu) B0/norm(B0);
23 Bw0 = [[zeros(n,1); my_2d_vector(n0, '(sin(pi*x)*sin(pi*y))')] zeros(N,1)];
24 Bw = @(mu) 20*Bw0/norm(Bw0);
25 
26 Dz = @(mu) 0.1*[0; 1];
27 C0 = [my_2d_vector(n0, '(0.2 <= x) & (x <= 0.8) & (0.3 <= y) & (y <= 0.4)')' zeros(1,n)];
28 C = @(mu) C0/norm(C0);
29 Cz0 = [C0; zeros(1,N)];
30 Cz = @(mu) mu(2)*[Cz0/norm(Cz0)];
31 Dw = @(mu) 0.05*[0 1];
32 
33 % Set up the RBmatlab models for the two algebraic Riccati equations:
34 descrP.A_comp = { [zeros(n) eye(n); diffA, zeros(n)], ...
35  blkdiag(zeros(n), -eye(n)) };
36 descrP.A_coeff = @(m)[ 1, m.mu(1) ];
37 descrP.B_comp = { B0/norm(B0) };
38 descrP.B_coeff = @(m) 1;
39 descrP.C_comp = { Cz0/norm(Cz0) };
40 descrP.C_coeff = @(m) m.mu(2);
41 descrP.R_comp = { (0.1*[0; 1])'*(0.1*[0; 1]) };
42 descrP.mu = [1 1];
43 descrP.mu_ranges = {[0.1, 2], [1, 100]};
44 descrP.mu_names = {'mu_damping', 'mu_Q'};
45 descrP.n = 2*n;
46 model_P = ARE.DefaultModel(descrP);
47 model_data_P = gen_model_data(model_P);
48 
49 descrQ = descrP;
50 descrQ.A_comp = {descrP.A_comp{1}', descrP.A_comp{2}'};
51 descrQ.B_comp = {C0'/norm(C0)};
52 descrQ.B_coeff = descrP.C_coeff;
53 descrQ.C_comp = { 20*Bw0'/norm(Bw0) };
54 descrQ.C_coeff = @(m) 1;
55 descrQ.R_comp = { (0.05*[0 1])*(0.05*[0 1])' };
56 model_Q = ARE.DefaultModel(descrQ);
57 model_data_Q = gen_model_data(model_Q);
58 
59 % Define the training parameters:
60 N_train = 40;
61 mu_train = gen_mus(N_train);
62 
63 % Make sure we have a valid folder for the results:
64 [~,~,~] = mkdir('result');
65 [~,~,~] = mkdir('plots');
66 
67 % For the simulation we still need some more data:
68 t = linspace(0,5,500)'; % The time span, final time T = 5
69 u = sin(t); % Just some input function
70 w = randn(length(t), 2); % Disturbance
71 x0 = [my_2d_vector(n0, '(sin(pi*x)*sin(pi*y))^2'); zeros(n,1)];
72 x0 = x0/norm(x0);
73 
74 mu_test = [1 1];
75 
76 for sk = 1:length(input)
77 switch input{sk}
78  case 'all-results'
79  display('Calculating all results')
80  paper({'basis-Q', 'basis-P', 'basis-xhat', 'basis-plot', 'approaching-z2-plot'})
81 
82  case 'basis-plot'
83  display('Plotting results for the basis generation plot')
84 
85  case 'full-dimensional'
86  display('Performing a full dimensional simulation')
87 
88  % Assemble the full system:
89  sys = get_ss(mu_test);
90  % Perform uncontrolled simulation:
91  [sim_uc,~,x_uc] = lsim(sys, [0*u, w], t, [x0]);
92 
93  % For comparison controlled simulation:
94  closed_loop = get_full_c_sys(mu_test);
95  [sim_c,~,x_c] = lsim(closed_loop, [0*u, w], t, [x0; 0*x0]);
96 
97  plot(t, [sim_uc(:,1) sim_c(:,1)])
98  title('Output z_1')
99  legend('Uncontrolled', 'Controlled')
100 
101  case 'full-solution-3d'
102  % Simulate the system with the noisy input and without any control.
103  % Output a surface plot
104  sys = get_ss([1 1 1]);
105  [~,~,x] = lsim(sys, [0*u, w], t, x0);
106 
107  for k = 1:length(t)
108  % Use the first n elements in the vector (they represent the
109  % system states)
110  X = x(k,1:n)';
111  surf(linspace(0,1,n0),linspace(0,1,n0),reshape(X,n0,n0));
112  axis([0 1 0 1 min(min(x)) max(max(x))])
113  pause
114  end
115 
116  case 'uncontrolled_output_plot'
117  sys = get_ss([1 10 1]);
118  sim1 = lsim(sys, [0*u, w], t, [x0]);
119  sys2 = get_full_c_sys([1 10 1]);
120  sim2 = lsim(sys2, [0*u, w], t, [x0; 0*x0]);
121  clf;hold on;
122  p11 = plot(t, [sim1(:,3)], 'LineWidth', 0.8, 'Color', [0.65 0.65 0.65]);
123  plot(t, [sim2(:,3)], 'LineWidth', 0.8, 'Color', [0.65 0.65 0.65]);
124  p1 = plot(t, [sim1(:,4)], 'Color', 'b', 'LineWidth', 1.2);
125  p2 = plot(t, [sim2(:,4)], 'Color', 'r', 'LineWidth', 1.2);
126  legend([p1, p2], 'Uncontrolled', 'Controlled')
127  xlabel('Time t')
128  ylabel('Output')
129  matlab2tikz( 'plots/output-example.tex', 'height', '\fheight', 'width', '\fwidth', 'parseStrings', false, 'extraAxisOptions','ylabel shift={-10pt}' )
130 
131  case 'svd_plot'
132  dsxhat = diag(Sxhat);
133  dsp = diag(SP);
134  dsq = diag(SQ);
135  r = 1:80;
136  clf;
137  p = semilogy(r, dsp(r), '-', r, dsq(r), '--', r, dsxhat(r), '-.');
138  p(1).LineWidth = [1.5];
139  p(2).LineWidth = [1.5];
140  p(3).LineWidth = [1.5];
141  hold on;
142  legend('$\bar P$', '$\bar Q$', '$\bar x$')
143 
144  ylim = get(gca, 'YLim')
145  plot(rp, dsp(rp), 'kd', 'MarkerFaceColor', 'k')
146  line([rp rp], ylim, 'Color', 'k')
147 
148  plot(rxhat, dsxhat(rxhat), 'kd', 'MarkerFaceColor', 'k')
149  line([rxhat rxhat], ylim, 'Color', 'k')
150 
151  plot(rq, dsq(rq), 'kd', 'MarkerFaceColor', 'k')
152  line([rq rq], ylim, 'Color', 'k')
153 
154  matlab2tikz( 'plots/svd-decay.tex', 'height', '\fheight', 'width', '\fwidth', 'parseStrings', false )
155 
156  case 'reduced_full_simulation'
157  for k = 1:10
158  mu_test = gen_mus(1);
159 
160  full_closed_loop = get_full_c_sys(mu_test);
161  reduced_closed_loop = get_reduced_c_sys(mu_test, rp, rq, rxhat);
162  [full_sim_c,~,full_x_c] = lsim(full_closed_loop, [0*u, w], t, [x0; 0*x0]);
163  [reduced_sim_c,~,x_c] = lsim(reduced_closed_loop, [0*u, w], t, [x0; Vxhat(:,1:rxhat)'*0*x0]);
164  max(real(eig(reduced_closed_loop.a)))
165  plot(t, [full_sim_c(:,1), reduced_sim_c(:,1)])
166  pause
167  end
168 
169  case 'approaching-z2-plot'
170  this_t = linspace(0,1);
171  this_w = randn(100, 2);
172  this_u = randn(100,1);
173 
174  %rxhat = size(Vxhat,2);
175  full_closed_loop = get_full_c_sys(mu_test);
176  this_u = 0*exp(2*this_t').^2;
177  [full_sim_c,~,full_x_c] = lsim(full_closed_loop, [this_u, this_w], this_t, [x0; 0*x0]);
178 
179  clf;
180  plot(this_t, full_sim_c(:,2), 'linewidth', 1.2, 'DisplayName', 'Full');
181  hold all
182  sizes = [5 10 size(VP,2)];
183  ls = {':', '--', '-.'};
184 
185  for r = 1:length(sizes)
186  reduced_closed_loop = get_reduced_c_sys(mu_test, sizes(r), size(VQ,2), rxhat);
187  reduced_sim_c = lsim(reduced_closed_loop, [this_u, this_w], this_t, [x0; Vxhat(:,1:rxhat)'*0*x0]);
188  plot(this_t, reduced_sim_c(:,2), 'LineStyle', ls{r}, 'LineWidth', 1);
189  end
190  legend({'Full', 'N_P=5', 'N_P=10', 'N_P=60'}, 'Location', 'southeast')
191  xlabel('Time t')
192  ylabel('Control u(t)')
193  matlab2tikz( 'plots/approaching_z2.tex', 'height', '\fheight', 'width', '\fwidth' )
194 
195 
196  case 'timings'
197  display('Performing some performance checks')
198  tsum = 0;
199  lmax = 0;
200  for l = 1:0
201  %mu = gen_mus(1);
202  a = A(mu);
203  b = B(mu);
204  c = Cz(mu);
205  d = Dz(mu);
206  f = tic;
207  p = care(full(a),b,c'*c,d'*d);
208  tsum = tsum + toc(f);
209  end
210  display(['Average time for solution of one ARE is ', num2str(tsum/lmax)])
211 
212  mut = gen_mus(1);
213  [q,L] = Q(mut);
214  [p,K] = P(mut);
215 
216  % Simulate the dynamical system
217  tmp_obs = ss(A(mut) - L*C(mut) - B(mut)*K, L, C(mut), []);
218  this_t = linspace(0,10);
219  this_y = sin(pi*this_t);
220  f = tic;
221  [asd,~,~] = lsim(tmp_obs, this_y, this_t);
222  toc(f)
223 
224  V = Vxhat(:, 1:rxhat);
225  tmp_obs = ss(V'*(A(mut) - L*C(mut) - B(mut)*K)*V, V'*L, C(mut)*V, []);
226  f = tic;
227  [asd,~,~] = lsim(tmp_obs, this_y, this_t);
228  toc(f)
229 
230  case 'error-surface'
231  N_test = 1;
232  mu_test = gen_mus(N_test);
233  full_sys = {};
234  for k = 1:N_test
235  full_sys{k} = get_full_c_sys(mu_test(k,:));
236  end
237 
238  NP = [5 10 20 40 60];
239  NX = [5 10 15 17];
240  data = zeros(length(NX),length(NP));
241  grcs = @(m,np,nq,nx) get_reduced_c_sys(m,np,nq,nx);
242 
243  for i = 1:length(NP)
244  parfor j = 1:length(NX)
245  [i, j]
246  r = [];
247  for k = 1:N_test
248  rsys = grcs(mu_test(k,:),NP(i),2,NX(j));
249  try
250  fs = norm(full_sys{k}(1:2, 2:3),2);
251  r(k) = abs(norm(full_sys{k}(1:2, 2:3),2) - norm(rsys(1:2, 2:3), 2))/fs;
252  catch e
253  end
254  end
255  data(j,i) = max(r)
256  end
257  end
258 
259  save('result/error_surface.mat', 'data');
260 
261  case 'error-surface-output'
262  load('result/error_surface.mat');
263  input = {};
264  input.data = data;
265  NP = [5 10 20 40 60];
266  NX = [5 10 15 17];
267  data
268  input.tableColLabels = arrayfun(@(i) num2str(NP(i)), 1:length(NP), 'UniformOutput', false);
269  input.tableRowLabels = arrayfun(@(i) num2str(NX(i)), 1:length(NX), 'UniformOutput', false);
270  input.dataFormat = {'%.2e'};
271  latexTable(input);
272 
273  % =======================================================================
274  % Functions for the basis generation
275  % =======================================================================
276  case 'basis-P'
277  % Implementation of the LRFG algorithm (not using rb matlab for
278  % portability). However, this should be moved to RBmatlab at a certain
279  % point!
280  if exist('result/basis_P.mat','file')
281  display('Loading P basis from file')
282  load('result/basis_P.mat')
283  continue
284  end
285  dd = gen_detailed_data(model_P, model_data_P);
286  [VP,P_info] = lrfg( @P, @Phat );
287  save('result/basis_P.mat', 'VP', 'P_info');
288 
289  case 'basis-Q'
290  if exist('result/basis_Q.mat','file')
291  display('Loading Q basis from file')
292  load('result/basis_Q.mat')
293  continue
294  end
295  [VQ,Q_info] = lrfg( @Q, @Qhat );
296  save('result/basis_Q.mat', 'VQ', 'Q_info');
297 
298  case 'clear_basis'
299  display('Clearing the files')
300  delete('result/basis_P.mat', 'result/basis_Q.mat', ...
301  'result/basis_xhat.mat')
302 
303  case 'load_basis'
304  display('Loading all data matrices')
305  load('result/basis_Q.mat')
306  load('result/basis_P.mat')
307  load('result/basis_xhat.mat')
308  rp = size(VP,2)
309 
310  s = diag(Sxhat);
311  for k = 1:length(s)
312  eE = sum(s(1:k).^2)/sum(s.^2);
313  if eE >= 0.99999
314  break
315  end
316  end
317  rxhat = k;
318  rxhat
319  rq = size(VQ,2)
320  display(' ... done loading the basis')
321 
322  case 'basis-xhat'
323  display('Performing basis generation for the regulator')
324  if exist('result/basis_xhat.mat')
325  load('result/basis_xhat.mat')
326  continue
327  end
328 
329  snapshots = [];
330  for j = 1:10
331  % Assemble the observer:
332  % sys = get_ss(mu_train(j,:));
333  display(['Iteration ', num2str(j)])
334  mut = mu_train(j,:);
335  [~,L] = Q(mut);
336  [~,K] = P(mut);
337 
338  tmp_obs = ss(A(mut) - L*C(mut) - B(mut)*K, L, C(mut), []);
339  this_t = linspace(0,2*pi,500);
340  this_y = sin(t);
341  [~,~,x] = lsim(tmp_obs, this_y, this_t);
342  snapshots = [snapshots x'];
343  end
344  [Vxhat,Sxhat,~] = svd(snapshots);
345 
346  save('result/basis_xhat.mat', 'Vxhat', 'Sxhat');
347 
348 end
349 end % End of the "for"-loop
350 
351 %% ========================================================================
352 function syscl = get_full_c_sys(mu)
353 % Get the system with full controller for comparison
354 syscl = get_reduced_c_sys(mu, [], [], []);
355 end
356 
357 function syscl = get_reduced_c_sys(mu, red_p, red_q, red_xhat)
358 % Get the reduced systems. If one of the sizes red_p,red_q or red_xhat is
359 % empty, the full matrix/equation will be used.
360 sys = get_ss(mu);
361 
362 % 1. Control gain:
363 if ~isempty(red_p) && red_p > 0
364  [~,khat] = detailedPhat(mu, red_p);
365 else
366  dsim = filecache_function(@detailed_simulation, model_P, model_data_P);
367  khat = -dsim.K;
368  %[~,khat] = P(mu);
369 end
370 
371 % 2. Estimator gain:
372 if ~isempty(red_q) && red_q > 0
373  [~,lhat] = Qhat(mu, red_q);
374 else
375  dsim = filecache_function(@detailed_simulation, model_Q, model_data_Q);
376  lhat = -dsim.K';
377  %[~,lhat] = Q(mu);
378 end
379 
380 % 3. Assemble the controller:
381 if ~isempty(red_xhat) && red_xhat > 0
382  Vx = Vxhat(:, 1:red_xhat);
383 else
384  Vx = eye(N);
385 end
386 
387 % Assemble the controller:
388 controller = ss(Vx'*(A(mu) - lhat*C(mu) - B(mu)*khat)*Vx, Vx'*lhat, -khat*Vx, []);
389 
390 % Now close the loop and return the system:
391 % 1 --> u (Input of system)
392 % 3 --> y (Noisy output of system)
393 % Positive feedback!
394 syscl = feedback(sys, controller, [1], [3], 1);
395 end
396 function p = gen_mus(N)
397 % Generate random parameter vectors (p in [4, 41])
398 p = repmat([1 1], N, 1) + repmat([3, 19], N, 1).*rand(N,2);
399 end
400 function sys = get_ss(mu)
401 % Assemble the system and return it
402 sys = ss(full(A(mu)), [B(mu) Bw(mu)], [Cz(mu); C(mu); C(mu)], [[Dz(mu)] ...
403  zeros(2); 0 Dw(mu); 0 0 0], 'inputname', {'u', ...
404  'w1', 'w2'}, 'outputname', {'z1', 'z2', 'y', 'ytrue',});
405 end
406 
407 function [p,k] = P(mu)
408 % This function produces the "true" feedback gain and the true
409 % solution to the control ARE
410 p = care(A(mu),B(mu),Cz(mu)'*Cz(mu),Dz(mu)'*Dz(mu));
411 k = B(mu)'*inv(Dz(mu)'*Dz(mu))*p;
412 end
413 
414 function [phat,khat,PN,res] = Phat(mu,N,v)
415 % Calculate the reduced solution to a control ARE
416 
417 % Use the global VP matrix or the local one:
418 if exist('v','var')
419  V = v;
420 else
421  V = VP(:, 1:N);
422 end
423 Ar = V'*A(mu)*V;
424 Br = V'*B(mu);
425 Qr = V'*Cz(mu)'*Cz(mu)*V;
426 R = Dz(mu)'* Dz(mu);
427 
428 PN = care(Ar, Br, Qr, R);
429 phat = V*PN*V';
430 khat = B(mu)'*inv(R)*phat;
431 
432 if nargout == 4
433  nPhat = norm(phat,'fro');
434  normA = norm(A(mu),'fro');
435  d = norm(Cz(mu)'*Cz(mu),'fro');
436 
437  res = norm( A(mu)'*phat + phat*A(mu) - phat*B(mu)*inv(R)*B(mu)'*phat + Cz(mu)'*Cz(mu), 'fro') / (2*nPhat*normA + d + nPhat^2*norm(inv(R),'fro')*norm(B(mu),'fro')^2);
438 end
439 end
440 
441 function [q,l] = Q(mu)
442 % Calculate the observer gain (full dimensional)
443 q = care(A(mu)', C(mu)', Bw(mu)*Bw(mu)', Dw(mu)*Dw(mu)');
444 l = q*C(mu)'*inv(Dw(mu)*Dw(mu)');
445 end
446 function [qhat,lhat,QN,res] = Qhat(mu,N,v)
447 % Calculate the reduced solution to a control ARE (reduced)
448 if exist('v','var')
449  V = v;
450 else
451  V = VQ(:, 1:N);
452 end
453 
454 QN = care(V'*A(mu)'*V, V'*C(mu)', V'*Bw(mu)*Bw(mu)'*V, Dw(mu)*Dw(mu)');
455 qhat = V*QN*V';
456 lhat = qhat*C(mu)'*inv(Dw(mu)*Dw(mu)');
457 if nargout == 4
458  nqhat = norm(qhat,'fro');
459  normA = norm(A(mu),'fro');
460  d = norm(Bw(mu)*Bw(mu)','fro');
461  res = norm( A(mu)*qhat + qhat*A(mu)' - qhat*C(mu)'*inv(Dw(mu)*Dw(mu)')*C(mu)*qhat + Bw(mu)*Bw(mu)', 'fro')/(d + 2*normA*nqhat + nqhat^2*norm(inv(Dw(mu)*Dw(mu)'),'fro')*norm(C(mu)','fro')^2);
462 end
463 end
464 
465 % ================= Implementation of the greedy procedure ================
466 function [V,info] = lrfg( f, fr )
467 % Initialize the basis the most dominant mode of the first vector:
468 [V,~,~] = svds( f(mu_train(1,:)), 1);
469 
470 errs = recalc_errs();
471 tol = 1e-6;
472 info.tol = tol;
473 
474 err_decay = [];
475 while max(errs) > tol
476  [~,i] = max(errs);
477  err_decay = [err_decay max(errs)];
478 
479  display([num2str(i), ' - ', num2str(max(errs))])
480  mu_train(i,:)
481  x = f(mu_train(i,:));
482 
483  % Mimick the low rank factor:
484  % Calculate the low rank factor:
485  [g,d] = eig(x);
486  x = real(g*sqrt(d));
487  x(:, find(diag(sqrt(d)) <= tol)) = [];
488 
489  [u,s,~] = svd(x - V*V'*x, 'econ');
490 
491  s = diag(s);
492  for k = 1:length(s)
493  eE = sum(s(1:k))/sum(s);
494  if eE >= 0.90
495  break
496  end
497  end
498  display([' adding ', num2str(k), ' elements to the basis'])
499  V = orthonormalize_gram_schmidt([V, u(:,1:k)]);
500  errs = recalc_errs();
501  plot(errs);drawnow
502 end
503 
504  function errs = recalc_errs()
505  % Nested function that calculates the "error":
506  errs = zeros(size(mu_train,1),1);
507  parfor i = 1:size(mu_train,1)
508  [x,~,~,res] = fr(mu_train(i,:),[],V);
509  errs(i) = res;
510  end
511  end
512 
513 info.err_decay = err_decay;
514 end
515 end
516 
517 
518 function v = my_2d_vector(n0,f_str)
519 %
520 % Generates a vector v which contains the values of a function f(x,y)
521 % on an equidistant grid in the interior of the unit square. The grid
522 % points are numbered consistently with those used in the function
523 % 'fdm_2d_matrix'.
524 %
525 % This function is just used as an easy way to generate test problems
526 % rather than to solve PDEs.
527 %
528 % Calling sequence:
529 %
530 % v = fdm_2d_vector( n0, f_str)
531 %
532 % Input:
533 %
534 % n0 number of inner grid points in each dimension;
535 % f_str string describing the function f in the space variables 'x'
536 % and 'y', e.g., f_str = 'sin(x+2*y)+3'.
537 %
538 % Output:
539 %
540 % v vector of order n = n0^2 containing the values of f(x,y).
541 %
542 %
543 % LYAPACK 1.0 (Thilo Penzl, May 1999)
544 
545 % Input data not completely checked!
546 
547 na = nargin;
548 
549 if na~=2
550  error('Wrong number of input parameters.');
551 end
552 
553 h = 1.0/(n0+1);
554 
555 n2 = n0*n0;
556 
557 v = zeros(n2,1);
558 
559 i = 0;
560 
561 for iy = 1:n0
562  y = iy*h;
563  for ix = 1:n0
564  x = ix*h;
565  i = i+1;
566  v(i) = eval(f_str);
567  end
568 end
569 end
570 
571 function [A, name] = my_2d_matrix(n0,fx_str,fy_str,g_str,mu_str)
572 %
573 % Generates the stiffness matrix A for the finite difference
574 % discretization (equidistant grid) of the PDE
575 %
576 % mu* laplace(u) - fx du/dx - fy du/dy - g u = r.h.s. on Omega
577 %
578 % u = 0 on dOmega
579 %
580 % Omega = (0,1)x(0,1) (unit square).
581 %
582 % This function is just used as an easy way to generate test problems
583 % rather than to solve PDEs.
584 %
585 % Calling sequence:
586 %
587 % [A, name] = fdm_2d_matrix( n0, fx_str, fy_str, g_str )
588 %
589 % Input:
590 %
591 % n0 number of inner grid points in each dimension;
592 % fx_str string describing the function fx in the space variables
593 % 'x' and 'y', e.g., fx_str = 'sin(x+2*y)+3';
594 % fy_str string describing the function fy in the space variables
595 % 'x' and 'y';
596 % g_str string describing the function g in the space variables
597 % 'x' and 'y'.
598 %
599 % Output:
600 %
601 % A n-x-n sparse stiffness matrix, where n = n0^2;
602 % name string describing the problem.
603 %
604 %
605 % LYAPACK 1.0 (Thilo Penzl, May 1999)
606 
607 % Input data not completely checked!
608 
609 na = nargin;
610 
611 if na~=5
612  mu_str = '1';
613 end
614 
615 name = ['FDM-2D: fx=',fx_str,'; fy=',fy_str,'; g=',g_str];
616 
617 n2 = n0*n0;
618 
619 h = 1.0/(n0+1);
620 
621 h2 = h*h;
622 
623 t1 = 4.0/h2;
624 t2 = -1.0/h2;
625 t3 = 1.0/(2.0*h);
626 
627 len = 5*n2-4*n0;
628 I = zeros(len,1);
629 J = zeros(len,1);
630 S = zeros(len,1);
631 ptr = 0; % Pointer
632 i = 0; % Row Number
633 
634 for iy = 1:n0
635  y = iy*h;
636  for ix = 1:n0
637  x = ix*h;
638 
639  i = i+1;
640  fxv = eval(fx_str);
641  fyv = eval(fy_str);
642  gv = eval(g_str);
643  mu = eval(mu_str);
644 
645  if iy>1
646  ptr = ptr+1; % A(i,i-n)
647  I(ptr) = i;
648  J(ptr) = i-n0;
649  S(ptr) = mu*t2-fyv*t3;
650  end
651 
652  if ix>1
653  ptr = ptr+1; % A(i,i-1)
654  I(ptr) = i;
655  J(ptr) = i-1;
656  S(ptr) = mu*t2-fxv*t3;
657  end
658 
659  ptr = ptr+1; % A(i,i)
660  I(ptr) = i;
661  J(ptr) = i;
662  S(ptr) = mu*t1+gv;
663 
664  if ix<n0
665  ptr = ptr+1; % A(i,i+1)
666  I(ptr) = i;
667  J(ptr) = i+1;
668  S(ptr) = mu*t2+fxv*t3;
669  end
670 
671  if iy<n0
672  ptr = ptr+1; % A(i,i+n0)
673  I(ptr) = i;
674  J(ptr) = i+n0;
675  S(ptr) = mu*t2+fyv*t3;
676  end
677 
678  end
679 end
680 
681 A = -sparse(I,J,S,n2,n2);
682 end
function varargout = filecache_function(funcptr, varargin)
function used for file-caching other function calls.
function cpde2016(input)
CPDE2016 This function reproduces all the results from the CPDE paper.
Definition: cpde2016.m:17
Implementation of the parametric algebraic Riccati equation.
DEFAULTMODEL.
Definition: DefaultModel.m:18