dynare_solve

 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

0001 function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
0002 % function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
0003 % proposes different solvers
0004 %
0005 % INPUTS
0006 %    func:             name of the function to be solved
0007 %    x:                guess values
0008 %    jacobian_flag=1:  jacobian given by the 'func' function
0009 %    jacobian_flag=0:  jacobian obtained numerically
0010 %    varargin:         list of arguments following jacobian_flag
0011 %
0012 % OUTPUTS
0013 %    x:                solution
0014 %    info=1:           the model can not be solved
0015 %
0016 % SPECIAL REQUIREMENTS
0017 %    none
0018 
0019 % Copyright (C) 2001-2012 Dynare Team
0020 %
0021 % This file is part of Dynare.
0022 %
0023 % Dynare is free software: you can redistribute it and/or modify
0024 % it under the terms of the GNU General Public License as published by
0025 % the Free Software Foundation, either version 3 of the License, or
0026 % (at your option) any later version.
0027 %
0028 % Dynare is distributed in the hope that it will be useful,
0029 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0030 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0031 % GNU General Public License for more details.
0032 %
0033 % You should have received a copy of the GNU General Public License
0034 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0035 
0036 global options_
0037 
0038 info = 0;
0039 if options_.solve_algo == 0
0040     if ~exist('OCTAVE_VERSION')
0041         if ~user_has_matlab_license('optimization_toolbox')
0042             error('You can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
0043         end
0044     end
0045     options=optimset('fsolve');
0046     options.MaxFunEvals = 50000;
0047     options.MaxIter = 2000;
0048     options.TolFun=1e-8;
0049     options.Display = 'iter';
0050     if jacobian_flag
0051         options.Jacobian = 'on';
0052     else
0053         options.Jacobian = 'off';
0054     end
0055     if ~exist('OCTAVE_VERSION')
0056         [x,fval,exitval,output] = fsolve(func,x,options,varargin{:});
0057     else
0058         % Under Octave, use a wrapper, since fsolve() does not have a 4th arg
0059         func2 = str2func(func);
0060         func = @(x) func2(x, varargin{:});
0061         % The Octave version of fsolve does not converge when it starts from the solution
0062         fvec = feval(func,x);
0063         if max(abs(fvec)) >= options_.solve_tolf
0064             [x,fval,exitval,output] = fsolve(func,x,options);
0065         else
0066             exitval = 3;
0067         end;
0068     end
0069     
0070     if exitval > 0
0071         info = 0;
0072     else
0073         info = 1;
0074     end
0075 elseif options_.solve_algo == 1
0076     nn = size(x,1);
0077     [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,varargin{:});
0078 elseif options_.solve_algo == 2 || options_.solve_algo == 4
0079     nn = size(x,1) ;
0080     tolf = options_.solve_tolf ;
0081 
0082     if jacobian_flag
0083         [fvec,fjac] = feval(func,x,varargin{:});
0084     else
0085         fvec = feval(func,x,varargin{:});
0086         fjac = zeros(nn,nn) ;
0087     end
0088 
0089     i = find(~isfinite(fvec));
0090     
0091     if ~isempty(i)
0092         disp(['STEADY:  numerical initial values incompatible with the following' ...
0093               ' equations'])
0094         disp(i')
0095         disp('Please check for example')
0096         disp('   i) if all parameters occurring in these equations are defined')
0097         disp('  ii) that no division by an endogenous variable initialized to 0 occurs') 
0098         error('exiting ...')
0099     end
0100     
0101     if max(abs(fvec)) < tolf
0102         return ;
0103     end
0104 
0105     if ~jacobian_flag
0106         fjac = zeros(nn,nn) ;
0107         dh = max(abs(x),options_.gstep(1)*ones(nn,1))*eps^(1/3);
0108         for j = 1:nn
0109             xdh = x ;
0110             xdh(j) = xdh(j)+dh(j) ;
0111             fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ;
0112         end
0113     end
0114 
0115     [j1,j2,r,s] = dmperm(fjac);
0116     
0117     if options_.debug
0118         disp(['DYNARE_SOLVE (solve_algo=2|4): number of blocks = ' num2str(length(r))]);
0119     end
0120 
0121     % Activate bad conditioning flag for solve_algo = 2, but not for solve_algo = 4
0122     bad_cond_flag = (options_.solve_algo == 2);
0123     
0124     for i=length(r)-1:-1:1
0125         if options_.debug
0126             disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
0127         end
0128         [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{:});
0129         if info
0130             return
0131         end
0132     end
0133     fvec = feval(func,x,varargin{:});
0134     if max(abs(fvec)) > tolf
0135         [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, bad_cond_flag, varargin{:});
0136     end
0137 elseif options_.solve_algo == 3
0138     if jacobian_flag
0139         [x,info] = csolve(func,x,func,1e-6,500,varargin{:});
0140     else
0141         [x,info] = csolve(func,x,[],1e-6,500,varargin{:});
0142     end
0143 else
0144     error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4]')
0145 end