Solution of Nonlinear Equations

Most equations that arise in engineering have no formula for their solution. There is no closed form for the root of x^{3} - x - 2 = 0 that you would want to use, and none at all for \cos x = x.

Numerical methods find an approximate answer to any wanted accuracy by repeating a simple calculation. This unit covers the errors that creep in and the classical root-finding methods.

All iteration tables in this book were produced by running the programs in bca/tools/check_numerical.py; nothing is typed by hand. Arithmetic is done in double precision and shown to six decimals.

1.1 Types of equations and locating roots

An algebraic equation is a polynomial equation (x^{3} - x - 2 = 0). A transcendental equation involves trigonometric, exponential or logarithmic functions (e^{x} = 3 x ,\ \cos x = x). A root (zero) of f is a number α with f (\alpha) = 0.

Intermediate value theorem. If f is continuous on [a ,\ b] and f (a) \cdot f (b) < 0, then f has at least one root in (a ,\ b). This sign change is how we bracket a root: tabulate f at convenient points and look for a sign change.

For f (x) = x^{3} - x - 2 we have f (1) = - 2 and f (2) = + 4, so a root lies in (1, 2). The exact root, correct to ten decimals, is \alpha = 1.5213797068.

1.2 Errors in computing

If the true value is \ast x \ast and the approximation is x̂:

• Absolute error E = |x − x̂|.
• Relative error Eᵣ = |x − x̂| / |x|; percentage error = Eᵣ × 100.
• A number is correct to \ast n \ast significant digits if its relative error is below 5 \times 10^{-n}.

Sources of error

1. Inherent (input) error: measurement or data error before the computation starts.
2. Round-off error: a computer stores numbers with finite digits. Adding 0.3333 three times gives 0.9999, not 1. Double precision carries about 16 significant digits; single precision about 7.
3. Truncation error: an infinite process (series, limit) is cut short. Approximating e^{x} at x = 0.5 by the first three Taylor terms gives 1.625000, and by four terms gives 1.645833, whereas the true value is 1.648721. The truncation errors are 0.023721 and 0.002888; each extra term cuts the error by roughly a factor of ten here.
4. Propagation of error: errors in the inputs grow or shrink through the operations. Subtracting nearly equal numbers (cancellation) is the worst case.

For a sum or difference the absolute errors add; for a product or quotient the relative errors add. For a function, the error is about |f' (x)| \cdot \Delta x.

1.3 The bisection method

Idea. Keep an interval [a ,\ b] with f (a) \cdot f (b) < 0. Take the midpoint c = \frac{a + b}{2}. Replace b by c if f (a) \cdot f (c) < 0, otherwise replace a by c. The interval halves each time.

Algorithm. (1) Choose a ,\ b with f (a) f (b) < 0 and a tolerance ε. (2) c = \frac{a + b}{2}. (3) If f (c) = 0 or \frac{b - a}{2} < \epsilon, stop: c is the root. (4) If f (a) f (c) < 0, set b = c; else set a = c. (5) Go to 2.

Error bound. After \ast n \ast steps the root lies in an interval of length \frac{b - a}{2^{n}}, so the number of steps needed for accuracy ε is n \geq \log_{2} \left( \frac{b - a}{\epsilon} \right).

For [1, 2] this is n \geq 10 for \epsilon = 10^{-3} and n \geq 20 for \epsilon = 10^{-6}.

Properties. Always converges if f is continuous and the bracket is valid; simple; linear convergence (one binary digit per step); slow; cannot find roots of even multiplicity (no sign change, such as x^{2} = 0).

1.4 The method of false position (regula falsi)

Instead of the midpoint, use the point where the chord joining (a ,\ f (a)) and (b ,\ f (b)) crosses the x-axis:

and keep the sub-interval whose ends still have opposite signs.

1.5 The Newton–Raphson method

Idea. Replace f by its tangent line at the current guess \ast x_{n} \ast and take the tangent's zero as the next guess. From the tangent y = f (x_{n}) + f' (x_{n}) (x - x_{n}) setting y = 0:

The same formula follows from the first-order Taylor expansion f (\alpha) \approx f (x_{n}) + (\alpha - x_{n}) f' (x_{n}) = 0.

Quadratic convergence. Near a simple root the error satisfies e_{n+1} \approx \left[ \frac{f'' (\alpha)}{2 f' (\alpha)} \right] \cdot e_{n}^{2}, i.e. the number of correct digits roughly doubles each step.

For our cubic the constant is \frac{f''}{2 f'} = 0.768.

The actual errors of the iterates x_{1} \ldots x_{4} were 2.14e-02, 3.59e-04, 9.92e-08, 7.55e-15 and the ratios \frac{e_{n+1}}{e_{n}^{2}} were 0.786, 0.768, 0.768, settling at the predicted constant.

Failure cases. (1) f' (x_{n}) = 0 makes the step undefined; a tiny f' throws the iterate far away. (2) A poor starting guess can diverge or wander to another root.

(3) Cycling: for f (x) = x^{3} - 2 x + 2 starting at 0 the iterates alternate between 0 and 1. (4) For f (x) = x^{\frac{1}{3}}, the step is x_{n+1} = x_{n} - 3 x_{n} = - 2 x_{n}, so it diverges from every starting point.

(5) At a multiple root convergence drops to linear. Always start near the root, found by a bracketing sweep or a few bisection steps.

1.6 The secant method

Newton needs the derivative. The secant method replaces f' (x_{n}) by the slope of the line through the last two points:

It needs two starting values but no bracket and only one new function evaluation per step.

The order of convergence is the golden ratio p \approx 1.618: between linear and quadratic. Per function evaluation it can beat Newton, because Newton needs two evaluations (f and f') per step.

1.7 Fixed-point iteration and convergence

Rewrite f (x) = 0 as x = g (x) and iterate x_{n+1} = g (x_{n}). A root is a fixed point of g.

Convergence theorem. If g and g' are continuous near the fixed point α and |g' (x)| < 1 on an interval around α containing the starting value, the iteration converges to α, and the smaller |g' (\alpha)|, the faster.

Convergence is linear with ratio |g' (\alpha)|. If |g' (\alpha)| > 1 it diverges.

1.8 Systems of nonlinear equations

For two unknowns, solve f_{1} (x ,\ y) = 0 ,\ f_{2} (x ,\ y) = 0 by Newton's method: at each step solve the linear system J (x_{n} ,\ y_{n}) \cdot (\Delta x ,\ \Delta y)^{T} = - (f_{1} ,\ f_{2})^{T}, where J is the Jacobian of partial derivatives, and update x_{n+1} = x_{n} + \Delta x ,\ y_{n+1} = y_{n} + \Delta y.

An alternative is fixed-point iteration on x = g_{1} (x ,\ y) ,\ y = g_{2} (x ,\ y), converging when the Jacobian of g has norm below 1.

A practical strategy: bracket the root, take a few bisection steps, then finish with Newton or secant. A good program stops when |x_{n+1} - x_{n}| < \epsilon or |f (x_{n})| < \epsilon, and also after a maximum number of iterations.

1.10 Seeing the ideas: Newton's tangent, bracketing and how fast each method converges

Most of the root-finding methods in this unit are geometric ideas before they are formulas. Plate 1.1 pairs the man behind the tangent method with a picture of his iteration, and the figures after it compare bracketing with fast iteration.

(a) Isaac Newton

(b) Tangent steps from x0 = 2

1.10.1 Eight reasons engineers need numerical methods

1. Most equations have no formula. A polynomial of degree five or higher has no general solving formula, and equations such as \cos x = x cannot be solved by algebra at all. Numerical methods still give the root to any accuracy we need. A calculator finds square roots in the same way.
2. Real data is a table, not a formula. Hydrologists in Nepal record the river flow at fixed hours, so they only have numbers. Interpolation, differentiation and integration work directly on such tables. No formula is needed.
3. Exact integrals are often impossible. The area under a curve such as e to the power minus x squared has no elementary answer. The trapezoidal and Simpson rules give it quickly and with a known error. Engineers use them to find the volume of water behind a d a m.
4. Large systems of equations need computers. A bridge model or an electrical network leads to hundreds of linear equations. Gauss elimination and iterative methods solve them in seconds. Solving them by hand would take days and would almost surely contain errors.
5. Differential equations rarely have closed forms. Cooling, population growth, vibration and circuit problems lead to equations that cannot be integrated exactly. Euler and Runge-Kutta methods step forward in time and give a usable answer. Weather forecasting is built on this idea.
6. Errors can be measured and controlled. Every method has an error bound or a stopping rule, such as the number of bisection steps n greater than or equal to log2((b - a) / tolerance). We know how accurate the answer is. This is better than a guess.
7. Speed and cost. Newton's method reaches six decimals in two or three steps, where bisection needs nineteen. Choosing a better method saves computer time in programs that run millions of times. This matters in games, simulations and mobile apps.
8. Computers use them anyway. Every scientific calculator, spreadsheet solver and programming library (NumPy, MATLAB) runs these algorithms inside. A student who knows them can choose the right tool and understand when the answer is unreliable. This is the base for data science and engineering software.

(a) A Patriot battery (Romanian Army): software clock matters

(b) 0.1 s stored in 23 bits: drift per hour

Formula and fact sheet

• Bisection: c = \frac{a + b}{2}; steps n \geq \log_{2} \left( \frac{b - a}{\epsilon} \right).
• False position: c = \frac{a f (b) - b f(a) }{f(b) - f (a)}.
• Newton: x_{n+1} = x_{n} - \frac{f (x_{n})}{f' (x_{n})} ;\ e_{n+1} \approx \left[ \frac{f''}{2 f'} \right] e_{n}^{2}.
• Secant: x_{n+1} = x_{n} - \frac{f (x_{n}) (x_{n} - x_{n-1})}{f (x_{n}) - f (x_{n-1})}; order 1.618.
• Fixed point: x_{n+1} = g (x_{n}), converges if |g'| < 1.
• Newton for systems: J \cdot \Delta = - F.
• Errors: absolute |x − x̂|; relative = absolute/|x|.

Chapter summary

• Numerical root-finding replaces an impossible formula by a convergent repetition; errors (round-off, truncation, cancellation) set the achievable accuracy.
• Bisection and false position are safe but linear; Newton is quadratic but needs f' and a good start; the secant method trades derivative for a 1.618 order; fixed-point iteration depends on |g'| < 1.
• Newton's method extends to systems through the Jacobian.

Practice

1. Define absolute, relative and percentage error. If \sqrt{2} is approximated by 1.414, find all three.
2. How many bisection steps guarantee an error below 10^{-3} for a root in [1, 2]? Below 10^{-6}?
3. Apply six steps of bisection to x^{2} - 10 = 0 on [3, 4]. What is the maximum possible error after six steps?
4. Apply Newton's method to x^{2} - 10 = 0 from x_{0} = 3. Why does it converge so fast?
5. Show that Newton's method applied to f (x) = x^{\frac{1}{3}} diverges.
6. Show that the rearrangement x = x^{3} - 2 of x^{3} - x - 2 = 0 diverges while x = (x + 2)^{\frac{1}{3}} converges.
7. Find the second iterate of the secant method for x^{3} - x - 2 = 0 from x_{0} = 1 ,\ x_{1} = 2.
8. Perform the first Newton step for x^{2} + y^{2} = 4 ,\ x y = 1 from (2, 0.5).
9. Explain why the false-position method may converge slowly.
10. Why is \sqrt{10001} - \sqrt{10000} computed badly in short arithmetic, and how can the formula be rearranged?

Answers and worked solutions

1. Absolute error = |x − x̂|; relative = absolute/|x|; percentage = relative × 100. For \sqrt{2} = 1.4142136 and 1.414: absolute 0.0002136, relative 0.0001510, percentage 0.0151 %.

2. n \geq \log_{2} \left( \frac{2 - 1}{\epsilon} \right): \epsilon = 10^{-3} needs n = 10 (2^{10} = 1024 > 1000) ;\ \epsilon = 10^{-6} needs n = 20 (2^{20} \approx 10^{6}).

3. The last midpoint is 3.171875.

After six steps the root lies in an interval of width 1 / 2^{6} = 0.015625, namely [3.156250, 3.171875], so the maximum possible error of the last midpoint is 0.015625 (the true \sqrt{10} = 3.162278, so the actual error is about 0.0096).

4. The iterates are 3.166667, 3.162281, 3.162278 (table in Section 1.5), so \sqrt{10} \approx 3.162278. Newton is quadratically convergent at a simple root: the correct digits double each step.

5. f (x) = x^{\frac{1}{3}} ,\ f' (x) = \left( \frac{1}{3} \right) x^{- \frac{2}{3}}, so \frac{f}{f'} = 3 x.

Then x_{n+1} = x_{n} - 3 x_{n} = - 2 x_{n}: the magnitude doubles each step and the sign alternates; the iterates run away from the root 0 for every x_{0} \neq 0.

6. For g (x) = x^{3} - 2 ,\ g' (\alpha) = 3 \alpha^{2} = 6.944 > 1; starting at 1.5 the iterates are 1.375, 0.5996, -1.784, -7.682.

For g (x) = (x + 2)^{\frac{1}{3}} ,\ g' (\alpha) = 0.1440 < 1, and the iterates converge as the table in Section 1.7 shows.

7. x_{2} = 2 - \frac{f (2) (2 - 1)}{f (2) - f (1)} = 2 - \frac{4}{4 + 2} = 1.333333. The next is x_{3} = 1.462687 (see Section 1.6 table).

8. J (2 ,\ 0.5) = \begin{bmatrix} 4 & 1 \\ 0.5 & 2 \end{bmatrix} ,\ f = (0.25 ,\ 0) ,\ \Delta = (- 0.066667 ,\ 0.016667), so (x_{1} ,\ y_{1}) = (1.933333 ,\ 0.516667).

9. For a convex or concave f on the bracket, the chord always crosses the axis on the same side of the root, so one end of the interval never moves and the other creeps up on the root linearly.

10. The two square roots are nearly equal, so their difference cancels the leading digits and leaves few significant ones. Multiply and divide by the conjugate: \sqrt{10001} - \sqrt{10000} = \frac{1}{\sqrt{10001} + \sqrt{10000}} = 0.0049999, which has no subtraction.

Model answers: medium (5 marks) and long (10 marks) questions

Q1 (5 marks). Find the root of x^{3} - 5 x + 1 = 0 lying in [0, 1] by the bisection method (seven iterations). State how many iterations guarantee an error below 10^{- 3}.

Check. f (0) = 1 > 0 and f (1) = - 3 < 0, so a root lies in [0, 1] (Figure 1.1 shows the same idea for another function).

n	a	b	c = \frac{a + b}{2}	f (c)
1	0	1	0.500000	-1.375000
2	0	0.5	0.250000	-0.234375
3	0	0.25	0.125000	+0.376953
4	0.125	0.25	0.187500	+0.069092
5	0.1875	0.25	0.218750	-0.083282
6	0.1875	0.21875	0.203125	-0.007244
7	0.1875	0.203125	0.195312	+0.030888

After seven steps the root lies in [0.195312, 0.203125]; the exact root is 0.201640.

Steps needed n ≥ log2((b - a)/epsilon) = \log 2 (1 / 0.001) = 9.97, so 10 iterations guarantee an error below 10^{- 3}.

Conclusion bisection is slow (linear) but always converges for a valid bracket.

Q2 (5 marks). Use Newton-Raphson to solve \cos x - x = 0 starting with x0 = 1 (three iterations).

Formula. x (n + 1) = x (n) - \frac{f (x (n))}{f' (x (n))} with f (x) = \cos x - x and f' (x) = - \sin x - 1.

n	x (n)	f (x (n))	f' (x (n))	x (n + 1)
0	1.000000	-0.459698	-1.841471	0.750364
1	0.750364	-0.018923	-1.681905	0.739113
2	0.739113	-0.000046	-1.673633	0.739085
3	0.739085	-0.000000	-1.673612	0.739085

The root is 0.739085 (the correct digits roughly double each step: 0, 1, 3, 6).

Why it is fast near a simple root the error satisfies e (n + 1) is about \left[ \frac{f''}{2 f'} \right] e (n)^{2} (quadratic convergence, Figure 1.2).

Cautions it fails when f' (x) = 0 and may diverge from a poor start.

Conclusion Newton is the preferred method when the derivative is available and the start is close.

Q3 (10 marks). Compare the bisection, false-position, secant, Newton-Raphson and fixed-point methods for solving f (x) = 0. Illustrate with x^{3} - x - 2 = 0 and say how you would choose a method.

Introduction. Root-finding methods fall into two families: bracketing methods (bisection, false position) that keep a root trapped inside an interval and always converge, and open methods (secant, Newton, fixed point) that start from one or two guesses, converge faster but can fail.

Plate 1.1 shows Newton's tangent idea; Figure 1.1 shows bisection; Figure 1.2 compares their speeds.

Methods. (1) Bisection: midpoint of a sign-changing interval; error halves each step; n ≥ log2((b - a)/epsilon).

(2) False position: the zero of the chord joining (a ,\ f (a)) and (b ,\ f (b)) replaces the midpoint; often faster but one end may stay fixed.

(3) Secant: like false position without keeping the bracket, using the last two points; order about 1.62. (4) Newton-Raphson: x - \frac{f}{f'}; order 2; needs f'.

(5) Fixed point: rewrite as x = g (x) and iterate; converges when |g' (x)| < 1 near the root, with linear rate |g'|.

Experiment. Solving x^{3} - x - 2 = 0 (root 1.521380) to within 10^{- 6} (error, or residual for false position) needs: bisection on [1, 2] 19 iterations; false position 13; fixed point x = (x + 2)^{\frac{1}{3}} from 1.5 6; secant from 1 and 2 6; Newton from 1.5 only 2.

Method	Order	Needs bracket	Needs f'	Iterations (this problem)	Main risk
Bisection	1 (linear, 0.5)	Yes	No	19	Slow
False position	1	Yes	No	13	One end sticks
Fixed point	1 (rate a b s g')	No	No	6	May diverge
Secant	about 1.62	No	No	6	May diverge
Newton-Raphson	2	No	Yes	2	f' = 0, bad start

Choice. Use bisection (or false position) to obtain a safe bracket, then switch to Newton or secant to polish the answer; this hybrid is what library routines such as Brent's method do. Use the secant method when f' is expensive or unavailable.

Conclusion there is no best method for every problem; safety and speed trade against each other, and good practice is to combine them.