Graphics Programs Reference
In-Depth Information
has beenreduced. Thiscompletes the first iteration cycle. In the nextcycle an-
other inverse quadratic interpolationis attempted and the process is repeatedun-
til the convergence criterion
|
x
x 3 |
issatisfied, where
ε
is aprescribed error
tolerance.
The inverse quadratic interpolationiscarried out with Lagrange'sthree-point
interpolant describedinArt. 3.2. Interchanging the roles of x and f , wehave
( f
f 2 )( f
f 3 )
( f
f 1 )( f
f 3 )
( f
f 1 )( f
f 2 )
x ( f )
=
f 3 ) x 1 +
f 3 ) x 2 +
f 2 ) x 3
( f 1
f 2 )( f 1
( f 2
f 1 )( f 2
( f 3
f 1 )( f 3
Setting f
=
0 and simplifying, weobtain for the estimate of the root
f 2 f 3 x 1 ( f 2
f 3 )
+
f 3 f 1 x 2 ( f 3
f 1 )
+
f 1 f 2 x 3 ( f 1
f 2 )
x
=
x (0)
=−
( f 1
f 2 )( f 2
f 3 )( f 3
f 1 )
The change in the root is
f 3 x 3 ( f 1
f 2 )( f 2
f 3 +
f 1 )
+
f 2 x 1 ( f 2
f 3 )
+
f 1 x 2 ( f 3
f 1 )
x
=
x
x 3 =
(4.2)
( f 2
f 1 )( f 3
f 1 )( f 2
f 3 )
brent
The function brent listedbelowis a simplifiedversion of the algorithmproposedby
Brent. It omits someofBrent'ssafeguards against slow convergence;it also uses a less
sophisticated convergence criterion.
function root = brent(func,a,b,tol)
%Findsarootoff(x)=0bycombiningquadratic
% interpolation with bisection (Brent's method).
% USAGE: root = brent(func,a,b,tol)
% INPUT:
% func = handle of function that returns f(x).
% a,b = limits of the interval containing the root.
% tol = error tolerance (default is 1.0e6*eps).
% OUTPUT:
% root = zero of f(x) (root = NaN if failed to converge).
if nargin < 4; tol = 1.0e6*eps; end
% First step is bisection
x1 = a; f1 = feval(func,x1);
if f1 == 0; root = x1; return; end
x2 = b; f2 = feval(func,x2);
if f2 == 0; root = x2; return; end
Search WWH ::




Custom Search