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 (c) 2009-2025 Jens Saak, Martin Koehler, Peter Benner and others. % All rights reserved. % License: BSD 2-Clause License (see COPYING) %
Set operations
opts = struct();
[oper, opts] = operatormanager(opts, '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 = false; 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, opts); eqn.haveE = true; n_ode = trace(eqn.E_); % Stokes is FDM discretized, so so this is % the dimension of the velocity space eqn.manifold_dim = n_ode; eqn.B = eqn.B(1:n_ode, :); eqn.C = eqn.C(:, 1:n_ode); case 'nse' [eqn, K0, ~] = mess_get_NSE(re, level); opts.nm.K0 = K0; opts.radi.K0 = K0; otherwise mess_err(opts, 'illegal_input', ... '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 = false; 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 = true; % in case you want to e.g. specify the factored Newton solver for % the projected equations uncomment the following % opts.nm.projection.meth = 'care_nwt_fac'; opts.nm.res = struct('maxiter', 10, ... 'tol', 1e-6, ... 'orth', 0); opts.nm.linesearch = true; opts.nm.inexact = 'superlinear'; opts.nm.tau = 0.1; opts.nm.accumulateRes = true;
use low-rank Newton-Kleinman-ADI
t_mess_lrnm = tic; outnm = mess_lrnm(eqn, opts, oper); t_elapsed1 = toc(t_mess_lrnm); mess_fprintf(opts, 'mess_lrnm took %6.2f seconds \n\n', t_elapsed1); if not(istest) figure(1); 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 [mZ, nZ] = size(outnm.Z); mess_fprintf(opts, 'size outnm.Z: %d x %d\n\n', mZ, nZ);
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.459719e-03 relative change in Z: 6.921313e-02 normalized outer residual: 1.459678e-03 ADI step: 7 normalized residual: 9.263129e-05 relative change in Z: 5.492719e-02 normalized outer residual: 9.253362e-05 outer tolerance: 3.863989e-04 NM step: 2 normalized residual: 9.253362e-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.832581e-05 relative change in Z: 5.135530e-03 normalized outer residual: 2.832580e-05 ADI step: 7 normalized residual: 2.435341e-05 relative change in Z: 4.159018e-03 normalized outer residual: 2.435341e-05 ADI step: 8 normalized residual: 2.563225e-07 relative change in Z: 1.651385e-03 normalized outer residual: 2.563217e-07 outer tolerance: 3.304772e-06 NM step: 3 normalized residual: 2.563217e-07 relative change in K: 1.145647e-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.406515e-05 relative change in Z: 1.521109e-02 normalized outer residual: 8.406341e-05 ADI step: 7 normalized residual: 9.211192e-06 relative change in Z: 1.267083e-02 normalized outer residual: 9.211192e-06 ADI step: 8 normalized residual: 3.101854e-06 relative change in Z: 1.302219e-03 normalized outer residual: 3.101854e-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.170173e-09 relative change in Z: 8.933790e-05 normalized outer residual: 6.170173e-09 ADI step: 11 normalized residual: 1.864738e-09 relative change in Z: 6.988233e-05 normalized outer residual: 1.864738e-09 ADI step: 12 normalized residual: 1.354432e-09 relative change in Z: 2.042626e-05 normalized outer residual: 1.354432e-09 ADI step: 13 normalized residual: 9.506907e-11 relative change in Z: 4.068722e-05 normalized outer residual: 9.506907e-11 outer tolerance: 3.943411e-09 NM step: 4 normalized residual: 9.506907e-11 relative change in K: 1.797487e-07 number of ADI steps: 13 mess_lrnm took 0.37 seconds size outnm.Z: 180 x 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 = true; % Turned on for numerical stability reasons. opts.radi.get_ZZt = false; 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); mess_fprintf(opts, '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: 1.190239e-01 RADI step: 2 normalized residual: 2.025592e-02 RADI step: 3 normalized residual: 5.326059e-03 RADI step: 4 normalized residual: 4.009789e-03 RADI step: 5 normalized residual: 1.105483e-04 RADI step: 6 normalized residual: 9.456828e-05 RADI step: 7 normalized residual: 4.772424e-06 RADI step: 8 normalized residual: 3.409402e-07 RADI step: 9 normalized residual: 2.053645e-07 RADI step: 10 normalized residual: 1.257609e-08 RADI step: 11 normalized residual: 7.092632e-09 RADI step: 12 normalized residual: 1.209978e-09 RADI step: 13 normalized residual: 5.526603e-11 mess_lrradi took 0.14 seconds

compare
if istest if min(outnm.res) >= opts.nm.res_tol mess_err(opts, 'TEST:accuracy', ... 'unexpectedly inaccurate result'); end if min(outradi.res) >= opts.radi.res_tol mess_err(opts, '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
