1function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,method,debug) % --*-- Unitary tests --*-- 2% Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices. 3% If a has some unit roots, the function computes only the solution of the stable subsystem. 4% 5% INPUTS: 6% a [double] n*n matrix. 7% b [double] n*n matrix. 8% qz_criterium [double] unit root threshold for eigenvalues 9% lyapunov_fixed_point_tol [double] convergence criteria for fixed_point algorithm. 10% lyapunov_complex_threshold [double] scalar, complex block threshold for the upper triangular matrix T. 11% method [integer] Scalar, if method=0 [default] then U, T, n and k are not persistent. 12% method=1 then U, T, n and k are declared as persistent 13% variables and the Schur decomposition is triggered. 14% method=2 then U, T, n and k are declared as persistent 15% variables and the Schur decomposition is not performed. 16% method=3 fixed point method 17% OUTPUTS 18% x: [double] m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem. 19% u: [double] Schur vectors associated with unit roots 20% 21% ALGORITHM 22% Uses reordered Schur decomposition (Bartels-Stewart algorithm) 23% [method<3] or a fixed point algorithm (method==3) 24% 25% SPECIAL REQUIREMENTS 26% None 27 28% Copyright (C) 2006-2017 Dynare Team 29% 30% This file is part of Dynare. 31% 32% Dynare is free software: you can redistribute it and/or modify 33% it under the terms of the GNU General Public License as published by 34% the Free Software Foundation, either version 3 of the License, or 35% (at your option) any later version. 36% 37% Dynare is distributed in the hope that it will be useful, 38% but WITHOUT ANY WARRANTY; without even the implied warranty of 39% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 40% GNU General Public License for more details. 41% 42% You should have received a copy of the GNU General Public License 43% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 44 45if nargin<6 || isempty(method) 46 method = 0; 47end 48 49if nargin<7 50 debug = 0; 51end 52 53if method == 3 54 persistent X method1; 55 if ~isempty(method1) 56 method = method1; 57 end 58 if debug 59 fprintf('lyapunov_symm:: [method=%d] \n',method); 60 end 61 if method == 3 62 %tol = 1e-10; 63 it_fp = 0; 64 evol = 100; 65 if isempty(X) || length(X)~=length(b) 66 X = b; 67 max_it_fp = 2000; 68 else 69 max_it_fp = 300; 70 end 71 at = a'; 72 %fixed point iterations 73 while evol > lyapunov_fixed_point_tol && it_fp < max_it_fp 74 X_old = X; 75 X = a * X * at + b; 76 evol = max(sum(abs(X - X_old))); %norm_1 77 %evol = max(sum(abs(X - X_old)')); %norm_inf 78 it_fp = it_fp + 1; 79 end 80 if debug 81 fprintf('lyapunov_symm:: lyapunov fixed_point iterations=%d norm=%g\n',it_fp,evol); 82 end 83 if it_fp >= max_it_fp 84 disp(['lyapunov_symm:: convergence not achieved in solution of Lyapunov equation after ' int2str(it_fp) ' iterations, switching method from 3 to 0']); 85 method1 = 0; 86 method = 0; 87 else 88 method1 = 3; 89 x = X; 90 return 91 end 92 end 93end 94 95if method 96 persistent U T k n 97else 98 % if exist('U','var') 99 % clear('U','T','k','n') 100 % end 101end 102 103u = []; 104 105if size(a,1) == 1 106 x=b/(1-a*a); 107 return 108end 109 110if method<2 111 [U,T] = schur(a); 112 e1 = abs(ordeig(T)) > 2-qz_criterium; 113 k = sum(e1); % Number of unit roots. 114 n = length(e1)-k; % Number of stationary variables. 115 if k > 0 116 % Selects stable roots 117 [U,T] = ordschur(U,T,e1); 118 T = T(k+1:end,k+1:end); 119 end 120end 121 122B = U(:,k+1:end)'*b*U(:,k+1:end); 123 124x = zeros(n,n); 125i = n; 126 127while i >= 2 128 if abs(T(i,i-1))<lyapunov_complex_threshold 129 if i == n 130 c = zeros(n,1); 131 else 132 c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ... 133 T(i,i)*T(1:i,i+1:end)*x(i+1:end,i); 134 end 135 q = eye(i)-T(1:i,1:i)*T(i,i); 136 x(1:i,i) = q\(B(1:i,i)+c); 137 x(i,1:i-1) = x(1:i-1,i)'; 138 i = i - 1; 139 else 140 if i == n 141 c = zeros(n,1); 142 c1 = zeros(n,1); 143 else 144 c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ... 145 T(i,i)*T(1:i,i+1:end)*x(i+1:end,i) + ... 146 T(i,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1); 147 c1 = T(1:i,:)*(x(:,i+1:end)*T(i-1,i+1:end)') + ... 148 T(i-1,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1) + ... 149 T(i-1,i)*T(1:i,i+1:end)*x(i+1:end,i); 150 end 151 q = [ eye(i)-T(1:i,1:i)*T(i,i) , -T(1:i,1:i)*T(i,i-1) ; ... 152 -T(1:i,1:i)*T(i-1,i) , eye(i)-T(1:i,1:i)*T(i-1,i-1) ]; 153 z = q\[ B(1:i,i)+c ; B(1:i,i-1) + c1 ]; 154 x(1:i,i) = z(1:i); 155 x(1:i,i-1) = z(i+1:end); 156 x(i,1:i-1) = x(1:i-1,i)'; 157 x(i-1,1:i-2) = x(1:i-2,i-1)'; 158 i = i - 2; 159 end 160end 161if i == 1 162 c = T(1,:)*(x(:,2:end)*T(1,2:end)') + T(1,1)*T(1,2:end)*x(2:end,1); 163 x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1)); 164end 165x = U(:,k+1:end)*x*U(:,k+1:end)'; 166u = U(:,1:k);