function [Ar, Br, Cr] = bt_mor_FDM_tol(tol, n0, shifts, istest)
% bt_mor_FDM_tol computes a reduced order model via the standard Lyapunov % balanced truncation (see e.g. [1]) for a finite difference discretized % convection diffusion model on the unit square described in [2]. % % Usage: % [Ar, Br, Cr] = bt_mor_FDM_tol(tol,max_ord,n0,test) % % Inputs % % tol truncation tolerance for the Hankel singular values % (optional; defaults to 1e-6) % % n0 n0^2 gives the dimension of the original model, i.e. n0 is % the number of degrees of freedom, i.e. grid points, per % spatial direction % (optional; defaults to 50) % % shifts shift selection used in ADI; possible choices: % 'heur' : Penzl heuristic shifts % 'projection' : projection shifts using the last columns % of the solution factor % (optional, defaults to 'heur') % % istest flag to determine whether this demo runs as a CI test or % interactive demo % (optional, defaults to 0, i.e. interactive demo) % % Outputs % % Ar, Br, Cr the reduced order system matrices. % % References % [1] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, Vol. % 6 of Adv. Des. Control, SIAM Publications, Philadelphia, PA, 2005. % https://doi.org/10.1137/1.9780898718713 % % [2] T. Penzl, LyaPack Users Guide, Tech. Rep. SFB393/00-33, % Sonderforschungsbereich 393 Numerische Simulation auf massiv % parallelen Rechnern, TU Chemnitz, 09107 Chemnitz, Germany, % available from http://www.tu-chemnitz.de/sfb393/sfb00pr.html (2000). % % This file is part of the M-M.E.S.S. project % (http://www.mpi-magdeburg.mpg.de/projects/mess). % Copyright © 2009-2022 Jens Saak, Martin Koehler, Peter Benner and others. % All rights reserved. % License: BSD 2-Clause License (see COPYING) %
narginchk(0, 4); % BT tolerance and maximum order for the ROM if nargin < 1, tol = 1e-6; end if nargin < 2, n0 = 50; end if nargin < 3, shifts = 'heur'; end if nargin < 4, istest = 0; end % ADI tolerance and maximum iteration number opts.adi.maxiter = 100; opts.adi.res_tol = 1e-9; opts.adi.rel_diff_tol = 1e-16; opts.adi.info = 1; opts.norm = 'fro';
operations
oper = operatormanager('default'); % Problem data eqn.A_ = fdm_2d_matrix(n0, '10*x','100*y','0'); eqn.B = fdm_2d_vector(n0, '.1<x<=.3'); eqn.C = fdm_2d_vector(n0, '.7<x<=.9'); eqn.C = eqn.C'; eqn.haveE = 0; % problem dimension n=oper.size(eqn, opts);
Heuristic Parameters via basic Arnoldi or Projection shifts
opts.shifts.method = shifts; switch opts.shifts.method case 'heur' opts.shifts.num_desired = 25; opts.shifts.num_Ritz = 50; opts.shifts.num_hRitz = 25; opts.shifts.b0 = ones(n, 1); case 'projection' opts.shifts.method = 'projection'; opts.shifts.num_desired = 10; otherwise warning('MESS:Test:invalid_parameter_fallback',... 'invalid shift selection, falling back to "heur"'); opts.shifts.num_desired = 25; opts.shifts.num_Ritz = 50; opts.shifts.num_hRitz = 25; opts.shifts.b0 = ones(n, 1); opts.shifts.method = 'heur'; end
controllability
eqn.type = 'N'; t_mess_lradi = tic; outB = mess_lradi(eqn, opts, oper); t_elapsed1 = toc(t_mess_lradi); fprintf(1,'mess_lradi took %6.2f seconds \n' ,t_elapsed1); if istest if min(outB.res)>=opts.adi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); semilogy(outB.res,'LineWidth',3); title('A X + X A^T = -BB^T'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('size outB.Z:'); disp(size(outB.Z));
ADI step: 1 normalized residual: 8.428057e-01 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 7.378016e-01 relative change in Z: 6.917507e-01 ADI step: 3 normalized residual: 6.508479e-01 relative change in Z: 5.637997e-01 ADI step: 4 normalized residual: 5.614350e-01 relative change in Z: 5.169589e-01 ADI step: 5 normalized residual: 4.639724e-01 relative change in Z: 4.948088e-01 ADI step: 6 normalized residual: 3.674147e-01 relative change in Z: 4.645217e-01 ADI step: 7 normalized residual: 1.830519e-01 relative change in Z: 5.907596e-01 ADI step: 8 normalized residual: 4.914772e-02 relative change in Z: 5.625659e-01 ADI step: 9 normalized residual: 1.439251e-02 relative change in Z: 3.708279e-01 ADI step: 10 normalized residual: 2.640914e-03 relative change in Z: 2.116408e-01 ADI step: 11 normalized residual: 1.583523e-04 relative change in Z: 7.793154e-02 ADI step: 13 normalized residual: 1.066545e-04 relative change in Z: 8.368038e-03 ADI step: 15 normalized residual: 7.634587e-05 relative change in Z: 7.751726e-03 ADI step: 17 normalized residual: 4.229441e-05 relative change in Z: 8.133399e-03 ADI step: 19 normalized residual: 1.675468e-05 relative change in Z: 6.178842e-03 ADI step: 21 normalized residual: 2.303907e-06 relative change in Z: 4.177396e-03 ADI step: 23 normalized residual: 1.919842e-07 relative change in Z: 1.872833e-03 ADI step: 25 normalized residual: 1.238977e-10 relative change in Z: 6.576621e-04 mess_lradi took 0.32 seconds size outB.Z: 2500 25

observability
eqn.type = 'T'; t_mess_lradi = tic; outC = mess_lradi(eqn, opts, oper); t_elapsed2 = toc(t_mess_lradi); fprintf(1, 'mess_lradi took %6.2f seconds \n', t_elapsed2); if istest if min(outC.res) >= opts.adi.res_tol error('MESS:TEST:accuracy', 'unexpectedly inaccurate result'); end else figure(2); semilogy(outC.res, 'LineWidth', 3); title('A^T X + X A = -C^T C'); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('size outC.Z:'); disp(size(outC.Z));
ADI step: 1 normalized residual: 9.739434e-01 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 9.600200e-01 relative change in Z: 7.149146e-01 ADI step: 3 normalized residual: 9.507972e-01 relative change in Z: 5.989262e-01 ADI step: 4 normalized residual: 9.434197e-01 relative change in Z: 5.644110e-01 ADI step: 5 normalized residual: 9.378728e-01 relative change in Z: 5.588782e-01 ADI step: 6 normalized residual: 9.355761e-01 relative change in Z: 5.491134e-01 ADI step: 7 normalized residual: 9.533528e-01 relative change in Z: 7.354621e-01 ADI step: 8 normalized residual: 1.080148e+00 relative change in Z: 8.167288e-01 ADI step: 9 normalized residual: 1.086516e+00 relative change in Z: 7.398403e-01 ADI step: 10 normalized residual: 4.755716e-01 relative change in Z: 5.745268e-01 ADI step: 11 normalized residual: 5.254634e-02 relative change in Z: 2.652358e-01 ADI step: 13 normalized residual: 3.985457e-02 relative change in Z: 4.242675e-02 ADI step: 15 normalized residual: 2.784426e-02 relative change in Z: 4.005181e-02 ADI step: 17 normalized residual: 1.387804e-02 relative change in Z: 4.044519e-02 ADI step: 19 normalized residual: 4.202789e-03 relative change in Z: 2.932766e-02 ADI step: 21 normalized residual: 2.110828e-04 relative change in Z: 1.612450e-02 ADI step: 23 normalized residual: 7.179830e-07 relative change in Z: 3.008242e-03 ADI step: 25 normalized residual: 2.305240e-07 relative change in Z: 1.640317e-04 ADI step: 26 normalized residual: 1.970927e-07 relative change in Z: 3.081957e-05 ADI step: 27 normalized residual: 1.682114e-07 relative change in Z: 2.941062e-05 ADI step: 28 normalized residual: 1.425346e-07 relative change in Z: 2.854049e-05 ADI step: 29 normalized residual: 1.163449e-07 relative change in Z: 2.983788e-05 ADI step: 30 normalized residual: 8.919103e-08 relative change in Z: 3.178993e-05 ADI step: 31 normalized residual: 6.415080e-08 relative change in Z: 3.227039e-05 ADI step: 32 normalized residual: 2.271836e-08 relative change in Z: 4.511281e-05 ADI step: 33 normalized residual: 2.096728e-09 relative change in Z: 3.126206e-05 ADI step: 34 normalized residual: 8.476203e-10 relative change in Z: 7.747729e-06 mess_lradi took 0.65 seconds size outC.Z: 2500 34

opts.srm.tol = tol; opts.srm.info = 1; [TL, TR, HSV, eqn, opts, ~] = mess_square_root_method(eqn, opts, oper,... outB.Z, outC.Z); Ar = TL' * (eqn.A_ * TR); Br = TL' * eqn.B; Cr = eqn.C * TR; opts.sigma.nsample = 200; % 200 frequency samples opts.sigma.fmin = -3; % min. frequency 1e-3 opts.sigma.fmax = 4; % max. frequency 1e4 if istest opts.sigma.info = 1; % no output else opts.sigma.info = 2; % show messages and plots end ROM.A = Ar; ROM.B = Br; ROM.C = Cr; ROM.E = eye(size(ROM.A,1)); out = mess_sigma_plot(eqn, opts, oper, ROM); err = out.err; if istest if max(err)>tol error('MESS:TEST:accuracy', 'unexpectedly inaccurate result'); end else figure; semilogy(HSV, 'LineWidth', 3); title('Computed Hankel singular values'); xlabel('index'); ylabel('magnitude'); end
reduced system order: 12 (max possible/allowed: 25/2500) Computing TFMs of original and reduced order systems and MOR errors Step 20 / 200 Step 40 / 200 Step 60 / 200 Step 80 / 200 Step 100 / 200 Step 120 / 200 Step 140 / 200 Step 160 / 200 Step 180 / 200 Step 200 / 200 ans = 1.0e+03 * Columns 1 through 7 -0.0883 -0.0010 -0.0945 -0.0390 -0.0027 -0.0201 -0.0079 0.0010 -0.0001 -0.2738 -0.0053 -0.0003 -0.0025 -0.0010 -0.0945 0.2738 -0.6449 -0.4721 -0.0354 -0.2609 -0.1039 -0.0390 0.0053 -0.4721 -1.1982 -0.2092 -1.1249 -0.4791 0.0027 -0.0003 0.0354 0.2092 -0.0249 -0.6521 -0.1604 -0.0201 0.0025 -0.2609 -1.1249 0.6521 -3.3650 -2.2758 -0.0079 0.0010 -0.1039 -0.4791 0.1604 -2.2758 -3.1028 -0.0007 0.0001 -0.0087 -0.0410 0.0121 -0.2386 -0.8482 -0.0019 0.0002 -0.0253 -0.1184 0.0364 -0.6412 -1.4211 0.0001 -0.0000 0.0012 0.0054 -0.0017 0.0298 0.0718 -0.0010 0.0001 -0.0129 -0.0604 0.0185 -0.3301 -0.7647 0.0001 -0.0000 0.0008 0.0036 -0.0011 0.0199 0.0465 Columns 8 through 12 0.0007 -0.0019 -0.0001 -0.0010 -0.0001 0.0001 -0.0002 -0.0000 -0.0001 -0.0000 0.0087 -0.0253 -0.0012 -0.0129 -0.0008 0.0410 -0.1184 -0.0054 -0.0604 -0.0036 0.0121 -0.0364 -0.0017 -0.0185 -0.0011 0.2386 -0.6412 -0.0298 -0.3301 -0.0198 0.8481 -1.4211 -0.0718 -0.7647 -0.0464 -0.0562 0.3902 0.0136 0.1698 0.0099 -0.3902 -2.9133 -0.6835 -2.6582 -0.1892 0.0136 0.6835 -0.0100 -0.2759 -0.0122 -0.1698 -2.6582 0.2759 -6.4306 -1.5055 0.0099 0.1896 -0.0123 1.5069 -0.0476


