Contents
function LQR_DAE2(problem, level, re, istest)
% Computes a stabilizing feedback by implicitly solving the generalized % Riccati equation for the equivalent projected system on the hidden manifold. % % Inputs: % problem either 'Stokes' or 'NSE' to choose the Stokes demo or the % linearized Navier-Stokes-Equation. % (optional, defaults to 'Stokes') % % level discretization level 1 through 5 % (optional, only used in 'NSE' case, default: 1) % % re Reynolds number 300, 400, or 500 % (optional, only used in 'NSE' case, default: 500) % % istest flag to determine whether this demo runs as a CI test or % interactive demo % (optional, defaults to 0, i.e. interactive demo) % % Note that the 'NSE' option requires additional data available in a % separate 270MB archive and at least the 5th discretization level needs a % considerable amount of main memory installed in your machine. % % 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) %
Set operations
oper = operatormanager('dae_2');
Problem data
if nargin<1, problem='stokes'; end if nargin<2, level=1; end if nargin<3, re=500; end if nargin<4, istest=0; end switch lower(problem) case 'stokes' nin = 5; nout = 5; nx = 10; ny = 10; [eqn.E_,eqn.A_,eqn.B,eqn.C]=stokes_ind2(nin,nout,nx,ny); eqn.haveE=1; st=trace(eqn.E_); % Stokes is FDM discretized, so so this is % the dimension of the velocity space eqn.st=st; eqn.B=eqn.B(1:st,:); eqn.C=eqn.C(:,1:st); case 'nse' [eqn, K0, ~] = mess_get_NSE(re, level); opts.nm.K0 = K0; opts.radi.K0 = K0; otherwise error('input ''problem'' must be either ''NSE'' or ''Stokes'''); end
degrees of freedom: ------------------------------ total : 280 velocity : 180 pressure : 100 n_finite : 81 ------------------------------ Generating FVM matrices... -> Laplacians... -> gradient and divergence operator... -> B1, C1 and v1 velocity nodes... -> B2, C2 and v2 velocity nodes... Setting up system matrices ...
First we run the Newton-ADI Method
opts.norm = 2; % ADI tolerances and maximum iteration number opts.adi.maxiter = 300; opts.adi.res_tol = 1e-12; opts.adi.rel_diff_tol = 1e-16; opts.adi.info = 1; opts.adi.LDL_T=0; eqn.type='T';
n=size(eqn.A_, 1); opts.shifts.num_desired=5;%*nout; opts.shifts.num_Ritz=50; opts.shifts.num_hRitz=25; opts.shifts.method = 'projection'; opts.shifts.b0=ones(n,1);
Newton tolerances and maximum iteration number
opts.nm.maxiter = 20; opts.nm.res_tol = 1e-10; opts.nm.rel_diff_tol = 1e-16; opts.nm.info = 1; opts.nm.projection.freq=0; opts.nm.projection.ortho=1; %opts.nm.projection.meth='care_nwt_fac'; opts.nm.res=struct('maxiter',10,'tol',1e-6,'orth',0); opts.nm.linesearch = 1; opts.nm.inexact = 'superlinear'; opts.nm.tau = 0.1; opts.nm.accumulateRes = 1;
use low-rank Newton-Kleinman-ADI
t_mess_lrnm = tic; outnm = mess_lrnm(eqn, opts, oper); t_elapsed1 =toc(t_mess_lrnm); fprintf(1,'mess_lrnm took %6.2f seconds \n' , t_elapsed1); if not(istest) figure(1); disp(outnm.res); semilogy(outnm.res,'LineWidth',3); title('0 = C^T C + A^T X E + E^T X A - E^T X BB^T X E'); xlabel('number of newton iterations'); ylabel('normalized residual norm'); pause(1); end disp('size outnm.Z:'); disp(size(outnm.Z));
ADI step: 1 normalized residual: 4.749604e-01 relative change in Z: 1.000000e+00 normalized outer residual: 4.748502e-01 ADI step: 2 normalized residual: 2.313480e-01 relative change in Z: 6.115557e-01 normalized outer residual: 2.309720e-01 ADI step: 3 normalized residual: 1.160895e-01 relative change in Z: 4.086903e-01 normalized outer residual: 1.154641e-01 ADI step: 4 normalized residual: 3.274932e-02 relative change in Z: 3.384842e-01 normalized outer residual: 3.186747e-02 ADI step: 5 normalized residual: 4.408328e-03 relative change in Z: 2.031583e-01 normalized outer residual: 3.477590e-03 outer tolerance: 1.188528e-02 NM step: 1 normalized residual: 3.477590e-03 relative change in K: 1.000000e+00 number of ADI steps: 5 ADI step: 1 normalized residual: 4.751369e-01 relative change in Z: 1.000000e+00 normalized outer residual: 4.746418e-01 ADI step: 2 normalized residual: 2.310230e-01 relative change in Z: 6.109172e-01 normalized outer residual: 2.308408e-01 ADI step: 3 normalized residual: 1.156490e-01 relative change in Z: 4.073025e-01 normalized outer residual: 1.155924e-01 ADI step: 4 normalized residual: 3.274385e-02 relative change in Z: 3.354002e-01 normalized outer residual: 3.273882e-02 ADI step: 5 normalized residual: 4.365684e-03 relative change in Z: 2.003457e-01 normalized outer residual: 4.365677e-03 ADI step: 6 normalized residual: 1.459723e-03 relative change in Z: 6.921308e-02 normalized outer residual: 1.459682e-03 ADI step: 7 normalized residual: 9.263155e-05 relative change in Z: 5.492726e-02 normalized outer residual: 9.253389e-05 outer tolerance: 3.863989e-04 NM step: 2 normalized residual: 9.253389e-05 relative change in K: 1.322216e-02 number of ADI steps: 7 ADI step: 1 normalized residual: 5.908896e-01 relative change in Z: 1.000000e+00 normalized outer residual: 5.902469e-01 ADI step: 2 normalized residual: 2.801793e-01 relative change in Z: 6.987959e-01 normalized outer residual: 2.799246e-01 ADI step: 3 normalized residual: 7.784030e-02 relative change in Z: 5.237500e-01 normalized outer residual: 7.780849e-02 ADI step: 4 normalized residual: 6.236314e-04 relative change in Z: 3.257720e-01 normalized outer residual: 6.234901e-04 ADI step: 5 normalized residual: 5.863780e-05 relative change in Z: 3.252883e-02 normalized outer residual: 5.863776e-05 ADI step: 6 normalized residual: 2.832583e-05 relative change in Z: 5.135531e-03 normalized outer residual: 2.832582e-05 ADI step: 7 normalized residual: 2.435343e-05 relative change in Z: 4.159017e-03 normalized outer residual: 2.435343e-05 ADI step: 8 normalized residual: 2.563349e-07 relative change in Z: 1.651384e-03 normalized outer residual: 2.563341e-07 outer tolerance: 3.304782e-06 NM step: 3 normalized residual: 2.563341e-07 relative change in K: 1.145651e-04 number of ADI steps: 8 ADI step: 1 normalized residual: 5.952306e-01 relative change in Z: 1.000000e+00 normalized outer residual: 5.945829e-01 ADI step: 2 normalized residual: 3.403191e-01 relative change in Z: 6.609782e-01 normalized outer residual: 3.399864e-01 ADI step: 3 normalized residual: 1.078645e-01 relative change in Z: 5.650374e-01 normalized outer residual: 1.078086e-01 ADI step: 4 normalized residual: 2.082073e-02 relative change in Z: 3.428347e-01 normalized outer residual: 2.081777e-02 ADI step: 5 normalized residual: 1.759708e-04 relative change in Z: 1.767522e-01 normalized outer residual: 1.759642e-04 ADI step: 6 normalized residual: 8.406516e-05 relative change in Z: 1.521109e-02 normalized outer residual: 8.406342e-05 ADI step: 7 normalized residual: 9.211191e-06 relative change in Z: 1.267083e-02 normalized outer residual: 9.211191e-06 ADI step: 8 normalized residual: 3.101855e-06 relative change in Z: 1.302219e-03 normalized outer residual: 3.101855e-06 ADI step: 9 normalized residual: 2.053260e-08 relative change in Z: 1.875062e-03 normalized outer residual: 2.053260e-08 ADI step: 10 normalized residual: 6.170189e-09 relative change in Z: 8.933796e-05 normalized outer residual: 6.170189e-09 ADI step: 11 normalized residual: 1.864742e-09 relative change in Z: 6.988243e-05 normalized outer residual: 1.864742e-09 ADI step: 12 normalized residual: 1.354435e-09 relative change in Z: 2.042627e-05 normalized outer residual: 1.354435e-09 ADI step: 13 normalized residual: 9.506928e-11 relative change in Z: 4.068726e-05 normalized outer residual: 9.506928e-11 outer tolerance: 3.943602e-09 NM step: 4 normalized residual: 9.506928e-11 relative change in K: 1.797500e-07 number of ADI steps: 13 mess_lrnm took 0.50 seconds 0.0035 0.0001 0.0000 0.0000 size outnm.Z: 180 54

Lets try the RADI method and compare
opts.norm = 2; % RADI-MESS settings opts.shifts.history = opts.shifts.num_desired*size(eqn.C,1); opts.shifts.num_desired = opts.shifts.num_desired; % choose either of the three shift methods, here opts.shifts.method = 'gen-ham-opti'; % opts.shifts.method = 'heur'; % opts.shifts.method = 'projection'; opts.shifts.naive_update_mode = false; % .. Suggest false (smart update is faster; convergence is the same). opts.shifts.info = 0; opts.radi.compute_sol_fac = 1; % Turned on for numerical stability reasons. opts.radi.get_ZZt = 0; opts.radi.maxiter = opts.adi.maxiter; opts.radi.res_tol = opts.nm.res_tol; opts.radi.rel_diff_tol = 0; opts.radi.info = 1; t_mess_lrradi = tic; outradi = mess_lrradi(eqn, opts, oper); t_elapsed2 = toc(t_mess_lrradi); fprintf(1,'mess_lrradi took %6.2f seconds \n' ,t_elapsed2); if not(istest) figure(); semilogy(outradi.res,'LineWidth',3); title('0 = C^TC + A^T X E + E^T X A - E^T X BB^T X E'); xlabel('number of iterations'); ylabel('normalized residual norm'); end
RADI step: 1 normalized residual: 7.035353e-02 RADI step: 2 normalized residual: 1.197302e-02 RADI step: 3 normalized residual: 3.148166e-03 RADI step: 4 normalized residual: 2.370135e-03 RADI step: 5 normalized residual: 6.534372e-05 RADI step: 6 normalized residual: 5.589811e-05 RADI step: 7 normalized residual: 2.820919e-06 RADI step: 8 normalized residual: 2.015254e-07 RADI step: 9 normalized residual: 1.213884e-07 RADI step: 10 normalized residual: 7.433566e-09 RADI step: 11 normalized residual: 4.192365e-09 RADI step: 12 normalized residual: 7.152029e-10 RADI step: 13 normalized residual: 3.266705e-11 mess_lrradi took 0.10 seconds

compare
if istest if min(outnm.res)>=opts.nm.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end if min(outradi.res)>=opts.radi.res_tol, error('MESS:TEST:accuracy','unexpectedly inaccurate result'); end else figure(); ls_nm=[outnm.adi.niter]; ls_radi=1:outradi.niter; semilogy(cumsum(ls_nm),outnm.res,'k--',ls_radi,outradi.res,'b-','LineWidth',3); title('0 = C^T C + A^T X E + E^T X A - E^T X BB^T X E'); xlabel('number of solves with A+p*M'); ylabel('normalized residual norm'); legend('LR-NM','RADI'); end
