Root finding - Numerical Analysis

#%config InlineBackend.figure_format = 'svg'
from pylab import *

Given a continuous function $f : [a,b] \to \re$ , find an $r \in [a,b]$ such that $f(r) = 0$ . Such an $r$ is called a root of $f$ or a zero of $f$ . There could be many roots for a given function.

At a root, the positive and negative parts of $f$ cancel one another. So we have to be careful about roundoff errors and care must be taken to evaluate the function so as to avoid loss of significance errors.

Since we work with finite precision on the computer, we cannot expect to find an $r$ that makes the function exactly zero. At best, we hope to keep relative error under control. In this sense, zeros far from the origin cannot be computed with great absolute accuracy. Since the purpose is to solve some problem that anyway involves approximations, it is not necessary to exactly locate the root, and an approximation is sufficient.

3.1Examples¶

Root finding problems arise in many problems as an intermediate step.

Example 2 (Solving ODE)

Consider the system of ODE

\dd{u(t)}{t} = f(u(t)), \qquad t \in [0,T]

(2)

where $u \in \re^N$ and $f : \re^N \to \re^N$ . We are given some initial condition $u(0) = u_0 \in \re^N$ . The backward Euler scheme is given by: $u^0 = u_0$ and

\frac{u^{n+1} - u^n}{\Delta t} = f(u^{n+1}), \qquad n=0,1,\ldots

(3)

This is a set of $N$ non-linear equations for the unknown solution $u^{n+1}$ . We can define

H : \re^N \to \re^N, \qquad H(u) = \frac{u - u^n}{\Delta t} - f(u)

(4)

The solution $u^{n+1}$ we seek is the root of the function $H$ , i.e., $H(u^{n+1}) = 0$ . We have to solve a multi-dimensional root problem for every time interval.

3.2Bracketing the root in 1-D¶

There are no general methods to find the roots of an arbitrary function. The user must have some knowledge of where the root might possibly lie. It is hence useful to first find some small interval in which we expect the root to lie.

Around a root, the function takes both positive and negative values. The IVT can be used to find an interval containing the root.

3Root finding

3.1Examples¶

3.2Bracketing the root in 1-D¶