/usr/share/dynare/matlab/lyapunov_symm.m is in dynare-common 4.4.1-1build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 | function [x,u] = lyapunov_symm(a,b,third_argument,lyapunov_complex_threshold,method, R)
% 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
% None
% Copyright (C) 2006-2012 Dynare Team
%
% This file is part of Dynare.
%
% Dynare 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 3 of the License, or
% (at your option) any later version.
%
% Dynare 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 Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<5
method = 0;
end
if method == 3
persistent X method1;
if ~isempty(method1)
method = method1;
end;
tol = third_argument;
fprintf(' [methode=%d] ',method);
if method == 3
%tol = 1e-10;
it_fp = 0;
evol = 100;
if isempty(X)
X = b;
max_it_fp = 2000;
else
max_it_fp = 300;
end;
at = a';
%fixed point iterations
while evol > tol && it_fp < max_it_fp;
X_old = X;
X = a * X * at + b;
evol = max(sum(abs(X - X_old))); %norm_1
%evol = max(sum(abs(X - X_old)')); %norm_inf
it_fp = it_fp + 1;
end;
fprintf('lyapunov it_fp=%d evol=%g\n',it_fp,evol);
if it_fp >= max_it_fp
disp(['convergence not achieved in solution of Lyapunov equation after ' int2str(it_fp) ' iterations, switching method from 3 to 0']);
method1 = 0;
method = 0;
else
method1 = 3;
x = X;
return;
end;
end;
elseif method == 4
% works only with Matlab System Control toolbox or octave the control package,
if isoctave
if ~user_has_octave_forge_package('control')
error('lyapunov=square_root_solver not available; you must install the control package from Octave Forge')
end
else
if ~user_has_matlab_license('control_toolbox')
error('lyapunov=square_root_solver not available; you must install the control system toolbox')
end
end
chol_b = R*chol(b,'lower');
Rx = dlyapchol(a,chol_b);
x = Rx' * Rx;
return;
end;
qz_criterium = third_argument;
if method
persistent U T k n
else
if exist('U','var')
clear('U','T','k','n')
end
end
u = [];
if size(a,1) == 1
x=b/(1-a*a);
return
end
if method<2
[U,T] = schur(a);
e1 = abs(ordeig(T)) > 2-qz_criterium;
k = sum(e1); % Number of unit roots.
n = length(e1)-k; % Number of stationary variables.
if k > 0
% Selects stable roots
[U,T] = ordschur(U,T,e1);
T = T(k+1:end,k+1:end);
end
end
B = U(:,k+1:end)'*b*U(:,k+1:end);
x = zeros(n,n);
i = n;
while i >= 2
if abs(T(i,i-1))<lyapunov_complex_threshold
if i == n
c = zeros(n,1);
else
c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ...
T(i,i)*T(1:i,i+1:end)*x(i+1:end,i);
end
q = eye(i)-T(1:i,1:i)*T(i,i);
x(1:i,i) = q\(B(1:i,i)+c);
x(i,1:i-1) = x(1:i-1,i)';
i = i - 1;
else
if i == n
c = zeros(n,1);
c1 = zeros(n,1);
else
c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ...
T(i,i)*T(1:i,i+1:end)*x(i+1:end,i) + ...
T(i,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1);
c1 = T(1:i,:)*(x(:,i+1:end)*T(i-1,i+1:end)') + ...
T(i-1,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1) + ...
T(i-1,i)*T(1:i,i+1:end)*x(i+1:end,i);
end
q = [ eye(i)-T(1:i,1:i)*T(i,i) , -T(1:i,1:i)*T(i,i-1) ; ...
-T(1:i,1:i)*T(i-1,i) , eye(i)-T(1:i,1:i)*T(i-1,i-1) ];
z = q\[ B(1:i,i)+c ; B(1:i,i-1) + c1 ];
x(1:i,i) = z(1:i);
x(1:i,i-1) = z(i+1:end);
x(i,1:i-1) = x(1:i-1,i)';
x(i-1,1:i-2) = x(1:i-2,i-1)';
i = i - 2;
end
end
if i == 1
c = T(1,:)*(x(:,2:end)*T(1,2:end)') + T(1,1)*T(1,2:end)*x(2:end,1);
x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1));
end
x = U(:,k+1:end)*x*U(:,k+1:end)';
u = U(:,1:k);
|