From 2d83b2216d7b9518718efd83808a1d7c47c488dd Mon Sep 17 00:00:00 2001
From: sebastien <sebastien@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Tue, 23 Sep 2008 09:06:03 +0000
Subject: [PATCH] 4.0: merged changesets 2078, 2089 and 2090 from trunk * new
debug messages in dynare_solve.m * solve_algo=2 exits immediately after
failure to solve one block * new solve_algo=4 with different way of dealing
with badly scaled or nearly singular Jacobian * solve_algo=0 fails if
Matlab's Optimization Toolbox not present * dynare_solve fails if solve_algo
not in [0,1,2,3,4]
git-svn-id: https://www.dynare.org/svn/dynare/branches/4.0@2095 ac1d8469-bf42-47a9-8791-bf33cf982152
---
matlab/dynare_solve.m | 80 +++++++++++---------
matlab/solve2.m | 171 ++++++++++++++++++++++++++++++++++++++++++
2 files changed, 215 insertions(+), 36 deletions(-)
create mode 100644 matlab/solve2.m
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index c376a4f1aa..2619c4cc25 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -38,33 +38,31 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
options_ = set_default_option(options_,'solve_algo',2);
info = 0;
if options_.solve_algo == 0
- if ~isempty(which('fsolve'))
- options=optimset('fsolve');
- options.MaxFunEvals = 50000;
- options.MaxIter = 2000;
- options.TolFun=1e-8;
- options.Display = 'iter';
- if jacobian_flag
- options.Jacobian = 'on';
- else
- options.Jacobian = 'off';
- end
- [x,fval,exitval,output] = fsolve(func,x,options,varargin{:});
- if exitval > 0
- info = 0;
- else
- info = 1;
- end
- return
- else
- options_.solve_algo = 1;
+ if exist('OCTAVE_VERSION') || isempty(ver('optim'))
+ % Note that fsolve() exists under Octave, but has a different syntax
+ % So we fail for the moment under Octave, until we add the corresponding code
+ error('DYNARE_SOLVE: you can''t use solve_algo=0 since you don''t have Matlab''s Optimization Toolbox')
end
- end
-
- if options_.solve_algo == 1
+ options=optimset('fsolve');
+ options.MaxFunEvals = 50000;
+ options.MaxIter = 2000;
+ options.TolFun=1e-8;
+ options.Display = 'iter';
+ if jacobian_flag
+ options.Jacobian = 'on';
+ else
+ options.Jacobian = 'off';
+ end
+ [x,fval,exitval,output] = fsolve(func,x,options,varargin{:});
+ if exitval > 0
+ info = 0;
+ else
+ info = 1;
+ end
+ elseif options_.solve_algo == 1
nn = size(x,1);
[x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
- elseif options_.solve_algo == 2
+ elseif options_.solve_algo == 2 || options_.solve_algo == 4
nn = size(x,1) ;
tolf = options_.solve_tolf ;
@@ -84,8 +82,6 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
error('exiting ...')
end
-% f = 0.5*fvec'*fvec ;
-
if max(abs(fvec)) < tolf
return ;
end
@@ -102,25 +98,37 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
[j1,j2,r,s] = dmperm(fjac);
+ if options_.debug
+ disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r))]);
+ end
+
for i=length(r)-1:-1:1
- [x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:});
- if info & options_.debug
- error(sprintf('Solve block = %d check = %d\n',i,info));
+ 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
+ if info
+ return
end
end
fvec = feval(func,x,varargin{:});
if max(abs(fvec)) > tolf
- [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
+ 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
end
elseif options_.solve_algo == 3
if jacobian_flag
[x,info] = csolve(func,x,func,1e-6,500,varargin{:});
else
[x,info] = csolve(func,x,[],1e-6,500,varargin{:});
- end
+ end
+ else
+ error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4]')
end
-% fvec1 = feval(func,x,varargin{:})
-
- % 08/28/03 MJ add a final call to solve1 for solve_algo == 1 in case
- % initvals generates 'false' zeros in the Jacobian
-
\ No newline at end of file
diff --git a/matlab/solve2.m b/matlab/solve2.m
new file mode 100644
index 0000000000..db0012fcd6
--- /dev/null
+++ b/matlab/solve2.m
@@ -0,0 +1,171 @@
+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 <http://www.gnu.org/licenses/>.
+
+ 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')
--
GitLab