Contents
function Lyapunov_rail_LDL_ADI(k,shifts,implicit,istest)
% Computes the solution of the generalized Lyapunov equation via the % low-rank ADI iteration in both ZZ^T [1,2] and LDL^T[3] formulation for the % selective cooling of Steel profiles application described in [4,5,6]. % % Usage: Lyapunov_rail_LDL_ADI(k,shifts,istest) % % Inputs: % % k refinement level of the model to use % (1-4, i.e. 1357-79841Dofs) % (optinal, defaults to 1) % % shifts ADI shift selection strategy. Possible values: % 'heur' Penzl's heuristic shifts % 'Wachspress' Wachspress shifts, optimally solving the dense % shift selection problem. % 'projection' projection shifts [7] % (optional, defaults to 'heur') % % implicit projection shifts with implicit computation of projected A % matrix % % istest flag to determine whether this demo runs as a CI test or % interactive demo % (optional, defaults to 0, i.e. interactive demo) % % References: % % [1] P. Benner, J.-R. Li, T. Penzl, Numerical solution of large-scale % Lyapunov equations, Riccati equations, and linear-quadratic optimal % control problems, Numer. Lin. Alg. Appl. 15 (9) (2008) 755–777. % https://doi.org/10.1002/nla.622. % % [2] P. Kürschner, Efficient low-rank solution of large-scale matrix % equations, Dissertation, Otto-von-Guericke-Universität, Magdeburg, % Germany, shaker Verlag, ISBN 978-3-8440-4385-3 (Apr. 2016). % URL http://hdl.handle.net/11858/00-001M-0000-0029-CE18-2 % % [3] N. Lang, H. Mena, J. Saak, On the benefits of the LDLT factorization % for largescale differential matrix equation solvers, Linear Algebra % Appl. 480 (2015) 44–71. % https://doi.org/10.1016/j.laa.2015.04.006. % % [4] J. Saak, Effiziente numerische Lösung eines % Optimalsteuerungsproblems für die Abkühlung von Stahlprofilen, % Diplomarbeit, Fachbereich 3/Mathematik und Informatik, Universität % Bremen, D-28334 Bremen (Sep. 2003). % % [5] P. Benner, J. Saak, A semi-discretized heat transfer model for % optimal cooling of steel profiles, in: P. Benner, V. Mehrmann, D. % Sorensen (Eds.), Dimension Reduction of Large-Scale Systems, Vol. 45 % of Lect. Notes Comput. Sci. Eng., Springer-Verlag, Berlin/Heidelberg, % Germany, 2005, pp. 353–356. https://doi.org/10.1007/3-540-27909-1_19. % % [6] J. Saak, Efficient numerical solution of large scale algebraic matrix % equations in PDE control and model order reduction, Dissertation, % Technische Universität Chemnitz, Chemnitz, Germany (Jul. 2009). % URL http://nbn-resolving.de/urn:nbn:de:bsz:ch1-200901642 % % [7] P. Benner, P. Kürschner, J. Saak, Self-generating and efficient shift % parameters in ADI methods for large Lyapunov and Sylvester equations, % Electron. Trans. Numer. Anal. 43 (2014) 142–162. % URL http://etna.mcs.kent.edu/volumes/2011-2020/vol43/abstract.php?vol=43&pages=142-162 % % 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); if nargin<1, k=1; end if nargin<2, shifts='wachspress'; end if nargin<3, implicit=0; end if nargin<4, istest=0; end
set operation
oper = operatormanager('default'); % Problem data eqn = getrail(k);
% ADI tolerances and maximum iteration number opts.adi.maxiter = 100; opts.adi.res_tol = 1e-12; opts.adi.rel_diff_tol = 1e-16; opts.adi.info = 1; eqn.type = 'T';
%Heuristic shift parameters via basic Arnoldi n=oper.size(eqn, opts); opts.shifts.num_Ritz=50; opts.shifts.num_hRitz=25; opts.shifts.num_desired = 6; opts.shifts.b0=ones(n,1); switch lower(shifts) case 'heur' opts.shifts.method = 'heur'; case 'wachspress' opts.shifts.method = 'wachspress'; opts.shifts.wachspress = 'T'; case 'projection' opts.shifts.method = 'projection'; opts.shifts.implicitVtAV = implicit; end opts.shifts.p=mess_para(eqn,opts,oper); opts.norm = 'fro';
fprintf('########################\n'); fprintf('# ADI with ZZ^T: \n'); fprintf('########################\n'); tic; out = mess_lradi(eqn, opts, oper); toc; if istest if min(out.res)>=opts.adi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(1); semilogy(out.res); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('size out.Z:'); disp(size(out.Z));
######################## # ADI with ZZ^T: ######################## ADI step: 1 normalized residual: 9.295015e-01 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 8.668178e-01 relative change in Z: 6.761705e-01 ADI step: 3 normalized residual: 8.051179e-01 relative change in Z: 5.486908e-01 ADI step: 4 normalized residual: 7.409268e-01 relative change in Z: 4.816916e-01 ADI step: 5 normalized residual: 6.727096e-01 relative change in Z: 4.393321e-01 ADI step: 6 normalized residual: 6.002287e-01 relative change in Z: 4.078640e-01 ADI step: 7 normalized residual: 5.241929e-01 relative change in Z: 3.813269e-01 ADI step: 8 normalized residual: 4.461882e-01 relative change in Z: 3.570955e-01 ADI step: 9 normalized residual: 3.688459e-01 relative change in Z: 3.338816e-01 ADI step: 10 normalized residual: 2.958357e-01 relative change in Z: 3.108154e-01 ADI step: 11 normalized residual: 2.311456e-01 relative change in Z: 2.874086e-01 ADI step: 12 normalized residual: 1.776061e-01 relative change in Z: 2.639929e-01 ADI step: 13 normalized residual: 1.355442e-01 relative change in Z: 2.419540e-01 ADI step: 14 normalized residual: 1.028509e-01 relative change in Z: 2.230617e-01 ADI step: 15 normalized residual: 7.660017e-02 relative change in Z: 2.079796e-01 ADI step: 16 normalized residual: 5.480322e-02 relative change in Z: 1.954637e-01 ADI step: 17 normalized residual: 3.697245e-02 relative change in Z: 1.833006e-01 ADI step: 18 normalized residual: 2.336999e-02 relative change in Z: 1.699510e-01 ADI step: 19 normalized residual: 1.399895e-02 relative change in Z: 1.554995e-01 ADI step: 20 normalized residual: 8.249446e-03 relative change in Z: 1.414048e-01 ADI step: 21 normalized residual: 5.025404e-03 relative change in Z: 1.287738e-01 ADI step: 22 normalized residual: 3.181876e-03 relative change in Z: 1.165575e-01 ADI step: 23 normalized residual: 2.001619e-03 relative change in Z: 1.024183e-01 ADI step: 24 normalized residual: 1.205763e-03 relative change in Z: 8.534499e-02 ADI step: 25 normalized residual: 7.172900e-04 relative change in Z: 6.705462e-02 ADI step: 26 normalized residual: 4.694384e-04 relative change in Z: 5.153971e-02 ADI step: 27 normalized residual: 3.527571e-04 relative change in Z: 4.274872e-02 ADI step: 28 normalized residual: 2.854815e-04 relative change in Z: 4.066369e-02 ADI step: 29 normalized residual: 2.370529e-04 relative change in Z: 4.160086e-02 ADI step: 30 normalized residual: 1.962232e-04 relative change in Z: 4.324807e-02 ADI step: 31 normalized residual: 1.568807e-04 relative change in Z: 4.518725e-02 ADI step: 32 normalized residual: 1.170897e-04 relative change in Z: 4.679393e-02 ADI step: 33 normalized residual: 7.879545e-05 relative change in Z: 4.699500e-02 ADI step: 34 normalized residual: 4.589439e-05 relative change in Z: 4.478056e-02 ADI step: 35 normalized residual: 2.187953e-05 relative change in Z: 3.950392e-02 ADI step: 36 normalized residual: 7.902182e-06 relative change in Z: 3.136163e-02 ADI step: 37 normalized residual: 1.961263e-06 relative change in Z: 2.177405e-02 ADI step: 38 normalized residual: 3.071648e-07 relative change in Z: 1.294172e-02 ADI step: 39 normalized residual: 3.094205e-08 relative change in Z: 6.510069e-03 ADI step: 40 normalized residual: 1.953366e-09 relative change in Z: 2.662204e-03 ADI step: 41 normalized residual: 4.715330e-11 relative change in Z: 7.698935e-04 ADI step: 42 normalized residual: 4.709828e-13 relative change in Z: 1.230127e-04 Elapsed time is 0.076513 seconds. size out.Z: 1357 252

Set LDL fields
opts.LDL_T = 1; eqn.S = diag([4,4,9,9,16,16]); eqn.G = eqn.C' * diag([4,4,9,9,16,16].^(-0.5));
fprintf('########################\n'); fprintf('# ADI with LDL^T: \n'); fprintf('########################\n'); opts.shifts.p=mess_para(eqn,opts,oper); tic; out1 = mess_lradi(eqn, opts, oper); toc; if istest if min(out.res)>=opts.adi.res_tol error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(2); semilogy(out1.res); xlabel('number of iterations'); ylabel('normalized residual norm'); pause(1); end disp('size out1.Z:'); disp(size(out1.Z));
######################## # ADI with LDL^T: ######################## ADI step: 1 normalized residual: 9.295015e-01 relative change in Z: 1.000000e+00 ADI step: 2 normalized residual: 8.668178e-01 relative change in Z: 6.845616e-01 ADI step: 3 normalized residual: 8.051179e-01 relative change in Z: 5.622064e-01 ADI step: 4 normalized residual: 7.409268e-01 relative change in Z: 4.996237e-01 ADI step: 5 normalized residual: 6.727096e-01 relative change in Z: 4.615544e-01 ADI step: 6 normalized residual: 6.002287e-01 relative change in Z: 4.342783e-01 ADI step: 7 normalized residual: 5.241929e-01 relative change in Z: 4.115193e-01 ADI step: 8 normalized residual: 4.461882e-01 relative change in Z: 3.900572e-01 ADI step: 9 normalized residual: 3.688459e-01 relative change in Z: 3.679366e-01 ADI step: 10 normalized residual: 2.958357e-01 relative change in Z: 3.438585e-01 ADI step: 11 normalized residual: 2.311456e-01 relative change in Z: 3.173376e-01 ADI step: 12 normalized residual: 1.776061e-01 relative change in Z: 2.891777e-01 ADI step: 13 normalized residual: 1.355442e-01 relative change in Z: 2.616553e-01 ADI step: 14 normalized residual: 1.028509e-01 relative change in Z: 2.376189e-01 ADI step: 15 normalized residual: 7.660017e-02 relative change in Z: 2.184736e-01 ADI step: 16 normalized residual: 5.480322e-02 relative change in Z: 2.028716e-01 ADI step: 17 normalized residual: 3.697245e-02 relative change in Z: 1.877678e-01 ADI step: 18 normalized residual: 2.336999e-02 relative change in Z: 1.707601e-01 ADI step: 19 normalized residual: 1.399895e-02 relative change in Z: 1.516709e-01 ADI step: 20 normalized residual: 8.249446e-03 relative change in Z: 1.325671e-01 ADI step: 21 normalized residual: 5.025404e-03 relative change in Z: 1.159418e-01 ADI step: 22 normalized residual: 3.181876e-03 relative change in Z: 1.021666e-01 ADI step: 23 normalized residual: 2.001619e-03 relative change in Z: 8.940539e-02 ADI step: 24 normalized residual: 1.205763e-03 relative change in Z: 7.627503e-02 ADI step: 25 normalized residual: 7.172900e-04 relative change in Z: 6.385690e-02 ADI step: 26 normalized residual: 4.694384e-04 relative change in Z: 5.519119e-02 ADI step: 27 normalized residual: 3.527571e-04 relative change in Z: 5.247123e-02 ADI step: 28 normalized residual: 2.854815e-04 relative change in Z: 5.461069e-02 ADI step: 29 normalized residual: 2.370529e-04 relative change in Z: 5.889037e-02 ADI step: 30 normalized residual: 1.962232e-04 relative change in Z: 6.357819e-02 ADI step: 31 normalized residual: 1.568807e-04 relative change in Z: 6.783239e-02 ADI step: 32 normalized residual: 1.170897e-04 relative change in Z: 7.069701e-02 ADI step: 33 normalized residual: 7.879545e-05 relative change in Z: 7.095376e-02 ADI step: 34 normalized residual: 4.589439e-05 relative change in Z: 6.735728e-02 ADI step: 35 normalized residual: 2.187953e-05 relative change in Z: 5.906471e-02 ADI step: 36 normalized residual: 7.902182e-06 relative change in Z: 4.646421e-02 ADI step: 37 normalized residual: 1.961263e-06 relative change in Z: 3.180240e-02 ADI step: 38 normalized residual: 3.071648e-07 relative change in Z: 1.851713e-02 ADI step: 39 normalized residual: 3.094205e-08 relative change in Z: 9.120345e-03 ADI step: 40 normalized residual: 1.953366e-09 relative change in Z: 3.697835e-03 ADI step: 41 normalized residual: 4.715330e-11 relative change in Z: 1.072806e-03 ADI step: 42 normalized residual: 4.709828e-13 relative change in Z: 1.717290e-04 Elapsed time is 0.074660 seconds. size out1.Z: 1357 252

Difference of Lyapunov solutions
if k<3 % This is mainly for consistency checking on our continuous % integration tests. % NEVER FORM SUCH DYADIC PRODUCTS IN PRODUCTION CODE!!! err = norm(out.Z * out.Z' - out1.Z * out1.D * out1.Z') /... norm(out.Z * out.Z'); fprintf(['Relative difference between solution with and without ' ... 'LDL^T: \t %g\n'], err); if err>1e-12 if implicit shifts = [shifts '(implicit)']; end error('MESS:TEST:accuracy',... ['unexpectedly inaccurate result relative difference',... ' %e > 1e-12 in case %s'],err,shifts); end end
Relative difference between solution with and without LDL^T: 4.11452e-16