2006-08-24 15:07:26 +02:00
|
|
|
function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
|
2008-01-03 13:12:07 +01:00
|
|
|
% function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
|
|
|
|
% proposes different solvers
|
|
|
|
%
|
|
|
|
% INPUTS
|
|
|
|
% func: name of the function to be solved
|
|
|
|
% x: guess values
|
|
|
|
% jacobian_flag=1: jacobian given by the 'func' function
|
|
|
|
% jacobian_flag=0: jacobian obtained numerically
|
|
|
|
% varargin: list of arguments following jacobian_flag
|
2014-09-10 16:30:17 +02:00
|
|
|
%
|
2008-01-03 13:12:07 +01:00
|
|
|
% OUTPUTS
|
|
|
|
% x: solution
|
|
|
|
% info=1: the model can not be solved
|
|
|
|
%
|
|
|
|
% SPECIAL REQUIREMENTS
|
|
|
|
% none
|
|
|
|
|
2014-08-31 17:52:05 +02:00
|
|
|
% Copyright (C) 2001-2014 Dynare Team
|
2008-08-01 14:40:33 +02:00
|
|
|
%
|
|
|
|
% 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/>.
|
2008-01-03 13:12:07 +01:00
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
global options_
|
|
|
|
|
2012-10-31 11:08:22 +01:00
|
|
|
tolf = options_.solve_tolf ;
|
2009-12-16 18:17:34 +01:00
|
|
|
info = 0;
|
2012-11-01 09:26:05 +01:00
|
|
|
nn = size(x,1);
|
2012-10-31 11:08:22 +01:00
|
|
|
|
|
|
|
% checking initial values
|
|
|
|
if jacobian_flag
|
|
|
|
[fvec,fjac] = feval(func,x,varargin{:});
|
2013-06-20 18:53:35 +02:00
|
|
|
if any(any(isinf(fjac) | isnan(fjac)))
|
|
|
|
[infrow,infcol]=find(isinf(fjac) | isnan(fjac));
|
|
|
|
M=evalin('base','M_'); %get variable names from workspace
|
|
|
|
fprintf('\nSTEADY: The Jacobian contains Inf or NaN. The problem arises from: \n\n')
|
2014-03-28 14:52:06 +01:00
|
|
|
display_problematic_vars_Jacobian(infrow,infcol,M,x,'static','STEADY: ')
|
2014-09-10 16:30:17 +02:00
|
|
|
error('An element of the Jacobian is not finite or NaN')
|
2013-06-20 18:53:35 +02:00
|
|
|
end
|
2012-10-31 11:08:22 +01:00
|
|
|
else
|
|
|
|
fvec = feval(func,x,varargin{:});
|
|
|
|
fjac = zeros(nn,nn) ;
|
|
|
|
end
|
|
|
|
|
|
|
|
i = find(~isfinite(fvec));
|
|
|
|
|
|
|
|
if ~isempty(i)
|
|
|
|
disp(['STEADY: numerical initial values or parameters incompatible with the following' ...
|
|
|
|
' equations'])
|
|
|
|
disp(i')
|
|
|
|
disp('Please check for example')
|
|
|
|
disp(' i) if all parameters occurring in these equations are defined')
|
2014-09-10 16:30:17 +02:00
|
|
|
disp(' ii) that no division by an endogenous variable initialized to 0 occurs')
|
2012-10-31 11:08:22 +01:00
|
|
|
info = 1;
|
2014-05-26 16:16:22 +02:00
|
|
|
x = NaN(size(fvec));
|
2012-10-31 11:08:22 +01:00
|
|
|
return;
|
|
|
|
end
|
|
|
|
|
|
|
|
if max(abs(fvec)) < tolf
|
|
|
|
return ;
|
|
|
|
end
|
|
|
|
|
2009-12-16 18:17:34 +01:00
|
|
|
if options_.solve_algo == 0
|
2013-11-04 10:54:45 +01:00
|
|
|
if ~isoctave
|
2012-01-09 11:06:26 +01:00
|
|
|
if ~user_has_matlab_license('optimization_toolbox')
|
2011-12-19 12:36:35 +01:00
|
|
|
error('You can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
|
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
end
|
2008-09-22 15:13:04 +02:00
|
|
|
options=optimset('fsolve');
|
|
|
|
options.MaxFunEvals = 50000;
|
2014-08-31 17:52:05 +02:00
|
|
|
options.MaxIter = options_.steady.maxit;
|
2014-09-10 16:30:17 +02:00
|
|
|
options.TolFun = tolf;
|
|
|
|
options.Display = 'iter';
|
2008-09-22 15:13:04 +02:00
|
|
|
if jacobian_flag
|
2009-12-16 18:17:34 +01:00
|
|
|
options.Jacobian = 'on';
|
2008-09-22 15:13:04 +02:00
|
|
|
else
|
2009-12-16 18:17:34 +01:00
|
|
|
options.Jacobian = 'off';
|
2008-09-22 15:13:04 +02:00
|
|
|
end
|
2013-11-04 10:54:45 +01:00
|
|
|
if ~isoctave
|
2010-10-28 12:07:33 +02:00
|
|
|
[x,fval,exitval,output] = fsolve(func,x,options,varargin{:});
|
|
|
|
else
|
|
|
|
% Under Octave, use a wrapper, since fsolve() does not have a 4th arg
|
|
|
|
func2 = str2func(func);
|
|
|
|
func = @(x) func2(x, varargin{:});
|
2011-09-19 12:36:37 +02:00
|
|
|
% The Octave version of fsolve does not converge when it starts from the solution
|
|
|
|
fvec = feval(func,x);
|
2012-10-31 11:08:22 +01:00
|
|
|
if max(abs(fvec)) >= tolf
|
2011-01-14 12:27:13 +01:00
|
|
|
[x,fval,exitval,output] = fsolve(func,x,options);
|
|
|
|
else
|
|
|
|
exitval = 3;
|
2011-02-04 17:17:48 +01:00
|
|
|
end;
|
2010-10-28 12:07:33 +02:00
|
|
|
end
|
2014-09-10 16:30:17 +02:00
|
|
|
|
2008-09-22 15:13:04 +02:00
|
|
|
if exitval > 0
|
2009-12-16 18:17:34 +01:00
|
|
|
info = 0;
|
2008-09-22 15:13:04 +02:00
|
|
|
else
|
2009-12-16 18:17:34 +01:00
|
|
|
info = 1;
|
2008-09-22 15:13:04 +02:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
elseif options_.solve_algo == 1
|
2014-05-07 21:25:10 +02:00
|
|
|
[x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,options_.gstep, ...
|
|
|
|
tolf,options_.solve_tolx, ...
|
|
|
|
options_.steady.maxit,options_.debug,varargin{:});
|
2014-05-09 15:33:26 +02:00
|
|
|
elseif options_.solve_algo == 9
|
2014-05-07 21:25:10 +02:00
|
|
|
[x,info]=trust_region(func,x,1:nn,1:nn,jacobian_flag,options_.gstep, ...
|
2012-10-31 11:08:22 +01:00
|
|
|
tolf,options_.solve_tolx, ...
|
2013-10-09 15:38:07 +02:00
|
|
|
options_.steady.maxit,options_.debug,varargin{:});
|
2009-12-16 18:17:34 +01:00
|
|
|
elseif options_.solve_algo == 2 || options_.solve_algo == 4
|
2014-02-04 17:55:55 +01:00
|
|
|
|
|
|
|
if options_.solve_algo == 2
|
|
|
|
solver = @solve1;
|
|
|
|
else
|
|
|
|
solver = @trust_region;
|
|
|
|
end
|
|
|
|
|
2006-01-18 17:50:33 +01:00
|
|
|
if ~jacobian_flag
|
2009-12-16 18:17:34 +01:00
|
|
|
fjac = zeros(nn,nn) ;
|
2012-03-09 16:36:26 +01:00
|
|
|
dh = max(abs(x),options_.gstep(1)*ones(nn,1))*eps^(1/3);
|
2009-12-16 18:17:34 +01:00
|
|
|
for j = 1:nn
|
|
|
|
xdh = x ;
|
|
|
|
xdh(j) = xdh(j)+dh(j) ;
|
|
|
|
fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ;
|
|
|
|
end
|
2006-01-18 17:50:33 +01:00
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
|
|
|
|
[j1,j2,r,s] = dmperm(fjac);
|
2014-09-10 16:30:17 +02:00
|
|
|
|
2008-09-16 19:24:11 +02:00
|
|
|
if options_.debug
|
2009-12-16 18:17:34 +01:00
|
|
|
disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r))]);
|
2008-09-16 19:24:11 +02:00
|
|
|
end
|
2009-08-21 14:05:34 +02:00
|
|
|
|
2005-02-18 20:54:39 +01:00
|
|
|
for i=length(r)-1:-1:1
|
2009-12-16 18:17:34 +01:00
|
|
|
if options_.debug
|
|
|
|
disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
|
|
|
|
end
|
2014-02-04 17:55:55 +01:00
|
|
|
[x,info]=solver(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, ...
|
2014-02-04 17:46:54 +01:00
|
|
|
options_.gstep, ...
|
2012-10-31 11:08:22 +01:00
|
|
|
tolf,options_.solve_tolx, ...
|
2013-10-09 15:38:07 +02:00
|
|
|
options_.steady.maxit,options_.debug,varargin{:});
|
2009-12-16 18:17:34 +01:00
|
|
|
if info
|
|
|
|
return
|
|
|
|
end
|
2005-02-18 20:54:39 +01:00
|
|
|
end
|
2007-01-04 15:42:27 +01:00
|
|
|
fvec = feval(func,x,varargin{:});
|
|
|
|
if max(abs(fvec)) > tolf
|
2014-02-04 17:55:55 +01:00
|
|
|
[x,info]=solver(func,x,1:nn,1:nn,jacobian_flag, ...
|
2012-10-31 11:08:22 +01:00
|
|
|
options_.gstep, tolf,options_.solve_tolx, ...
|
2013-10-09 15:38:07 +02:00
|
|
|
options_.steady.maxit,options_.debug,varargin{:});
|
2007-01-04 15:42:27 +01:00
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
elseif options_.solve_algo == 3
|
2007-11-27 15:16:13 +01:00
|
|
|
if jacobian_flag
|
2009-12-16 18:17:34 +01:00
|
|
|
[x,info] = csolve(func,x,func,1e-6,500,varargin{:});
|
2007-11-27 15:16:13 +01:00
|
|
|
else
|
2009-12-16 18:17:34 +01:00
|
|
|
[x,info] = csolve(func,x,[],1e-6,500,varargin{:});
|
2008-09-22 15:13:04 +02:00
|
|
|
end
|
2014-05-12 17:15:49 +02:00
|
|
|
elseif options_.solve_algo == 10
|
|
|
|
olmmcp = options_.lmmcp;
|
|
|
|
[x,fval,exitflag] = lmmcp(func,x,olmmcp.lb,olmmcp.ub,olmmcp,varargin{:});
|
|
|
|
if exitflag == 1
|
|
|
|
info = 0;
|
|
|
|
else
|
|
|
|
info = 1;
|
|
|
|
end
|
2009-12-16 18:17:34 +01:00
|
|
|
else
|
2014-05-09 15:33:26 +02:00
|
|
|
error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9]')
|
2009-12-16 18:17:34 +01:00
|
|
|
end
|