Mercurial > hg > octave-lyh

comparison scripts/optimization/fsolve.m @ 11587:c792872f8942

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
all script files: untabify and strip trailing whitespace
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:35:29 -0500
parents fd0a3ac60b0e
children b9a89ca0fb75
comparison
equal deleted inserted replaced
11586:12df7854fa7c 11587:c792872f8942
23 ## @deftypefnx {Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}] =} fsolve (@var{fcn}, @dots{}) 23 ## @deftypefnx {Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{fjac}] =} fsolve (@var{fcn}, @dots{})
24 ## Solve a system of nonlinear equations defined by the function @var{fcn}. 24 ## Solve a system of nonlinear equations defined by the function @var{fcn}.
25 ## @var{fcn} should accept a vector (array) defining the unknown variables, 25 ## @var{fcn} should accept a vector (array) defining the unknown variables,
26 ## and return a vector of left-hand sides of the equations. Right-hand sides 26 ## and return a vector of left-hand sides of the equations. Right-hand sides
27 ## are defined to be zeros. 27 ## are defined to be zeros.
28 ## In other words, this function attempts to determine a vector @var{x} such 28 ## In other words, this function attempts to determine a vector @var{x} such
29 ## that @code{@var{fcn} (@var{x})} gives (approximately) all zeros. 29 ## that @code{@var{fcn} (@var{x})} gives (approximately) all zeros.
30 ## @var{x0} determines a starting guess. The shape of @var{x0} is preserved 30 ## @var{x0} determines a starting guess. The shape of @var{x0} is preserved
31 ## in all calls to @var{fcn}, but otherwise it is treated as a column vector. 31 ## in all calls to @var{fcn}, but otherwise it is treated as a column vector.
32 ## @var{options} is a structure specifying additional options. 32 ## @var{options} is a structure specifying additional options.
33 ## Currently, @code{fsolve} recognizes these options: 33 ## Currently, @code{fsolve} recognizes these options:
34 ## @code{"FunValCheck"}, @code{"OutputFcn"}, @code{"TolX"}, 34 ## @code{"FunValCheck"}, @code{"OutputFcn"}, @code{"TolX"},
35 ## @code{"TolFun"}, @code{"MaxIter"}, @code{"MaxFunEvals"}, 35 ## @code{"TolFun"}, @code{"MaxIter"}, @code{"MaxFunEvals"},
36 ## @code{"Jacobian"}, @code{"Updating"}, @code{"ComplexEqn"} 36 ## @code{"Jacobian"}, @code{"Updating"}, @code{"ComplexEqn"}
37 ## @code{"TypicalX"}, @code{"AutoScaling"} and @code{"FinDiffType"}. 37 ## @code{"TypicalX"}, @code{"AutoScaling"} and @code{"FinDiffType"}.
38 ## 38 ##
39 ## If @code{"Jacobian"} is @code{"on"}, it specifies that @var{fcn}, 39 ## If @code{"Jacobian"} is @code{"on"}, it specifies that @var{fcn},
40 ## called with 2 output arguments, also returns the Jacobian matrix 40 ## called with 2 output arguments, also returns the Jacobian matrix
41 ## of right-hand sides at the requested point. @code{"TolX"} specifies 41 ## of right-hand sides at the requested point. @code{"TolX"} specifies
42 ## the termination tolerance in the unknown variables, while 42 ## the termination tolerance in the unknown variables, while
43 ## @code{"TolFun"} is a tolerance for equations. Default is @code{1e-7} 43 ## @code{"TolFun"} is a tolerance for equations. Default is @code{1e-7}
44 ## for both @code{"TolX"} and @code{"TolFun"}. 44 ## for both @code{"TolX"} and @code{"TolFun"}.
45 ## 45 ##
46 ## If @code{"AutoScaling"} is on, the variables will be automatically scaled 46 ## If @code{"AutoScaling"} is on, the variables will be automatically scaled
47 ## according to the column norms of the (estimated) Jacobian. As a result, 47 ## according to the column norms of the (estimated) Jacobian. As a result,
48 ## TolF becomes scaling-independent. By default, this option is off, because 48 ## TolF becomes scaling-independent. By default, this option is off, because
49 ## it may sometimes deliver unexpected (though mathematically correct) results. 49 ## it may sometimes deliver unexpected (though mathematically correct) results.
50 ## 50 ##
51 ## If @code{"Updating"} is "on", the function will attempt to use Broyden 51 ## If @code{"Updating"} is "on", the function will attempt to use Broyden
52 ## updates to update the Jacobian, in order to reduce the amount of Jacobian 52 ## updates to update the Jacobian, in order to reduce the amount of Jacobian
53 ## calculations. 53 ## calculations.
54 ## If your user function always calculates the Jacobian (regardless of number 54 ## If your user function always calculates the Jacobian (regardless of number
55 ## of output arguments), this option provides no advantage and should be set to 55 ## of output arguments), this option provides no advantage and should be set to
56 ## false. 56 ## false.
57 ## 57 ##
58 ## @code{"ComplexEqn"} is @code{"on"}, @code{fsolve} will attempt to solve 58 ## @code{"ComplexEqn"} is @code{"on"}, @code{fsolve} will attempt to solve
59 ## complex equations in complex variables, assuming that the equations possess a 59 ## complex equations in complex variables, assuming that the equations possess a
60 ## complex derivative (i.e., are holomorphic). If this is not what you want, 60 ## complex derivative (i.e., are holomorphic). If this is not what you want,
61 ## should unpack the real and imaginary parts of the system to get a real 61 ## should unpack the real and imaginary parts of the system to get a real
62 ## system. 62 ## system.
63 ## 63 ##
64 ## For description of the other options, see @code{optimset}. 64 ## For description of the other options, see @code{optimset}.
65 ## 65 ##
66 ## On return, @var{fval} contains the value of the function @var{fcn} 66 ## On return, @var{fval} contains the value of the function @var{fcn}
67 ## evaluated at @var{x}, and @var{info} may be one of the following values: 67 ## evaluated at @var{x}, and @var{info} may be one of the following values:
68 ## 68 ##
69 ## @table @asis 69 ## @table @asis
70 ## @item 1 70 ## @item 1
71 ## Converged to a solution point. Relative residual error is less than 71 ## Converged to a solution point. Relative residual error is less than
72 ## specified by TolFun. 72 ## specified by TolFun.
73 ## 73 ##
74 ## @item 2 74 ## @item 2
75 ## Last relative step size was less that TolX. 75 ## Last relative step size was less that TolX.
76 ## 76 ##
77 ## @item 3 77 ## @item 3
78 ## Last relative decrease in residual was less than TolF. 78 ## Last relative decrease in residual was less than TolF.
79 ## 79 ##
80 ## @item 0 80 ## @item 0
81 ## Iteration limit exceeded. 81 ## Iteration limit exceeded.
82 ## 82 ##
83 ## @item -3 83 ## @item -3
84 ## The trust region radius became excessively small. 84 ## The trust region radius became excessively small.
85 ## @end table 85 ## @end table
86 ## 86 ##
87 ## Note: If you only have a single nonlinear equation of one variable, using 87 ## Note: If you only have a single nonlinear equation of one variable, using
88 ## @code{fzero} is usually a much better idea. 88 ## @code{fzero} is usually a much better idea.
89 ## @seealso{fzero, optimset} 89 ## @seealso{fzero, optimset}
90 ## 90 ##
91 ## Note about user-supplied Jacobians: 91 ## Note about user-supplied Jacobians:
92 ## As an inherent property of the algorithm, Jacobian is always requested for a 92 ## As an inherent property of the algorithm, Jacobian is always requested for a
120 ## ## maybe output iteration status, etc. 120 ## ## maybe output iteration status, etc.
121 ## endif 121 ## endif
122 ## endfunction 122 ## endfunction
123 ## 123 ##
124 ## ## @dots{}. 124 ## ## @dots{}.
125 ## 125 ##
126 ## fsolve (@@user_func, x0, optimset ("OutputFcn", @@user_func, @dots{})) 126 ## fsolve (@@user_func, x0, optimset ("OutputFcn", @@user_func, @dots{}))
127 ## @end example 127 ## @end example
128 ## @end deftypefn 128 ## @end deftypefn
129 129
130 ## PKG_ADD: __all_opts__ ("fsolve"); 130 ## PKG_ADD: __all_opts__ ("fsolve");
141 return; 141 return;
142 endif 142 endif
143 143
144 if (nargin < 2 || nargin > 3 || ! ismatrix (x0)) 144 if (nargin < 2 || nargin > 3 || ! ismatrix (x0))
145 print_usage (); 145 print_usage ();
146 endif 146 endif
147 147
148 if (ischar (fcn)) 148 if (ischar (fcn))
149 fcn = str2func (fcn, "global"); 149 fcn = str2func (fcn, "global");
150 elseif (iscell (fcn)) 150 elseif (iscell (fcn))
151 fcn = @(x) make_fcn_jac (x, fcn{1}, fcn{2}); 151 fcn = @(x) make_fcn_jac (x, fcn{1}, fcn{2});
373 endif 373 endif
374 endif 374 endif
375 375
376 ## Tests for termination conditions. A mysterious place, anything 376 ## Tests for termination conditions. A mysterious place, anything
377 ## can happen if you change something here... 377 ## can happen if you change something here...
378 378
379 ## The rule of thumb (which I'm not sure M*b is quite following) 379 ## The rule of thumb (which I'm not sure M*b is quite following)
380 ## is that for a tolerance that depends on scaling, only 0 makes 380 ## is that for a tolerance that depends on scaling, only 0 makes
381 ## sense as a default value. But 0 usually means uselessly long 381 ## sense as a default value. But 0 usually means uselessly long
382 ## iterations, so we need scaling-independent tolerances wherever 382 ## iterations, so we need scaling-independent tolerances wherever
383 ## possible. 383 ## possible.
406 break; 406 break;
407 endif 407 endif
408 408
409 ## Compute the scaled Broyden update. 409 ## Compute the scaled Broyden update.
410 if (useqr) 410 if (useqr)
411 u = (fvec1 - q*w) / sn; 411 u = (fvec1 - q*w) / sn;
412 v = dg .* ((dg .* s) / sn); 412 v = dg .* ((dg .* s) / sn);
413 413
414 ## Update the QR factorization. 414 ## Update the QR factorization.
415 [q, r] = qrupdate (q, r, u, v); 415 [q, r] = qrupdate (q, r, u, v);
416 else 416 else
442 fx = fun (x); 442 fx = fun (x);
443 jx = []; 443 jx = [];
444 endif 444 endif
445 445
446 if (! complexeqn && ! (isreal (fx) && isreal (jx))) 446 if (! complexeqn && ! (isreal (fx) && isreal (jx)))
447 error ("fsolve:notreal", "fsolve: non-real value encountered"); 447 error ("fsolve:notreal", "fsolve: non-real value encountered");
448 elseif (complexeqn && ! (isnumeric (fx) && isnumeric(jx))) 448 elseif (complexeqn && ! (isnumeric (fx) && isnumeric(jx)))
449 error ("fsolve:notnum", "fsolve: non-numeric value encountered"); 449 error ("fsolve:notnum", "fsolve: non-numeric value encountered");
450 elseif (any (isnan (fx(:)))) 450 elseif (any (isnan (fx(:))))
451 error ("fsolve:isnan", "fsolve: NaN value encountered"); 451 error ("fsolve:isnan", "fsolve: NaN value encountered");
452 endif 452 endif
453 endfunction 453 endfunction
454 454
455 function [fx, jx] = make_fcn_jac (x, fcn, fjac) 455 function [fx, jx] = make_fcn_jac (x, fcn, fjac)
456 fx = fcn (x); 456 fx = fcn (x);
457 if (nargout == 2) 457 if (nargout == 2)
458 jx = fjac (x); 458 jx = fjac (x);
459 endif 459 endif
460 endfunction 460 endfunction
461 461
462 %!function retval = f (p) 462 %!function retval = f (p)
463 %! x = p(1); 463 %! x = p(1);
464 %! y = p(2); 464 %! y = p(2);
465 %! z = p(3); 465 %! z = p(3);
466 %! retval = zeros (3, 1); 466 %! retval = zeros (3, 1);
467 %! retval(1) = sin(x) + y**2 + log(z) - 7; 467 %! retval(1) = sin(x) + y**2 + log(z) - 7;
493 %! [x, fval, info] = fsolve (@f, [-1, 1, 2, -1]); 493 %! [x, fval, info] = fsolve (@f, [-1, 1, 2, -1]);
494 %! assert (info > 0); 494 %! assert (info > 0);
495 %! assert (norm (x - x_opt, Inf) < tol); 495 %! assert (norm (x - x_opt, Inf) < tol);
496 %! assert (norm (fval) < tol); 496 %! assert (norm (fval) < tol);
497 497
498 %!function retval = f (p) 498 %!function retval = f (p)
499 %! x = p(1); 499 %! x = p(1);
500 %! y = p(2); 500 %! y = p(2);
501 %! z = p(3); 501 %! z = p(3);
502 %! retval = zeros (3, 1); 502 %! retval = zeros (3, 1);
503 %! retval(1) = sin(x) + y**2 + log(z) - 7; 503 %! retval(1) = sin(x) + y**2 + log(z) - 7;
512 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]); 512 %! [x, fval, info] = fsolve (@f, [ 0.5; 2.0; 2.5 ]);
513 %! assert (info > 0); 513 %! assert (info > 0);
514 %! assert (norm (x - x_opt, Inf) < tol); 514 %! assert (norm (x - x_opt, Inf) < tol);
515 %! assert (norm (fval) < tol); 515 %! assert (norm (fval) < tol);
516 516
517 %!function retval = f (p) 517 %!function retval = f (p)
518 %! x = p(1); 518 %! x = p(1);
519 %! y = p(2); 519 %! y = p(2);
520 %! z = p(3); 520 %! z = p(3);
521 %! retval = zeros (3, 1); 521 %! retval = zeros (3, 1);
522 %! retval(1) = sin(x) + y**2 + log(z) - 7; 522 %! retval(1) = sin(x) + y**2 + log(z) - 7;
539 %! x = 0:.5:5; 539 %! x = 0:.5:5;
540 %! noise = 1e-5 * sin (100*x); 540 %! noise = 1e-5 * sin (100*x);
541 %! y = exp (-a0*x) + b0 + noise; 541 %! y = exp (-a0*x) + b0 + noise;
542 %! c_opt = [a0, b0]; 542 %! c_opt = [a0, b0];
543 %! tol = 1e-5; 543 %! tol = 1e-5;
544 %! 544 %!
545 %! [c, fval, info, output] = fsolve (@(c) (exp(-c(1)*x) + c(2) - y), [0, 0]); 545 %! [c, fval, info, output] = fsolve (@(c) (exp(-c(1)*x) + c(2) - y), [0, 0]);
546 %! assert (info > 0); 546 %! assert (info > 0);
547 %! assert (norm (c - c_opt, Inf) < tol); 547 %! assert (norm (c - c_opt, Inf) < tol);
548 %! assert (norm (fval) < norm (noise)); 548 %! assert (norm (fval) < norm (noise));
549 549
554 %! y(3) = x(1) * x(2) - x(3)^2 + (3+2i); 554 %! y(3) = x(1) * x(2) - x(3)^2 + (3+2i);
555 555
556 %!test 556 %!test
557 %! x_opt = [-1+i, 1-i, 2+i]; 557 %! x_opt = [-1+i, 1-i, 2+i];
558 %! x = [i, 1, 1+i]; 558 %! x = [i, 1, 1+i];
559 %! 559 %!
560 %! [x, f, info] = fsolve (@cfun, x, optimset ("ComplexEqn", "on")); 560 %! [x, f, info] = fsolve (@cfun, x, optimset ("ComplexEqn", "on"));
561 %! tol = 1e-5; 561 %! tol = 1e-5;
562 %! assert (norm (f) < tol); 562 %! assert (norm (f) < tol);
563 %! assert (norm (x - x_opt, Inf) < tol); 563 %! assert (norm (x - x_opt, Inf) < tol);
564 564