0001 function t_nleqs_master(quiet)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 if nargin < 1
0013 quiet = 0;
0014 end
0015
0016 core_sp_newton = struct( ...
0017 'alg', 'NEWTON', ...
0018 'name', 'Newton''s', ...
0019 'default_max_it', 10, ...
0020 'need_jac', 1, ...
0021 'update_fcn', @(x, f, J)newton_update_fcn(x, f, J, '') );
0022 core_sp_gs = struct( ...
0023 'alg', 'GS', ...
0024 'name', 'Gauss-Seidel', ...
0025 'default_max_it', 1000, ...
0026 'need_jac', 0, ...
0027 'update_fcn', @(x, f)x_update_fcn2(x, f) );
0028
0029 if have_feature('matlab')
0030
0031 cfg = {
0032 {'DEFAULT', 'default', [] [] },
0033 {'NEWTON', 'Newton', [], [] },
0034 {'FD', 'fast-decoupled Newton',[],[] },
0035 {'FSOLVE', 'fsolve-1', 'fsolve', [] },
0036 {'FSOLVE', 'fsolve-2', 'fsolve', struct('Algorithm', 'trust-region-dogleg') },
0037 {'FSOLVE', 'fsolve-3', 'fsolve', struct('Algorithm', 'trust-region-reflective') },
0038 {'FSOLVE', 'fsolve-4', 'fsolve', struct('Algorithm', 'levenberg-marquardt', 'TolX', 1e-11) },
0039 {'GS', 'Gauss-Seidel',[], [] },
0040 {'CORE-N', 'Newton-CORE', [], core_sp_newton },
0041 {'CORE-GS', 'Gauss-Seidel-CORE',[], core_sp_gs },
0042 };
0043 if have_feature('matlab', 'vnum') <= 7.010
0044 cfg([6]) = [];
0045 end
0046 else
0047
0048 cfg = {
0049 {'DEFAULT', 'default', [] [] },
0050 {'NEWTON', 'Newton', [], [] },
0051 {'FD', 'fast-decoupled Newton',[],[] },
0052 {'FSOLVE', 'fsolve', 'fsolve', struct('TolX', 1e-11) },
0053 {'GS', 'Gauss-Seidel',[], [] },
0054 {'CORE-N', 'Newton-CORE', [], core_sp_newton },
0055 {'CORE-GS', 'Gauss-Seidel-CORE',[], core_sp_gs },
0056 };
0057 end
0058
0059 n = 15;
0060
0061 t_begin(n*length(cfg), quiet);
0062
0063 for k = 1:length(cfg)
0064 alg = cfg{k}{1};
0065 name = cfg{k}{2};
0066 check = cfg{k}{3};
0067 opts = cfg{k}{4};
0068 if ~isempty(check) && ~have_feature(check)
0069 t_skip(n, sprintf('%s not installed', name));
0070 else
0071 opt = struct('verbose', 0, 'alg', alg, 'tol', 1e-11);
0072 switch alg
0073 case {'DEFAULT', 'NEWTON'}
0074 case {'FD'}
0075 opt.fd_opt.jac_approx_fcn = @jac_approx_fcn2;
0076 case 'FSOLVE'
0077 opt.fsolve_opt = opts;
0078 case 'GS'
0079 opt.gs_opt.x_update_fcn = @(x, f)x_update_fcn2(x, f);
0080 case {'CORE-N', 'CORE-GS'}
0081 opt.core_sp = opts;
0082 opt.alg = 'CORE';
0083 end
0084
0085 switch alg
0086 case {'DEFAULT', 'NEWTON', 'FSOLVE', 'CORE-N'}
0087 t = sprintf('%s - 2-d function : ', name);
0088 x0 = [-1;0];
0089 [x, f, e, out, jac] = nleqs_master(@f1, x0, opt);
0090 t_is(e, 1, 12, [t 'success']);
0091 t_is(x, [-3; 4], 8, [t 'x']);
0092 t_is(f, 0, 10, [t 'f']);
0093 switch alg
0094 case {'DEFAULT', 'CORE-N'}
0095 out_alg = 'NEWTON';
0096 otherwise
0097 out_alg = alg;
0098 end
0099 t_ok(strcmp(out.alg, out_alg), [t 'out.alg']);
0100 eJ = [1 1; 6 1];
0101 t_is(jac, eJ, 5.8, [t 'jac']);
0102
0103 t = sprintf('%s - 2-d function (struct) : ', name);
0104 p = struct('fcn', @f1, 'x0', [1;0], 'opt', opt);
0105 [x, f, e] = nleqs_master(p);
0106 t_is(e, 1, 12, [t 'success']);
0107 t_is(x, [2; -1], 8, [t 'x']);
0108 t_is(f, 0, 10, [t 'f']);
0109
0110 p.opt.max_it = 3;
0111 t = sprintf('%s - 2-d function (max_it) : ', name);
0112 [x, f, e, out] = nleqs_master(p);
0113 t_is(e, 0, 12, [t 'no success']);
0114 t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
0115 otherwise
0116 t_skip(10, sprintf('not implemented for solver ''%s''', alg));
0117 end
0118
0119 t = sprintf('%s - 2-d function2 (struct) : ', name);
0120 p = struct('fcn', @f2, 'x0', [1;2], 'opt', opt);
0121 [x, f, e] = nleqs_master(p);
0122 t_is(e, 1, 12, [t 'success']);
0123 t_is(x, [2; 3], 8, [t 'x']);
0124 t_is(f, 0, 10, [t 'f']);
0125
0126 p.opt.max_it = 3;
0127 t = sprintf('%s - 2-d function2 (max_it) : ', name);
0128 [x, f, e, out] = nleqs_master(p);
0129 t_is(e, 0, 12, [t 'no success']);
0130 t_ok(out.iterations == 3 || out.iterations == 4, [t 'iterations']);
0131 end
0132 end
0133
0134 t_end;
0135
0136
0137
0138
0139 function [f, J] = f1(x)
0140 f = [ x(1) + x(2) - 1;
0141 -x(1)^2 + x(2) + 5 ];
0142 if nargout > 1
0143 J = sparse([1 1; -2*x(1) 1]);
0144 end
0145
0146
0147
0148 function [f, J] = f2(x)
0149 f = [ x(1)^2 + x(1)*x(2) - 10;
0150 x(2) + 3*x(1)*x(2)^2 - 57 ];
0151 if nargout > 1
0152 J = sparse([ 2*x(1)+x(2) x(1);
0153 3*x(2)^2 6*x(1)*x(2)+1 ]);
0154 end
0155
0156 function JJ = jac_approx_fcn2()
0157
0158 J = [7 2; 27 37];
0159 JJ = {J(1,1), J(2,2)};
0160
0161 function x = x_update_fcn2(x, f)
0162
0163 x(1) = sqrt(10 - x(1)*x(2));
0164 x(2) = sqrt((57-x(2))/3/x(1));
0165
0166 function x = newton_update_fcn(x, f, J, lin_solver)
0167 dx = mplinsolve(J, -f, lin_solver);
0168 x = x + dx;