1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
function [x, errorflag, errorcode] = newton_solve(func, x, jacobian_flag, gstep, tolf, tolx, maxit, solve_algo, incomplete_lu_opts, varargin)
% Solves systems of non linear equations of several variables using a Newton solver, with three
% variants for the inner linear solver:
% solve_algo = 6: use a sparse LU
% solve_algo = 7: use GMRES
% solve_algo = 8: use BiCGStab
%
% INPUTS
% func: name of the function to be solved
% x: guess values
% jacobian_flag=true: jacobian given by the 'func' function
% jacobian_flag=false: jacobian obtained numerically
% gstep increment multiplier in numerical derivative
% computation
% tolf tolerance for residuals
% tolx tolerance for solution variation
% maxit maximum number of iterations
% solve_algo from options_
% incomplete_lu_opts structure of options for ilu function
% varargin: list of extra arguments to the function
%
% OUTPUTS
% x: results
% errorflag=true: the model can not be solved
% Copyright © 2024-2025 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 <https://www.gnu.org/licenses/>.
nn = length(x);
errorflag = false;
% Declaration needed for sharing with nested function
fvec = NaN(nn);
fjac = NaN(nn, nn);
function update_fvec_fjac
if jacobian_flag
[fvec, fjac] = feval(func, x, varargin{:});
else
fvec = feval(func, x, varargin{:});
dh = max(abs(x),gstep(1)*ones(nn,1))*eps^(1/3);
for j = 1:nn
xdh = x;
xdh(j) = xdh(j)+dh(j);
t = feval(func, xdh, varargin{:});
fjac(:, j) = (t - fvec) ./ dh(j);
end
end
end
update_fvec_fjac;
idInf = isinf(fvec);
idNan = isnan(fvec);
idCpx = ~isreal(fvec);
if any(idInf)
disp('newton_solve: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number:')
disp(idInf')
errorcode = 0;
errorflag = true;
return
end
if any(idNan)
disp('newton_solve: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a nan:')
disp(idNan')
errorcode = 0;
errorflag = true;
return
end
if any(idNan)
disp('newton_solve: during the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a complex number:')
disp(idCpx')
errorcode = 0;
errorflag = true;
return
end
if norm(fvec, 'Inf') < tolf
% Initial guess is a solution
errorcode = -1;
return
end
for it = 1:maxit
if solve_algo == 6
p = fjac\fvec;
else
% Should be the same options as in solve_one_boundary.m and bytecode/Interpreter.cc (static case)
[L1, U1] = ilu(fjac, incomplete_lu_opts);
if solve_algo == 7
p = gmres(fjac, fvec, [], [], [], L1, U1);
elseif solve_algo == 8
p = bicgstab(fjac, fvec, [], [], L1, U1);
else
error('newton_solve: invalid value for solve_algo')
end
end
x = x - p;
update_fvec_fjac;
if norm(fvec, 'Inf') < tolf
errorcode = 1;
return
end
end
errorflag = true;
errorcode = 2;
end
|
