Roots of non-linear equations

This is the first practical problem we will treat. It is required for the solution of higher order ordinary differential equations with constant coefficients. It also occurs in many other problems during their analytical solution.

It is best to visualize this problem geometrically. We are looking for a value x = a for which a function y = f(x) vanishes, i.e., a root, where y = f(a) = 0.

We observe that, as a rule, the sign of y changes at a root. However, exceptions do occur, namely when the root is complex or multiple

Let us for the moment disregard these abnormal situations and concentrate on methods for finding simple roots. The simplest method, called bisection, employs the fact that the sign of f(x) changes at x=a, i.e, $f$ (x₁)*f(x₂)<0, x₁ < a < x₂. We then say that the values of the function at x_i bracket the root. However, this implies that the straight line between the points (x_i ,f(x_i)) must cross the x-axis at a point in between which is nearer to the root. While it is simple enough to write down the equation of this straight line and actually find the point, where it cuts the x-axis, bisection chooses the simpler way of halving the distance between these points, which yields x₃ = (x₁ + x₂)/2 , and computing the corresponding value f(x₃).

Obviously, if f(x₃) = 0, the root has been found. If not, it combines one of the first two points with the third point, ensuring that the selected couple of points again bracket the root . Since at each step the distance between the points x_i is halved, we have immediately at step n the error estimate (x₁-x₂)^-n. This is the only method in which you can set beforehand the number of cycles n.

However, note that the distance between the initial, bracketing values is also important. Often, it may be worthwhile to play with these values before starting the algorithm which, in pseudo-code, can be given the form:

Find a root of f(x), accurate within a fixed tolerance t, when initial, bracketing values x₁, x₂ are known:

x₃ = (x₁ + x₂)/2,
if f(x₁)´ f(x₃)/2< 0, then x₂ = x₃,
else x₁ = x₃,
end if x₁ = x₃,
repeat until x₁ - x₃ < t of f(x₃) = 0..

^{Note that this algorithm may fail if the
tolerance t is too small or if f(x) is
discontinuous.}

^Example

^{Find the positive root of the function}

^{Always try to}^sketch^{the curves in a root finding problem. It allows you
to understand what may happen later on. In this case, we see that
there must be a root}

^{between x = 0 and x = 1,
since f(0) = 1, f(1) = e - 4 < 0,
whence we can start:}

^{The root is at x=0.70148. You see
that the value of f(x) oscillates, because the half points
lie above and below at different distances from the root. The
objective is to find an accurate value of x rather than a}^{prescribed accuracy}^-^{difference
from zero}^{- of f(x)
at the root.}

^Example

^{Find for the tolerance t = 0.001
one of the roots of the polynomial}

^{We obtain readily}

^{In this manner, we find the roots
approximately by a table, or by a graph. Since we have already
located the root at -1, we can show by factorization that the
roots are 3}^{^1/2}^{. Locate the
positive root by bisection:}

^{We have stopped because both x and f(x)
satisfy the tolerance. In general, this will not be the case and
we must decide which criterion is to be applied. Note that 3}^{^1/2}^{=1.7320508
( or recall the value we found from the}^{continued
fraction}^!).

^{It is now natural to ask whether we have
made the best use of our routine calculations:}

^{Evaluations of f(x).}^{We arrive at the conclusion that}^bisection^{is really very arbitrary. It would be more
productive to use a more logical argument for the calculation of
the next point. I mentioned already that we could select instead
of the half point the intersection of the straight line, joining
the points f(x}1^),f(x2^):

^{We see from this figure that}

^{whence the}^{secant algorithm}^or^Regula
Falsi^becomes

^{,

if f(a)f(c') < 0, b
= c',

else

a = c',

end if,

REPEAT until a - b <}^e^{or f(c') = 0.}

^{Note that this algorithm may
fail if}^e^{is too small or f(x) is
discontinuous.}

^EXAMPLES

^{Reconsider the positive root of}

^{Start with the same values 0 and 1:}

^{Recompute the root 3 of the equation x³ +
x² - 3x - 3 = 0 by bracketing and linear interpolation:}

^{You see that we got there much faster!}

^{There has been developed an even faster
method by Van Wijngaarden-Dekker-Brent. It argues that once you
have computed three values of the function, you might just as
well use a quadratic rather than a straight line to find the
point of intersection with the x-axis. It uses for this
purpose}^{inverse
interpolation}^{. What is that,
you will ask.}

^{The Van
Wijngaarden-Dekker-Brent method}

^The^{pseudo code}^:

^{STEP 1 Select x}0^,x1^{so that y}0^.y1^{< 0

STEP 2 Bisect to obtain a third point (x}3^{, y}3^{)

STEP 3 Use inverse interpolation by a quadratic to find a new
point.

STEP 4 Select from these four points three which bracket the root
and lie closest together

STEP 5 Go to STEP 3}

^The^{inverse interpolation formula}^is

^{which becomes for the point of
intersection where y = 0}

^{. It is probably Brent who gave the
algorithm its present form by setting}

^{In general, P/Q will be small and
the root remains bracketed. If this ratio is not small, there
occurs another bisection between the points, which are closests.}

^EXAMPLE

^{Use again}

^{NEWTON'S METHOD}

^{You have been taught this method in the
first year. It employs the function as well as its first
derivative:}

^{i.e., we use the tangent to the curve to
reach the intersection with the x-axis.}^{If the value of x is close
to the root, Newton's method is very good.}

^EXAMPLE

^{Consider again}

^{Start at}

^ITERATION

^{In this case, we rewrite the equation y
= f(x) in the form}

^{and convert it into an}^{iteration algorithm}^{by writing}

^{When will this method
yield a root?}^{Consider the
difference}

^{and continue this operation, so that
eventually we find}

^{Thus, the method will converge as long as}

^{i.e., this derivative never
equals 1.}

^{Let us try to do our old example by this
method. If we write}

^{in the region of the root. On the other
hand, as we are looking for x such that e}^{^x}^{equals 2, take the inverse of this equality and form the}^{iteration relationship}

^{for which}

^{We find:}

^{RATES OF CONVERGENCE}

^{We have already done one such study for
iteration, in order to find the conditions under which it will
converge. Now we want to compare different methods and perhaps
find an explanation of the fact that, for example,}^{linear lnterpolation}^{yielded the result more quickly than bisection. It
would be nice to establish this fact}^analytically^.

^{For this purpose, we study the difference
between two consecutive values of x. As our approximations
improve, this difference will decrease. What is their rate of
decrease?}

^BISECTION

^{The determining quantity is the difference
x}1^{- x}0^{which is halved at each step. Since therefore}

^{denoting by r the root itself and
by e}j^{the error at step n, we find}

^{i.e., the error reduces linearly.}

^{LINEAR INTERPOLATION or REGULA
FALSI}

^{We recall the formula}

^{Transforming it in terms of errors, we
find}

^Using^{Taylor expansions}^{and taking into account the fact that y(r)
= 0, we find}

^{Now let e}n^=e^{^k}n-1^{, where k is
to be determined. Then the last formula yields:}

^i.e.,

^{Such a method is said to be}^superlinear^{, and that is why we have arrived more quickly at
the result!}

^{NEWTON'S METHOD}

^{We have found the formula}

^i.e.,

^{Thus, Newton's method is}^quadratic^{. However, as has been mentioned earlier, it is not
fool proof. The figure shows a situation when one has to be
closer to the root}

^{This is a consequence of the points not
bracketing the root. Newton's method takes you away from the
root!}

^EXAMPLES

^{Find n roots of the equation}

^{Both functions are well known to you. The
exponential function is always positive and the tan function
repeats itself along the x-axis periodically going from
minus infinity to plus infinity. Hence the roots will lie on the
positive branches of the tan function, gradually approaching the
asymptotes of the tan function.}

^BISECTION

^{We find}

^{Note the slow movement towards the root
which is brought about by the initial choice of x}1^{. The final value is x=0.303
with y(0.303) = -0.04795.}

^{REGULA FALSI}

^{Show that we find after three steps:}

^NEWTON

^{Show that we reach at the second step the
same accuracy as with the}^{regula falsi}^{in 3 steps!.}

^ITERATION

^{We must reorganize our equation for
iteration, keeping in mind that the}^{derivative of the iteration
function must be less than unity}^{. We have here the options}

^{Obviously, the first derivative is larger
than 1 at x = 1/4, whence we will use the second version,
starting with x = 1/4:}

^{where the last value yields y(0.29649)
= -0.000187895.}

^BRENT

^{We recall}

^{Bisection yielded:}

^{OTHER
TYPES OF ROOTS}

^{We will use now only iteration, because we
know that tan x is periodic, so that later roots on the
positive side will occur ahead of}

^{We use as starting value an approximate
value of the previous root and add to it}

^{The iteration formula was}

^{the value of the next root is}

^{and of the third root}

^{The second example is a polynomial which
could, for example, arise in the solution of a third order,
ordinary differential equation with constant coefficients:}

^{I have taken a simple polynomial with
simple integer roots to make our work easier. Since the last term
is then the product of the roots, we know that they are 4, -3 and
10! First of all, we must find two values of x for which
the signs of y are opposite. This is often most easily
done by plotting the curve. We use for the evaluations, of
course,}^{Horner's
scheme}^:

^{when bisection yields 3 which is a
root.???}

^{We now obtain from Horner a}^bonus^{, because (1,-14,40) are the coefficients of the
quadratic obtained by division of the original polynomial by (x +
3), and this equation yields immediately:}

^{We have used here the argument:}^{Consider the general polynomial}

^{, If we want to evaluate the polynomial
for x=g , we form in the Horner scheme the coefficients b}0^=a0^{, b}j^{= a}j^{+ b}j-1^{, j = 1,2, .
. . ,n. This operation yields b}n^{= y}0^{. However, the coefficients b}j^{are the coefficients in the polynomial obtained by division of y(x)
by (x - g) so that}

^{Obviously, y(g) = $b$}n^{.
Multiplying out the last equation and collecting the coefficients
of powers of x, we obtain Horner's scheme. Moreover, we
can repeat the operation of Horner's scheme on the $b$}j^{,
j = 0,1, . . .,n-1, to obtain C(x), D(x), . . . ,
until there is only one number left. Then $b$}n^,
cn-1^{, d}n-2^{, . . .
are the coefficients in y(x) after setting x =z + g.
However, that is the same thing as using the Taylor series to
expand y(x) about the point g, so that the
coefficients $b$}n^{, c}n-1^,
dn-2^{, . . . equal the Taylor coefficients
at the point :}

^{Thus, when we want to apply Newton's
method, we can find y'(g) by repeating it for the
coefficients $b$}j^{. Let us start
from the point $x$}0^{= -4, y}0^{= -112.}

^{Horner's scheme yields}

^whence

^x1^{= -4 - (-112)/134 = -3.16}

^x2^{= -3.16 - (-15.076)/97.4768 = -3.00533.}

last next