From fec9486fb9ae409da6d1b78775e8cf621fb13806 Mon Sep 17 00:00:00 2001 From: sebastien Date: Fri, 21 Aug 2009 12:05:34 +0000 Subject: [PATCH] Change to option solve_algo=4 of steady: Even when the Jacobian is very badly conditioned, continue to use a Newton step git-svn-id: https://www.dynare.org/svn/dynare/trunk@2862 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/dynare_solve.m | 17 ++--- matlab/solve1.m | 13 ++-- matlab/solve2.m | 171 ------------------------------------------ 3 files changed, 13 insertions(+), 188 deletions(-) delete mode 100644 matlab/solve2.m diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m index 2619c4cc2..2dfc5735b 100644 --- a/matlab/dynare_solve.m +++ b/matlab/dynare_solve.m @@ -61,7 +61,7 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin) end elseif options_.solve_algo == 1 nn = size(x,1); - [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:}); + [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,varargin{:}); elseif options_.solve_algo == 2 || options_.solve_algo == 4 nn = size(x,1) ; tolf = options_.solve_tolf ; @@ -101,27 +101,22 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin) if options_.debug disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r))]); end + + % Activate bad conditioning flag for solve_algo = 2, but not for solve_algo = 4 + bad_cond_flag = (options_.solve_algo == 2); for i=length(r)-1:-1:1 if options_.debug disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]); end - if options_.solve_algo == 2 - [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:}); - else % solve_algo=4 - [x,info]=solve2(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:}); - end + [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, bad_cond_flag, varargin{:}); if info return end end fvec = feval(func,x,varargin{:}); if max(abs(fvec)) > tolf - if options_.solve_algo == 2 - [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:}); - else % solve_algo=4 - [x,info]=solve2(func,x,1:nn,1:nn,jacobian_flag,varargin{:}); - end + [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, bad_cond_flag, varargin{:}); end elseif options_.solve_algo == 3 if jacobian_flag diff --git a/matlab/solve1.m b/matlab/solve1.m index 9f6f84cbe..e451eccbf 100644 --- a/matlab/solve1.m +++ b/matlab/solve1.m @@ -1,5 +1,5 @@ -function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin) -% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin) +function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin) +% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin) % Solves systems of non linear equations of several variables % % INPUTS @@ -9,7 +9,9 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin) % j2: unknown variables index % jacobian_flag=1: jacobian given by the 'func' function % jacobian_flag=0: jacobian obtained numerically -% varargin: list of arguments following jacobian_flag +% bad_cond_flag=1: when Jacobian is badly conditionned, use an +% alternative formula to Newton step +% varargin: list of arguments following bad_cond_flag % % OUTPUTS % x: results @@ -18,7 +20,7 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin) % SPECIAL REQUIREMENTS % none -% Copyright (C) 2001-2008 Dynare Team +% Copyright (C) 2001-2009 Dynare Team % % This file is part of Dynare. % @@ -113,8 +115,7 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin) fvec = q'*fvec; p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)]; end -% elseif cond(fjac) > 10*sqrt(eps) - elseif cond(fjac) > 1/sqrt(eps) + elseif bad_cond_flag && cond(fjac) > 1/sqrt(eps) fjac2=fjac'*fjac; p=-(fjac2+sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec); else diff --git a/matlab/solve2.m b/matlab/solve2.m deleted file mode 100644 index db0012fcd..000000000 --- a/matlab/solve2.m +++ /dev/null @@ -1,171 +0,0 @@ -function [x,check] = solve2(func,x,j1,j2,jacobian_flag,varargin) -% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin) -% Solves systems of non linear equations of several variables -% -% This solver is used for solve_algo = ... -% -% This is a modified version of solve1: -% In solve1, before proceeding to the Newton step, we first check -% that the condition number is reasonable, and we use an alternate -% step formula if it is not the case. -% Here, we first do the Newton step, then we check if the left -% matrix division returned a warning (in case of badly scale or -% nearly singular jacobian) in which case we use the alternate -% step formula. -% -% INPUTS -% func: name of the function to be solved -% x: guess values -% j1: equations index for which the model is solved -% j2: unknown variables index -% jacobian_flag=1: jacobian given by the 'func' function -% jacobian_flag=0: jacobian obtained numerically -% varargin: list of arguments following jacobian_flag -% -% OUTPUTS -% x: results -% check=1: the model can not be solved -% -% SPECIAL REQUIREMENTS -% none - -% Copyright (C) 2001-2008 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 . - - global M_ options_ fjac - - nn = length(j1); - fjac = zeros(nn,nn) ; - g = zeros(nn,1) ; - - tolf = options_.solve_tolf ; - tolx = options_.solve_tolx; - tolmin = tolx ; - - stpmx = 100 ; - maxit = options_.solve_maxit ; - - check = 0 ; - - fvec = feval(func,x,varargin{:}); - fvec = fvec(j1); - - i = find(~isfinite(fvec)); - - if ~isempty(i) - disp(['STEADY: numerical initial values incompatible with the following' ... - ' equations']) - disp(j1(i)') - end - - f = 0.5*fvec'*fvec ; - - if max(abs(fvec)) < tolf - return ; - end - - stpmax = stpmx*max([sqrt(x'*x);nn]) ; - first_time = 1; - for its = 1:maxit - if jacobian_flag - [fvec,fjac] = feval(func,x,varargin{:}); - fvec = fvec(j1); - fjac = fjac(j1,j2); - else - dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3); - - for j = 1:nn - xdh = x ; - xdh(j2(j)) = xdh(j2(j))+dh(j) ; - t = feval(func,xdh,varargin{:}); - fjac(:,j) = (t(j1) - fvec)./dh(j) ; - g(j) = fvec'*fjac(:,j) ; - end - end - - g = (fvec'*fjac)'; - if options_.debug - disp(['cond(fjac) ' num2str(cond(fjac))]) - end - M_.unit_root = 0; - if M_.unit_root - if first_time - first_time = 0; - [q,r,e]=qr(fjac); - n = sum(abs(diag(r)) < 1e-12); - fvec = q'*fvec; - p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)]; - disp(' ') - disp('STEADY with unit roots:') - disp(' ') - if n > 0 - disp([' The following variable(s) kept their value given in INITVAL' ... - ' or ENDVAL']) - disp(char(e(:,end-n+1:end)'*M_.endo_names)) - else - disp(' STEADY can''t find any unit root!') - end - else - [q,r]=qr(fjac*e); - fvec = q'*fvec; - p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)]; - end - else - lastwarn(''); - p = -fjac\fvec; - if ~isempty(lastwarn) - fjac2=fjac'*fjac; - p=-(fjac2+sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec); - end - end - xold = x ; - fold = f ; - - [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin{:}); - - if options_.debug - disp([its f]) - disp([xold x]) - end - - if check > 0 - den = max([f;0.5*nn]) ; - if max(abs(g).*max([abs(x(j2)') ones(1,nn)])')/den < tolmin - return - else - disp (' ') - disp (['SOLVE: Iteration ' num2str(its)]) - disp (['Spurious convergence.']) - disp (x) - return - end - - if max(abs(x(j2)-xold(j2))./max([abs(x(j2)') ones(1,nn)])') < tolx - disp (' ') - disp (['SOLVE: Iteration ' num2str(its)]) - disp (['Convergence on dX.']) - disp (x) - return - end - elseif max(abs(fvec)) < tolf - return - end - end - - check = 1; - disp(' ') - disp('SOLVE: maxit has been reached')