Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
If a has some unit roots, the function computes only the solution of the stable subsystem.
INPUTS:
a [double] n*n matrix.
b [double] n*n matrix.
third_argument [double] scalar, if method <= 2 :
qz_criterium = third_argument unit root threshold for eigenvalues in a,
elseif method = 3 :
tol =third_argument the convergence criteria for fixed_point algorithm.
lyapunov_complex_threshold [double] scalar, complex block threshold for the upper triangular matrix T.
method [integer] Scalar, if method=0 [default] then U, T, n and k are not persistent.
method=1 then U, T, n and k are declared as persistent
variables and the schur decomposition is triggered.
method=2 then U, T, n and k are declared as persistent
variables and the schur decomposition is not performed.
method=3 fixed point method
OUTPUTS
x: [double] m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
u: [double] Schur vectors associated with unit roots
ALGORITHM
Uses reordered Schur decomposition
SPECIAL REQUIREMENTS
None0001 function [x,u] = lyapunov_symm(a,b,third_argument,lyapunov_complex_threshold,method, R) 0002 % Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices. 0003 % If a has some unit roots, the function computes only the solution of the stable subsystem. 0004 % 0005 % INPUTS: 0006 % a [double] n*n matrix. 0007 % b [double] n*n matrix. 0008 % third_argument [double] scalar, if method <= 2 : 0009 % qz_criterium = third_argument unit root threshold for eigenvalues in a, 0010 % elseif method = 3 : 0011 % tol =third_argument the convergence criteria for fixed_point algorithm. 0012 % lyapunov_complex_threshold [double] scalar, complex block threshold for the upper triangular matrix T. 0013 % method [integer] Scalar, if method=0 [default] then U, T, n and k are not persistent. 0014 % method=1 then U, T, n and k are declared as persistent 0015 % variables and the schur decomposition is triggered. 0016 % method=2 then U, T, n and k are declared as persistent 0017 % variables and the schur decomposition is not performed. 0018 % method=3 fixed point method 0019 % OUTPUTS 0020 % x: [double] m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem. 0021 % u: [double] Schur vectors associated with unit roots 0022 % 0023 % ALGORITHM 0024 % Uses reordered Schur decomposition 0025 % 0026 % SPECIAL REQUIREMENTS 0027 % None 0028 0029 % Copyright (C) 2006-2010 Dynare Team 0030 % 0031 % This file is part of Dynare. 0032 % 0033 % Dynare is free software: you can redistribute it and/or modify 0034 % it under the terms of the GNU General Public License as published by 0035 % the Free Software Foundation, either version 3 of the License, or 0036 % (at your option) any later version. 0037 % 0038 % Dynare is distributed in the hope that it will be useful, 0039 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0040 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0041 % GNU General Public License for more details. 0042 % 0043 % You should have received a copy of the GNU General Public License 0044 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0045 if nargin<5 0046 method = 0; 0047 end 0048 0049 if method == 3 0050 persistent X method1; 0051 if ~isempty(method1) 0052 method = method1; 0053 end; 0054 tol = third_argument; 0055 fprintf(' [methode=%d] ',method); 0056 if method == 3 0057 %tol = 1e-10; 0058 it_fp = 0; 0059 evol = 100; 0060 if isempty(X) 0061 X = b; 0062 max_it_fp = 2000; 0063 else 0064 max_it_fp = 300; 0065 end; 0066 at = a'; 0067 %fixed point iterations 0068 while evol > tol && it_fp < max_it_fp; 0069 X_old = X; 0070 X = a * X * at + b; 0071 evol = max(sum(abs(X - X_old))); %norm_1 0072 %evol = max(sum(abs(X - X_old)')); %norm_inf 0073 it_fp = it_fp + 1; 0074 end; 0075 fprintf('lyapunov it_fp=%d evol=%g\n',it_fp,evol); 0076 if it_fp >= max_it_fp 0077 disp(['convergence not achieved in solution of Lyapunov equation after ' int2str(it_fp) ' iterations, switching method from 3 to 0']); 0078 method1 = 0; 0079 method = 0; 0080 else 0081 method1 = 3; 0082 x = X; 0083 return; 0084 end; 0085 end; 0086 elseif method == 4 0087 % works only with Matlab System Control toolbox or octave the control package, 0088 if exist('OCTAVE_VERSION') 0089 if ~user_has_octave_forge_package('control') 0090 error('lyapunov=square_root_solver not available; you must install the control package from Octave Forge') 0091 end 0092 else 0093 if ~user_has_matlab_license('control_toolbox') 0094 error('lyapunov=square_root_solver not available; you must install the control system toolbox') 0095 end 0096 end 0097 chol_b = R*chol(b,'lower'); 0098 Rx = dlyapchol(a,chol_b); 0099 x = Rx' * Rx; 0100 return; 0101 end; 0102 0103 qz_criterium = third_argument; 0104 if method 0105 persistent U T k n 0106 else 0107 if exist('U','var') 0108 clear('U','T','k','n') 0109 end 0110 end 0111 0112 u = []; 0113 0114 if size(a,1) == 1 0115 x=b/(1-a*a); 0116 return 0117 end 0118 0119 if method<2 0120 [U,T] = schur(a); 0121 e1 = abs(ordeig(T)) > 2-qz_criterium; 0122 k = sum(e1); % Number of unit roots. 0123 n = length(e1)-k; % Number of stationary variables. 0124 if k > 0 0125 % Selects stable roots 0126 [U,T] = ordschur(U,T,e1); 0127 T = T(k+1:end,k+1:end); 0128 end 0129 end 0130 0131 B = U(:,k+1:end)'*b*U(:,k+1:end); 0132 0133 x = zeros(n,n); 0134 i = n; 0135 0136 while i >= 2 0137 if abs(T(i,i-1))<lyapunov_complex_threshold 0138 if i == n 0139 c = zeros(n,1); 0140 else 0141 c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ... 0142 T(i,i)*T(1:i,i+1:end)*x(i+1:end,i); 0143 end 0144 q = eye(i)-T(1:i,1:i)*T(i,i); 0145 x(1:i,i) = q\(B(1:i,i)+c); 0146 x(i,1:i-1) = x(1:i-1,i)'; 0147 i = i - 1; 0148 else 0149 if i == n 0150 c = zeros(n,1); 0151 c1 = zeros(n,1); 0152 else 0153 c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ... 0154 T(i,i)*T(1:i,i+1:end)*x(i+1:end,i) + ... 0155 T(i,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1); 0156 c1 = T(1:i,:)*(x(:,i+1:end)*T(i-1,i+1:end)') + ... 0157 T(i-1,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1) + ... 0158 T(i-1,i)*T(1:i,i+1:end)*x(i+1:end,i); 0159 end 0160 q = [ eye(i)-T(1:i,1:i)*T(i,i) , -T(1:i,1:i)*T(i,i-1) ; ... 0161 -T(1:i,1:i)*T(i-1,i) , eye(i)-T(1:i,1:i)*T(i-1,i-1) ]; 0162 z = q\[ B(1:i,i)+c ; B(1:i,i-1) + c1 ]; 0163 x(1:i,i) = z(1:i); 0164 x(1:i,i-1) = z(i+1:end); 0165 x(i,1:i-1) = x(1:i-1,i)'; 0166 x(i-1,1:i-2) = x(1:i-2,i-1)'; 0167 i = i - 2; 0168 end 0169 end 0170 if i == 1 0171 c = T(1,:)*(x(:,2:end)*T(1,2:end)') + T(1,1)*T(1,2:end)*x(2:end,1); 0172 x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1)); 0173 end 0174 x = U(:,k+1:end)*x*U(:,k+1:end)'; 0175 u = U(:,1:k);