Help-Root

NumericalMethods[Root] - approximates roots of functions

Calling Sequence
    Root(f, xinit=a, method=m, tolerance=t, options)
    Root(f, xinit=[a,b], method=m, tolerance=t, options)

Parameters
    f - algebraic; expression with variable x representing a continuous function
    a - numeric; initial approximate root or the first of two if required by the method
    b - numeric; second approximate root if required by the method
    m - name; method to be used
    t - positive; tolerance for the stopping condition
    options - (optional) parameters of the form keyword = value, where keyword is one of stoppingcriterion, output, maxiterations, maxnumber, digits.

Description
The Root command approximates a root of the given function (bisection, regula falsi, secant, Newton, and fixed point approach). For comparison, some other methods are included: Illinois, Brent, Chandrupatla (1997). Its main purpose is to exhibit behaviour of basic methods and be accessible to students not versed in Maple. Acordingly, the main focus is on ease of use and informative output, not performace. Controls are similar to those of Student[NumericalAnalysis][Roots], with some changes targeted at convenience of use and compatibility with other commands from this package.

There is no default method or default tolerance, they must be user specified.

This command is part of the NumericalMethods package, so it can be used in the form Root(..) only after executing the command with(NumericalMethods). However, it can always be accessed through the long form of the command by using NumericalMethods[Root](..).

Options

method = bisection, regfalsi, secant, newton, [newton,numeric], brent, new, fixedpoint, or [fixedpoint,algebraic]
The method used to approximate the root(s) of f numerically.
- bisection: the bisection metod. Requires two initial values xinit = [a,b], where a<b and f changes its sign at these points. The midpoint is taken as the outcome of an iteration.
- regfalsi: the regula falsi metod. Requires two initial values xinit = [a,b], where a<b and f changes its sign at these points.
- secant: the secant metod. Requires two initial values xinit = [a,b] taken as x₀, x₁.
- newton: the Newton method (tangents). Requires one initial value xinit = a.
- [newton,m] where m = numeric: the modified Newton method x_{k +1} = x_k - m f (x_k) / f '(x_k) for roots of multiplicity m. Requires one initial value xinit = a.
- illinois: the Illinois metod. Requires two initial values xinit = [a,b], where a<b and f changes its sign at these points. Iteration number is preceded by "M" if it was a modified step.
- brent: the Brent metod. Requires two initial values xinit = [a,b], where a<b and f changes its sign at these points. Iteration number is preceded by "B" (bisection), "S" (secant) or "Q" (inverse quadratic) to indicate the method used by the algorithm.
- new: the Chandrupatla metod (1997). Requires two initial values xinit = [a,b], where a<b and f changes its sign at these points. Iteration number is preceded by "B" (bisection) or "Q" (inverse quadratic) to indicate the method used by the algorithm.
- fixedpoint: converts f(x) = 0 into a fixed-point problem of the form g(x) = x solved by plain iteration. The default conversion uses iterator function g(x) = f(x)+x. Requires one initial value xinit = a.
- [fixedpoint,g] where g = algebraic: solves the fixed-point problem g(x) = x by plain iteration, with user-supplied iterator function g given as expression with variable x. The procedure does not check on consistency, that is, whether g(x) = x is equivalent to f(x) = 0. Requires one initial value xinit = a.

The options argument can be one or more of the following equations.

output = information, iterates, value, or sequence
Specifies what information is given by the procedure. In all cases, the iteration index k is printed at every iteration step, so that the user can see the progress and knows how many iterations it took to produce the final approximation.
- information returns the final numerical approximation of the root. During iterations, data are shown at each step: the generated approximation, function value at the current approximation and value test used in stoppingcriterion to determine whether the algorithm should stop. For bracketing methods also the brackets are shown.
- iterates returns the final numerical approximation of the root. During iterations, generated approximations are shown. For bracketing methods also the brackets are shown.
- value returns the final numerical approximation of the root.
- sequence: returns all approximations produced by iterations as a sequence.
The default value is output = information.

maxiterations = posint
The maximum number of iterations to perform. If it is reached, the algorithm is interrupted and a message is displayed. The return as defined by output is still produced.
The default value is maxiterations = 50.

maxnumber = positive
The maximal number that is not to be exceeded by approximations or function values during iterations. If it does happen, the algorithm is interrupted and a message is displayed. The return as defined by output is still produced.
The default value is maxnumber = 1000000.

digits = posint
Specifies how many decimal places should be shown when output = information or output = iterates is specified.
If digits is set to a value greater than Digits, the latter is set to digits for the duration of this procedure.
The default value is digits = Digits.

Examples

> f:=exp(1/2-x)-1;

This function has root 1/2.
> Root(f, xinit=[-2,1], method=bisection, tolerance=0.01);
k=00 a= -2.0000000000 b= 1.0000000000   x= -0.5000000000   f(x)= 1.7182818280 test= 1.5000000000
k=01 a= -0.5000000000 b= 1.0000000000   x= 0.2500000000   f(x)= 0.2840254170 test= 0.7500000000
k=02 a= 0.2500000000 b= 1.0000000000   x= 0.6250000000   f(x)= -0.1175030974 test= 0.3750000000
k=03 a= 0.2500000000 b= 0.6250000000   x= 0.4375000000   f(x)= 0.0644944590 test= 0.1875000000
k=04 a= 0.4375000000 b= 0.6250000000   x= 0.5312500000   f(x)= -0.0307667655 test= 0.0937500000
k=05 a= 0.4375000000 b= 0.5312500000   x= 0.4843750000   f(x)= 0.0157477090 test= 0.0468750000
k=06 a= 0.4843750000 b= 0.5312500000   x= 0.5078125000   f(x)= -0.0077820617 test= 0.0234375000
k=07 a= 0.4843750000 b= 0.5078125000   x= 0.4960937500   f(x)= 0.0039138890 test= 0.0117187500
k=08 a= 0.4960937500 b= 0.5078125000   x= 0.5019531250   f(x)= -0.0019512189 test= 0.0058593750

If the output turns out to be too wide, one can use the digits option, digits = 5 is often enough for educational purposes. In particular it does not make much sense to have digits significantly higher than tolerance. We will also try to cut the iteration short, leave out extra information and use the other bracketing method - regula falsi (false position).
> Root(f, xinit=[-2,1], method=regfalsi, tolerance=0.01, digits=3, maxiterations=5, output=iterates);
k=00 a= -2.000 b= 1.000   x= 0.898
k=01 a= -2.000 b= 0.898   x= 0.815
k=02 a= -2.000 b= 0.815   x= 0.749
k=03 a= -2.000 b= 0.749   x= 0.696
k=04 a= -2.000 b= 0.696   x= 0.654
k=05 a= -2.000 b= 0.654   x= 0.620
Maximal number of iterations reached.

The return is in the precision given by Digits. Note that the last approximation is returned even when the procedure is interrupted prematurely, in case you want to work further with it.
For the secant method we will try the minimal output.
> Root(f, xinit=[1,2], method=secant, tolerance=0.001, output=value);
k=00 k=01 k=02 k=03 k=04 k=05 k=06 k=07 k=08

We still get to see how many iterations it took. For our next example we use the Newton method, a different stopping criterion and ask for a sequence.
> Root(f, xinit=2, method=newton, tolerance=0.0001, stoppingcriterion=value, output=sequence);
k=00 k=01 k=02 k=03 k=04 k=05 k=06

We store them and supress the output using colon.
> sol:=Root(f, xinit=2, method=newton, tolerance=0.0001, stoppingcriterion=relative, output=sequence):
k=00 k=01 k=02 k=03 k=04 k=05 k=06 k=07

Note that when we set the stopping condition to relative, it took one more iteration.
> sol[7];

We can also ask for more digits to show than the default, the precision of calculations is reset accordingly for the duration.
> sol:=Root(f, xinit=2, method=newton, tolerance=1E-17, output=sequence, digits=25):
k=00 k=01 k=02 k=03 k=04 k=05 k=06 k=07 k=08 k=09 Root (probably) found

The procedure tells us that f (x) = 0 in the computer arithmetics, so we perhaps found the root.
Now we can do some analysis, for instance try to estimate the order of this method. Since Digits is now back to its default, many iterates would be rounded to the precise value of the root, so we need to increase precision also for the next calculation.
> Digits:=25:
last:=8:
for k from last-2 to last do err[k]:=abs(1/2-sol[k]) end do;
q:=ln(err[last-1]/err[last])/ln(err[last-2]/err[last-1]);
# estimate of the order of the method
Digits:=10:

Now for some fixed point iteration. We attempt to solve the equation x³ + x - 2 = 0 whose solution is 1. Given that fixed point iteration often fails, we set the limit on large numbers a bit lower. It turns out we were right to be careful.
> Root(x^3+x-2, xinit=0.5, method=fixedpoint, tolerance=0.01, maxnumber=1000);
k=00 x= 0.5000000000   f(x)= -1.3750000000 test= 1.3750000000
k=01 x= -0.8750000000   f(x)= -3.5449218750 test= 1.3750000000
k=02 x= -4.4199218750   f(x)=-92.7662311200 test= 3.5449218750
k=03 x=-97.1861529900   f(x)=-918036.8171000000 test= 92.7662311200
Numbers in iteration too large.

We try our own conversion, solving the given equation for the x in the cubic power.
> Root(x^3+x-2, xinit=0.5, method=[fixedpoint,(2-x)^(1/3)], tolerance=0.01);
k=00 x= 0.5000000000   f(x)= -1.3750000000 test= 1.3750000000
k=01 x= 1.1447142430   f(x)= 0.6447142450 test= 0.6447142430
k=02 x= 0.9492277221   f(x)= -0.1954865210 test= 0.1954865209
k=03 x= 1.0166454830   f(x)= 0.0674177600 test= 0.0674177609
k=04 x= 0.9944204320   f(x)= -0.0222250510 test= 0.0222250510
k=05 x= 1.0018564080   f(x)= 0.0074359770 test= 0.0074359760