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; defalts 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 orde 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 program is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation; either version 2 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, see <http://www.gnu.org/licenses/>. % % Copyright (C) Jens Saak, Martin Koehler, Peter Benner and others % 2009-2020 %
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 opts.shifts.p=mess_para(eqn,opts,oper); disp('(Initial) ADI shifts'); disp(opts.shifts.p);
(Initial) ADI shifts -2.065258190191611e+04 + 0.000000000000000e+00i -1.937292407469257e+04 + 0.000000000000000e+00i -1.749376770662731e+04 + 0.000000000000000e+00i -1.330511046698838e+04 + 0.000000000000000e+00i -9.263930276803778e+03 + 0.000000000000000e+00i -6.671163618267795e+03 + 0.000000000000000e+00i -1.698501218018387e+03 + 0.000000000000000e+00i -4.649900486215604e+02 + 0.000000000000000e+00i -2.603557887074633e+02 + 0.000000000000000e+00i -1.918393347707617e+02 + 0.000000000000000e+00i -1.112802038431261e+02 + 0.000000000000000e+00i -4.453715272674594e+03 + 9.448801986471070e+02i -4.453715272674594e+03 - 9.448801986471070e+02i -3.511454699403188e+03 + 1.017535965450999e+03i -3.511454699403188e+03 - 1.017535965450999e+03i -1.941507072350244e+03 + 7.426412079916519e+02i -1.941507072350244e+03 - 7.426412079916519e+02i -9.796882754126804e+02 + 9.470372612793864e+02i -9.796882754126804e+02 - 9.470372612793864e+02i -8.124319370443204e+02 + 4.557291003561135e+02i -8.124319370443204e+02 - 4.557291003561135e+02i -3.712015209114301e+02 + 2.393803494728078e+02i -3.712015209114301e+02 - 2.393803494728078e+02i -1.392548614286558e+02 + 2.618919038961343e+01i -1.392548614286558e+02 - 2.618919038961343e+01i
%observability eqn.type='N'; tic; outB = mess_lradi(eqn, opts, oper); toc; if istest if min(outB.res)>=opts.adi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); semilogy(outB.res); title('AX + XA^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 Elapsed time is 0.110372 seconds. size outB.Z: 2500 25

%controllability eqn.type = 'T'; tic; outC = mess_lradi(eqn, opts, oper); toc; if istest if min(outC.res)>=opts.adi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(2); semilogy(outC.res); title('A^TX + XA = -C^TC'); 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 Elapsed time is 0.159134 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); 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 = Columns 1 through 3 -8.829057200166768e+01 -1.019807420789014e+00 -9.445425120873438e+01 1.019809700641548e+00 -6.032814822854384e-02 -2.738130320613079e+02 -9.445425124221941e+01 2.738130476572438e+02 -6.449423333932912e+02 -3.895979695274298e+01 5.269097785916864e+00 -4.720956862483179e+02 2.657259610089089e+00 -3.229743304497948e-01 3.544604944344591e+01 -2.007136244229511e+01 2.505374132701651e+00 -2.609151425920568e+02 -7.941125513916729e+00 9.846878060831060e-01 -1.038824212154909e+02 -6.650291192945058e-01 8.230925023975408e-02 -8.715254979717750e+00 -1.933456615594407e+00 2.394417042909453e-01 -2.532332230991804e+01 8.840248681283010e-02 -1.094672911079475e-02 1.157997721452221e+00 -9.859754967001821e-01 1.220951580802673e-01 -1.291468089680958e+01 5.933504024874833e-02 -7.346286808055215e-03 7.772053996853828e-01 Columns 4 through 6 -3.895979674387746e+01 -2.657256449973046e+00 -2.007136239791246e+01 -5.269087839516895e+00 -3.229719701910339e-01 -2.505368460790542e+00 -4.720956816023990e+02 -3.544601011733786e+01 -2.609151414699180e+02 -1.198154522101892e+03 -2.091589392939177e+02 -1.124891937041254e+03 2.091591199188712e+02 -2.488425796556966e+01 -6.521257009031148e+02 -1.124891937831908e+03 6.521260437433831e+02 -3.364971173014024e+03 -4.790746202810799e+02 1.603658658652016e+02 -2.275787585528893e+03 -4.100898102733039e+01 1.212810319522633e+01 -2.385831669273557e+02 -1.183748028283387e+02 3.640648114391952e+01 -6.412041706538694e+02 5.421138574529287e+00 -1.652442001380978e+00 2.982625078550474e+01 -6.041914951522367e+01 1.849139038149388e+01 -3.300697867036573e+02 3.636575821763832e+00 -1.112034912791116e+00 1.990042621069726e+01 Columns 7 through 9 -7.941124327467590e+00 6.650193716657503e-01 -1.933454879063501e+00 -9.846840714258981e-01 8.230685197982039e-02 -2.394416150829136e-01 -1.038824070795128e+02 8.715128575640250e+00 -2.532329868077609e+01 -4.790745158538591e+02 4.100836444178677e+01 -1.183747120669450e+02 -1.603657970798759e+02 1.212800776525438e+01 -3.640633077964564e+01 -2.275786982086108e+03 2.385800241693510e+02 -6.412036399833823e+02 -3.102754701000881e+03 8.481464517018636e+02 -1.421087701775793e+03 -8.481508882672374e+02 -5.616941834441856e+01 3.902170010744658e+02 -1.421090573132168e+03 -3.902234082395105e+02 -2.913330645316168e+03 7.184392791618765e+01 1.358213962932752e+01 6.834788023398814e+02 -7.647495486709025e+02 -1.698252188863975e+02 -2.658153174937496e+03 4.653399280207658e+01 9.935020544111936e+00 1.896467105169808e+02 Columns 10 through 12 -8.839465760239280e-02 -9.859739496303747e-01 -5.917363441388186e-02 -1.094541017677039e-02 -1.220951463707010e-01 -7.327480813125683e-03 -1.157895568422158e+00 -1.291466004288380e+01 -7.750897632465272e-01 -5.420651529498241e+00 -6.041906542795733e+01 -3.626712140989104e+00 -1.652327004379968e+00 -1.849128882615374e+01 -1.108880341142835e+00 -2.982361743621728e+01 -3.300693060137281e+02 -1.984645671967601e+01 -7.183851286384869e+01 -7.647474632338972e+02 -4.640645751367969e+01 1.358163344690188e+01 1.698209970079201e+02 9.903309425863881e+00 -6.834610280793736e+02 -2.658152870223816e+03 -1.892032426369110e+02 -9.978790743996990e+00 -2.758794823644630e+02 -1.220734908508392e+01 2.759007134037782e+02 -6.430581526719097e+03 -1.505522035506696e+03 -1.226347412651648e+01 1.506870839135303e+03 -4.757197654714582e+01


