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
%
% OUTPUTS
% x: solution
% info=1: the model can not be solved
%
% SPECIAL REQUIREMENTS
% none
2012-01-09 11:06:26 +01:00
% Copyright (C) 2001-2012 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' )
for ii = 1 : length ( infrow )
if infcol ( ii ) < = M . orig_endo_nbr
2013-07-18 15:00:41 +02:00
fprintf ( ' STEADY: Derivative of Equation %d with respect to Variable %s (initial value of %s: %g) \n' , infrow ( ii ) , deblank ( M . endo_names ( infcol ( ii ) , : ) ) , deblank ( M . endo_names ( infcol ( ii ) , : ) ) , x ( infcol ( ii ) ) )
2013-06-20 18:53:35 +02:00
else %auxiliary vars
orig_var_index = M . aux_vars ( 1 , infcol ( ii ) - M . orig_endo_nbr ) . orig_index ;
2013-07-18 15:00:41 +02:00
fprintf ( ' STEADY: Derivative of Equation %d with respect to Variable %s (initial value of %s: %g) \n' , infrow ( ii ) , deblank ( M . endo_names ( orig_var_index , : ) ) , deblank ( M . endo_names ( orig_var_index , : ) ) , x ( infcol ( ii ) ) )
2013-06-20 18:53:35 +02:00
end
end
fprintf ( ' \nSTEADY: The problem most often occurs, because a variable with\n' )
fprintf ( ' STEADY: exponent smaller than 1 has been initialized to 0. Taking the derivative\n' )
fprintf ( ' STEADY: and evaluating it at the steady state then results in a division by 0.\n' )
error ( ' An element of the Jacobian is not finite or NaN' )
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' )
disp ( ' ii) that no division by an endogenous variable initialized to 0 occurs' )
info = 1 ;
x = NaN ;
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 ;
options . MaxIter = 2000 ;
options . TolFun = 1e-8 ;
options . Display = ' iter' ;
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
2011-02-04 17:17:48 +01: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
2012-08-06 23:32:15 +02:00
[ x , info ] = solve1 ( func , x , 1 : nn , 1 : nn , jacobian_flag , 1 , 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
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 ) ;
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
% Activate bad conditioning flag for solve_algo = 2, but not for solve_algo = 4
bad_cond_flag = ( options_ . solve_algo == 2 ) ;
2008-09-16 19:24:11 +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
2012-08-06 23:32:15 +02:00
[ x , info ] = solve1 ( func , x , j1 ( r ( i ) : r ( i + 1 ) - 1 ) , j2 ( r ( i ) : r ( i + 1 ) - 1 ) , jacobian_flag , ...
bad_cond_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
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
2012-08-06 23:32:15 +02:00
[ x , info ] = solve1 ( func , x , 1 : nn , 1 : nn , jacobian_flag , bad_cond_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
2009-12-16 18:17:34 +01:00
else
2008-09-22 15:13:04 +02:00
error ( ' DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4]' )
2009-12-16 18:17:34 +01:00
end