Root finding algorithm
Contents
Root finding algorithm#
Where is the function crossing 0 ? How to find these points effectively ?
Introduction#
\(F(x)\) is a non linear function. We want to find points where this function reach zero if they exist. Root finding refers to the general problem of searching for a solution of an equation \(F(x) = 0\).
I - Naive method#
\(f\) is a function:
continuous on \( \mathopen{[}a\,,b\mathclose{]} \)
monotonous on \( \mathopen{[}a\,,b\mathclose{]} \)
such as \(f(a)f(b) < 0 \)
Method :
we choose a step \(dx\)
we test \(f(a + k.dx)\) until \(f\) changes its sign. \(k\) is a positive integer
## input : function f and interval [a,b] and our step is dx
import math
def naive():
x = a + dx
while (f(x) * f(a) > 0):
x = x + dx
return x
Be aware the result is an interval \( \mathopen{[}x - dx\,,x\mathclose{]} \). You can launch again the algorithm on this interval with a smaller step.
It is actually a really bad method.
II - Bisection#
Bisection is the simplest root finding algorithm. It can be applied to any type of continuous function \(f(x)\) on the interval \( \mathopen{[}a\,,b\mathclose{]} \) where the value of the function \(f(x)\) changes sign from \(a\) to \(b\). The idea is simple: divide the interval in two, a solution must exist within one subinterval, select the subinterval where the sign of changes and repeat.
\(f\) is a function:
continuous on \( \mathopen{[}a\,,b\mathclose{]} \)
monotonous on \( \mathopen{[}a\,,b\mathclose{]} \)
such as \(f(a)f(b) < 0 \)
Method :
We’re looking for a small interval containing the zero.
We split the start interval in half
We keep the interval that contains the 0
And we do it all over again
def bisection(f,a,b,N):
'''
Parameters
----------
f : function
The function for which we are trying to approximate a solution f(x)=0.
a,b : numbers
The interval in which to search for a solution. The function returns
None if f(a)*f(b) >= 0 since a solution is not guaranteed.
N : (positive) integer
The number of iterations to implement.
'''
if f(a)*f(b) >= 0:
print("Bisection method fails.")
return None
for i in range(1,N):
c = (a+b)/2
if(f(a) * f(c) < 0):
b = c
else:
a = c
return (a+b)/2
f = lambda x: x**2 - x - 1
bisection(f,1,2,25)
1.6180339753627777
Note : To obtain a 2 times higher precision, make another interation. Bisection method only used sign, not the values of the function. Bisection method is an improvement of the recursive naive method, where the step equals half the interval.
III - False position#
This method is the trial and error technique of using test (“false”) values for the variable and then adjusting the test value according to the outcome.
\(f\) is a function:
continuous on \( \mathopen{[}a\,,b\mathclose{]} \)
monotonous on \( \mathopen{[}a\,,b\mathclose{]} \)
such as \(f(a)f(b) < 0 \)
Method :
we do a linear approximation of \(f\) on \(\mathopen{[}a\,,b\mathclose{]}\)
we approximate the \(f\) curve by its chord
we solve the linear equation
we calculate the equation of the line passing through \((a, f(a)), (b, f(b))\)
we calculate its intersection with the abscissa axis
we get a new interval
we start again
You stop when \(f(c)\) is small enough
# The function is x^3 - x^2 + 2
def function( x ):
return (x * x * x - x * x + 2)
def falsePosition(a,b,N):
if function(a) * function(b) >= 0:
print("You have not assumed right a and b")
return -1
c = a # Initialize result
for i in range(N):
# Find the point that touches x axis
c = (a * function(b) - b * function(a))/ (function(b) - function(a))
# Check if the above found point is root
if function(c) == 0:
break
# Decide the side to repeat the steps
elif function(c) * function(a) < 0:
b = c
else:
a = c
print("The value of root is : " , '%.4f' %c, "with", '%i' %N ,"iterations")
falsePosition(-200, 300, 1000)
falsePosition(-200, 300, 100000)
The value of root is : -6.2309 with 1000 iterations
The value of root is : -1.0033 with 100000 iterations
This method can be unsuitable for some functions. The false position method is adapted to the functions whose behaviour is linear on \( \mathopen{[}a\,,b\mathclose{]} \), which allows to limit \( \mathopen{[}a\,,b\mathclose{]} \) more quickly than with the dichotomy method. It uses the values of \(f\) (unlike the bisection)
IV - Newton#
Newton’s method is a root finding method that uses linear approximation. In particular, we guess a solution \(x_0\) of the equation \(f(x)=0\), compute the linear approximation of \(f(x)\) at \(x_0\) and then find the \(x\)-intercept (\(x_1\)) of the linear approximation.
Hypothesis : \(x_0\) close to the result Iterative method
Method:
we calculate the tangent to the point \((x_i,f(x_i))\)
we calculate the intersection of the tangent with the abscissa axis
we get \(x_{i+1}\)
we start again
Tangent : $\(y = f'(x_0)x + y_0 - f'(x_0)x_0\)$
If \(y_1 = 0\) and \(y_0 = f(x_0)\), then
We stop when \(f(x)\) is small enough.
def newton(f,Df,x0,epsilon,max_iter):
'''Parameters
----------
f : function
Function for which we are searching for a solution f(x)=0.
Df : function
Derivative of f(x).
x0 : number
Initial guess for a solution f(x)=0.
epsilon : number
Stopping criteria is abs(f(x)) < epsilon.
max_iter : integer
Maximum number of iterations of Newton's method.
'''
xn = x0
for n in range(0,max_iter):
fxn = f(xn)
if abs(fxn) < epsilon:
print('Found solution after',n,'iterations.')
return xn
Dfxn = Df(xn)
if Dfxn == 0:
print('Zero derivative. No solution found.')
return None
xn = xn - fxn/Dfxn
print('Exceeded maximum iterations. No solution found.')
return None
f = lambda x: x**3 - x**2 - 1
Df = lambda x: 3*x**2 - 2*x
newton(f,Df,1,1e-10,10)
Found solution after 6 iterations.
1.4655712318767877
V - Secant#
We want to use Newton’s method but we don’t know the derivative \(f'\) of \(f\). We make an estimation of the derivative of \(f\) in \(x_0\) that we reinject in Newton’s method.
**Method : **
We’re evaluating the derivative of \(f(x_i)\).
We compute the tangent to the point \((x_i, f(x_i))\)
Xe calculate the intersection of the tangent with the abscissa axis
We get x_{i+1}
We do it all over again
We stop when \(f(x)\) is small enough.
Derivative : $\(f'(x_0) = lim_{x->x_0}\frac{f(x) - f(x_0)}{x - x_0} \)$
Numerical Estimation of the derivative (back derivative) : $\(f'(x) \simeq \frac{f(x) - f(x - \Delta x)}{\Delta x}\)$
Iteration : $\( x_{n+1} = x_n − \Delta x \frac{f(x_n)}{f(x_n) - f(x_n - \Delta x)} \)$
def function(x):
# we are taking equation x^3+x-1
f = pow(x, 3) + x - 1;
return f;
def secant(x1, x2, E):
n = 0; xm = 0; x0 = 0; c = 0;
if (function(x1) * function(x2) < 0):
while True:
# calculate the intermediate value
x0 = ((x1 * function(x2) - x2 * function(x1)) / (function(x2) - function(x1)));
# check if x0 is root of
# equation or not
c = function(x1) * function(x0);
# update the value of interval
x1 = x2;
x2 = x0;
# update number of iteration
n += 1;
# if x0 is the root of equation
# then break the loop
if (c == 0):
break;
xm = ((x1 * function(x2) - x2 * function(x1)) / (function(x2) - function(x1)));
if(abs(xm - x0) < E):
break;
print("Root of the given equation =",round(x0, 6));
print("Number of iterations = ", n);
else:
print("Can not find a root in ",
"the given inteval");
x1 = 0;
x2 = 1;
E = 0.0001;
secant(x1, x2, E);
Root of the given equation = 0.682326
Number of iterations = 5
Conclusion#
The root-finding problem is one of the most important computational problems. It arises in wide variety of practical applications.
📖 Sources#
Websites#
Name |
Author |
Date |
Description |
|---|---|---|---|
Patrick Walls |
2019 |
Explanation of mathematical computing |
Math classes#
Teacher |
Curriculum |
|---|---|
IMAC |