Description of solve1

0001 function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
0002 % function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
0003 % Solves systems of non linear equations of several variables
0004 %
0005 % INPUTS
0006 %    func:            name of the function to be solved
0007 %    x:               guess values
0008 %    j1:              equations index for which the model is solved
0009 %    j2:              unknown variables index
0010 %    jacobian_flag=1: jacobian given by the 'func' function
0011 %    jacobian_flag=0: jacobian obtained numerically
0012 %    bad_cond_flag=1: when Jacobian is badly conditionned, use an
0013 %                     alternative formula to Newton step
0014 %    varargin:        list of arguments following bad_cond_flag
0015 %
0016 % OUTPUTS
0017 %    x:               results
0018 %    check=1:         the model can not be solved
0019 %
0020 % SPECIAL REQUIREMENTS
0021 %    none
0022 
0023 % Copyright (C) 2001-2010 Dynare Team
0024 %
0025 % This file is part of Dynare.
0026 %
0027 % Dynare is free software: you can redistribute it and/or modify
0028 % it under the terms of the GNU General Public License as published by
0029 % the Free Software Foundation, either version 3 of the License, or
0030 % (at your option) any later version.
0031 %
0032 % Dynare is distributed in the hope that it will be useful,
0033 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0034 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0035 % GNU General Public License for more details.
0036 %
0037 % You should have received a copy of the GNU General Public License
0038 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0039 
0040 global M_ options_ fjac  
0041 
0042 nn = length(j1);
0043 fjac = zeros(nn,nn) ;
0044 g = zeros(nn,1) ;
0045 
0046 tolf = options_.solve_tolf ;
0047 tolx = options_.solve_tolx;
0048 tolmin = tolx ;
0049 
0050 stpmx = 100 ;
0051 maxit = options_.solve_maxit ;
0052 
0053 check = 0 ;
0054 
0055 fvec = feval(func,x,varargin{:});
0056 fvec = fvec(j1);
0057 
0058 i = find(~isfinite(fvec));
0059 
0060 if ~isempty(i)
0061     disp(['STEADY:  numerical initial values incompatible with the following' ...
0062           ' equations'])
0063     disp(j1(i)')
0064 end
0065 
0066 f = 0.5*fvec'*fvec ;
0067 
0068 if max(abs(fvec)) < tolf
0069     return ;
0070 end
0071 
0072 stpmax = stpmx*max([sqrt(x'*x);nn]) ;
0073 first_time = 1;
0074 for its = 1:maxit
0075     if jacobian_flag
0076         [fvec,fjac] = feval(func,x,varargin{:});
0077         fvec = fvec(j1);
0078         fjac = fjac(j1,j2);
0079     else
0080         dh = max(abs(x(j2)),options_.gstep(1)*ones(nn,1))*eps^(1/3);
0081         
0082         for j = 1:nn
0083             xdh = x ;
0084             xdh(j2(j)) = xdh(j2(j))+dh(j) ;
0085             t = feval(func,xdh,varargin{:});
0086             fjac(:,j) = (t(j1) - fvec)./dh(j) ;
0087             g(j) = fvec'*fjac(:,j) ;
0088         end
0089     end
0090 
0091     g = (fvec'*fjac)';
0092     if options_.debug
0093         disp(['cond(fjac) ' num2str(cond(fjac))])
0094     end
0095     M_.unit_root = 0;
0096     if M_.unit_root
0097         if first_time
0098             first_time = 0;
0099             [q,r,e]=qr(fjac);
0100             n = sum(abs(diag(r)) < 1e-12);
0101             fvec = q'*fvec;
0102             p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
0103             disp(' ')
0104             disp('STEADY with unit roots:')
0105             disp(' ')
0106             if n > 0
0107                 disp(['   The following variable(s) kept their value given in INITVAL' ...
0108                       ' or ENDVAL'])
0109                 disp(char(e(:,end-n+1:end)'*M_.endo_names))
0110             else
0111                 disp('   STEADY can''t find any unit root!')
0112             end
0113         else
0114             [q,r]=qr(fjac*e);
0115             fvec = q'*fvec;
0116             p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
0117         end     
0118     elseif bad_cond_flag && cond(fjac) > 1/sqrt(eps)
0119         fjac
0120         fjac2=fjac'*fjac
0121         p=-(fjac2+1e6*sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec);
0122     else
0123         p = -fjac\fvec ;
0124     end
0125     xold = x ;
0126     fold = f ;
0127 
0128     [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin{:});
0129 
0130     if options_.debug
0131         disp([its f])
0132         disp([xold x])
0133     end
0134     
0135     if check > 0
0136         den = max([f;0.5*nn]) ;
0137         if max(abs(g).*max([abs(x(j2)') ones(1,nn)])')/den < tolmin
0138             return
0139         else
0140             disp (' ')
0141             disp (['SOLVE: Iteration ' num2str(its)])
0142             disp (['Spurious convergence.'])
0143             disp (x)
0144             return
0145         end
0146 
0147         if max(abs(x(j2)-xold(j2))./max([abs(x(j2)') ones(1,nn)])') < tolx
0148             disp (' ')
0149             disp (['SOLVE: Iteration ' num2str(its)])
0150             disp (['Convergence on dX.'])
0151             disp (x)
0152             return
0153         end
0154     elseif max(abs(fvec)) < tolf
0155         return
0156     end
0157 end
0158 
0159 check = 1;
0160 disp(' ')
0161 disp('SOLVE: maxit has been reached')
0162 
0163 % 01/14/01 MJ lnsearch is now a separate function
0164 % 01/16/01 MJ added varargin to function evaluation
0165 % 04/13/01 MJ added test  f < tolf !!
0166 % 05/11/01 MJ changed tests for 'check' so as to remove 'continue' which is
0167 %             an instruction which appears only in version 6
0168 
0169 
0170 
0171 
0172 
0173 
0174 
0175 
0176 
0177
solve1

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
