2. Solution of Nonlinear Equations#
In this chapter we will see the implementation in Python of a number of numerical algorithms to solve nonlinear equations. The code provided is meant to be of inspiration for the students to build their own version of the implementation.
2.1. Bisection method#
We begin with an implementation of the Bisection Method.
import math
def interval_bisection(f,a:float,b:float,epsilon:float):
'''
Implements the Interval Bisection method.
Takes the following arguments:
- f: the function for which we want to find the root.
- a: the lower bound of the interval.
- b: the upper bound of the interval.
- epsilon: the desired precision.
'''
# Computes the number of iterations.
n_iterations = ((math.log(b-a, 2)- math.log(epsilon, 2))/(math.log(2,2)))-1
n_iterations = math.ceil(n_iterations)
print("Performing ",n_iterations, " iterations. ")
# Prints the header for the log.
print(f"{'Iteration':<10} {'a_k':<20} {'bk':<20} {'mk':<20} {'f(m_k)':<20} {'|f(m_k)|< epsilon':<20}")
# Sets the current interval. Initially it is the whole interval.
ak = a
bk = b
# Sets the counter of the number of iterations.
k = 1
# Initializes a boolean variable which states whether the method should terminate.
terminate = False
while not terminate:
# Finds the approximate root.
mk = (ak + bk)/2
# Computes the value of f at the approximate root.
fmk = f(mk)
# Prints a log.
print(f"{k:<10} {ak:<20.10f} {bk:<20.10f} {mk:<20.10f} {fmk:<20.10f} {abs(fmk)<epsilon:<20} ")
# Updates the interval.
# By checking whether the signs are the same we effectively avoid multiplying f(a_k) and f(m_k) which can cause numerical errors.
if math.copysign(1,f(ak)) == math.copysign(1, fmk):
ak = mk
else:
bk = mk
# Updates the iteration number.
k = k+1
# Checks whether we can terminate.
if k > n_iterations:
terminate = True
2.1.1. Example#
Let us test the method to find the zero of the function \(f(x)=\frac{1}{x}-\frac{1}{2}\). We look in the interval \([a,b]\) with \(a=\frac{3}{2}\) and \(b=3\).
def f(x):
return (1/x) - 0.5
a = 3/2
b = 3
precision = pow(10,-8)
interval_bisection(f,1.5,3,precision)
Performing 27 iterations.
Iteration a_k bk mk f(m_k) |f(m_k)|< epsilon
1 1.5000000000 3.0000000000 2.2500000000 -0.0555555556 0
2 1.5000000000 2.2500000000 1.8750000000 0.0333333333 0
3 1.8750000000 2.2500000000 2.0625000000 -0.0151515152 0
4 1.8750000000 2.0625000000 1.9687500000 0.0079365079 0
5 1.9687500000 2.0625000000 2.0156250000 -0.0038759690 0
6 1.9687500000 2.0156250000 1.9921875000 0.0019607843 0
7 1.9921875000 2.0156250000 2.0039062500 -0.0009746589 0
8 1.9921875000 2.0039062500 1.9980468750 0.0004887586 0
9 1.9980468750 2.0039062500 2.0009765625 -0.0002440215 0
10 1.9980468750 2.0009765625 1.9995117188 0.0001221001 0
11 1.9995117188 2.0009765625 2.0002441406 -0.0000610277 0
12 1.9995117188 2.0002441406 1.9998779297 0.0000305194 0
13 1.9998779297 2.0002441406 2.0000610352 -0.0000152583 0
14 1.9998779297 2.0000610352 1.9999694824 0.0000076295 0
15 1.9999694824 2.0000610352 2.0000152588 -0.0000038147 0
16 1.9999694824 2.0000152588 1.9999923706 0.0000019074 0
17 1.9999923706 2.0000152588 2.0000038147 -0.0000009537 0
18 1.9999923706 2.0000038147 1.9999980927 0.0000004768 0
19 1.9999980927 2.0000038147 2.0000009537 -0.0000002384 0
20 1.9999980927 2.0000009537 1.9999995232 0.0000001192 0
21 1.9999995232 2.0000009537 2.0000002384 -0.0000000596 0
22 1.9999995232 2.0000002384 1.9999998808 0.0000000298 0
23 1.9999998808 2.0000002384 2.0000000596 -0.0000000149 0
24 1.9999998808 2.0000000596 1.9999999702 0.0000000075 1
25 1.9999999702 2.0000000596 2.0000000149 -0.0000000037 1
26 1.9999999702 2.0000000149 1.9999999925 0.0000000019 1
27 1.9999999925 2.0000000149 2.0000000037 -0.0000000009 1
2.2. Fixed point method#
Here is a basic implementation of the fixed point (or functional iteration) method. The method takes as input the function \(f(x)\) for which we want to find the root, the iteration function \(g(x)\), the starting point \(x_0\), the desired precision (default \(\epsilon = 10^{-5}\)) and the maximum number of iterations it can run for (default \(100\)). It implements \(|f(x)|<\epsilon\) as the termination criteria.
def functional_iteration(f,g,x_0:float,epsilon:float=pow(10,-5),K:int=100):
'''
Implements the functional iteration method.
Takes the following arguments:
- f: the function for which we want to find the root.
- g: the iteration function.
- x_0: the initial guess.
- K: the maximum number of iterations.
'''
k = 1
x_k = x_0
function_value = f(x_k)
terminate = False
while not terminate:
x_kplus1 = g(x_k)
function_value = f(x_kplus1)
print(f"{k:d} {x_k:.10f} {x_kplus1:.10f} {function_value:.10f}")
x_k = x_kplus1
k = k + 1
if abs(function_value) < epsilon:
print("Function value smaller than epsilon")
terminate = True
if k > K:
print("Maximum number of iterations reached")
terminate = True
2.2.1. Example#
Let us use the method to find the root of \(f(x)=x-\sqrt{2x+3}\). We use the iteration function \(g(x)=\sqrt{2x+3}\), see Section 2.2 in the lecture notes. We arbitrarily start from \(x_0=4\) and require a precision of \(10^{-8}\).
import math
def f(x):
return x - math.sqrt(2*x+3)
def g(x):
return math.sqrt(2*x+3)
functional_iteration(f,g,4,pow(10,-8))
1 4.0000000000 3.3166247904 0.2128771233
2 3.3166247904 3.1037476670 0.0693621717
3 3.1037476670 3.0343854953 0.0229454759
4 3.0343854953 3.0114400194 0.0076291001
5 3.0114400194 3.0038109193 0.0025408817
6 3.0038109193 3.0012700376 0.0008467216
7 3.0012700376 3.0004233160 0.0002822140
8 3.0004233160 3.0001411020 0.0000940684
9 3.0001411020 3.0000470336 0.0000313558
10 3.0000470336 3.0000156778 0.0000104519
11 3.0000156778 3.0000052259 0.0000034840
12 3.0000052259 3.0000017420 0.0000011613
13 3.0000017420 3.0000005807 0.0000003871
14 3.0000005807 3.0000001936 0.0000001290
15 3.0000001936 3.0000000645 0.0000000430
16 3.0000000645 3.0000000215 0.0000000143
17 3.0000000215 3.0000000072 0.0000000048
Function value smaller than epsilon
2.3. Newton’s method#
Here is a basic implementation of Newton’s method. The method takes as input the function \(f\) for which we wish to find a root, its derivative \(f'\), the starting point \(x_0\), the tolerance \(\epsilon\) (default \(10^{-5}\)) and the maximum number of iterations \(K\) we want to perform (default \(K=100\)). It stops when either \(|f(x_k)|<\epsilon\), \(|x_{k+1}-x_k|<\epsilon\), or \(k=K\).
def newtons_method(f,fprime,x_0:float,epsilon:float=pow(10,-5),K:int=100):
'''
Implements the Newton's method.
Takes the following arguments:
- f: the function for which we want to find the root.
- fprime: the first derivative of the function f.
- x_0: the initial guess.
- epsilon: the desired precision.
- K: the maximum number of iterations.
'''
terminate = False
k = 1
x_k = x_0
function_value = f(x_k)
while not terminate:
x_kplus1 = x_k - f(x_k)/fprime(x_k)
function_value = f(x_kplus1)
print(f"{k:d} {x_k:.10f} {x_kplus1:.10f} {function_value:.10f}")
if abs(x_kplus1-x_k) < epsilon:
print("Approximate solutions are closer than epsilon.")
terminate = True
x_k = x_kplus1
k = k + 1
if abs(function_value) < epsilon:
print("Function value smaller than epsilon")
terminate = True
if k > K:
print("Maximum number of iterations reached")
terminate = True
2.3.1. Example#
As an example, let us compute the root of the function \(f(x)=\frac{1}{x}-\frac{1}{2}\).
def f(x):
return 1/x - 1/2
def fprime(x):
return -1/pow(x,2)
newtons_method(f,fprime,1,pow(10,-8))
1 1.0000000000 1.5000000000 0.1666666667
2 1.5000000000 1.8750000000 0.0333333333
3 1.8750000000 1.9921875000 0.0019607843
4 1.9921875000 1.9999694824 0.0000076295
5 1.9999694824 1.9999999995 0.0000000001
Function value smaller than epsilon
2.4. Horner’s method#
We provide here a basic implementation of Horner’s method to evaluate a polynomial $\(p(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0\)\( at a point \)x_0\( using nested multiplication. The method takes as input a dictionary 'a' containing the coefficients of the polynomial and \)x_0\(, the value at which the polynomial must be evaluated. In the dictionary the keys are the integers \)0,\ldots,n\( and the values are \)a_0,\ldots,a_n\(. It returns \)b_0=p(x_0)\( and the coefficients of the polynomial \)q(x)\(o f degree \)n-1\( such that \)\(p(x)=(x-x_0)q(x)+b_0\)$
def horners_method(a:dict, x0):
'''
Evaluates a polynomial using Horner's method.
- a: A dictionary containing the coefficients of the polynomial:
{0:a_0,1:a_1, ...,n:a_n}
- x0: The point at which to evaluate the polynomial.
Returns:
- The value of the polynomial at x0.
- The coefficients b of the polynomial q.
'''
# Determines the value of n
n = len(a) - 1
# Creates the dictionary of coefficients b
b = {}
# Sets b_n = a_n
b[n] = a[n]
# for k = n-1,...,0
for k in range(n-1,-1,-1):
# Finds b_k
b[k] = a[k] + b[k+1]*x0
# Returns p(x_0)= b_0 and b
return b[0], b
2.4.1. Example#
Let us use the method to evaluate the polynomial $\(p(x)=2x^3-18x^2+52x-48\)\( at \)x_0=5\( and find the coefficients of \)q(x)\( such that \)\(p(x)=(x-x_0)q(x)+p(x_0)\)$
coeffs = {0:-48,1:52,2:-18,3:2}
x0 = 5
b0,b = horners_method(coeffs, x0)
print(f"The value of the polynomial at x={x0} is: {b0}\n")
q = ' + '.join([f"{b[k]}x^{k-1}" for k in range(len(b)-1,0,-1)])
print(f"The polynomial q(x) is: {q}")
The value of the polynomial at x=5 is: 12
The polynomial q(x) is: 2x^2 + -8x^1 + 12x^0
2.5. Horner-Newton’s method#
We provide here a basic implementation of Horner-Newton’s method. The method takes as input a dictionary ‘a’ containing the coefficients of the polynomial, the starting point \(x_0\), a tolerance \(\epsilon\) (default \(10^{-5}\)), and the maximum number of iterations allowed (default \(100\)). It implements three stopping criteria: \(|f(x_k)|<\epsilon\), \(|x_{k+1}-x_k|<\epsilon\), \(k=\)’max_n_iterations’. One can optionally remove the termination criteria which are not considered appropriate for the given function.
def newton_horners_method(a:dict, x0,epsilon:pow(10,-5),max_n_iterations:100):
'''
Computes a root of the polynomial defined by the coefficients
in a starting from x0.
- a: A dictionary containing the coefficients of the polynomial:
{0:a_0,1:a_1, ...,n:a_n}
- x0: The starting point
Returns:
- The approximate solution at termination.
'''
# This variable will be set to true when the termination criteria is satisfied.
terminate = False
# Initializes the iteration counter.
i = 1
# Initializes the solution.
xi = x0
while not terminate:
# Within each Newton's method iteration we use Horner's method
# simultaneously on p(x_i) and q(x_i)
# Determines the value of n
n = len(a) - 1
# Creates the dictionary of coefficients b to evaluate p(xi)
b = {}
# Creates the dictionary of coefficients c to evaluate q(xi)
c = {}
# Sets b_n = a_n
b[n] = a[n]
c[n] = b[n]
# for k = n-1,...,1
for k in range(n-1,0,-1):
# Finds b_k
b[k] = a[k] + b[k+1]*xi
# Finds c_k
c[k] = b[k] + c[k+1]*xi
# Computes b_0
b0 = a[0] + b[1]*xi
# Performs Newton's method iteration
xiplus1 = xi - b0/c[1]
# Prints some information
print(f"{i:d} {xi:.10f} {xiplus1:.10f} {b0:.10f} {abs(xiplus1-xi):.10f}")
# Checks the termination criteria
if abs(xiplus1-xi) < epsilon:
print("Approximate solutions closer than epsilon")
terminate = True
if abs(b0) < epsilon:
print("Function value smaller than epsilon")
terminate = True
if k > max_n_iterations:
print("Maximum number of iterations reached")
terminate = True
# Updates the iteration number
i = i + 1
# Updates the solution
xi = xiplus1
return xi
2.5.1. Example#
Let us use the method to find a root of the polynomial $\(p(x)=2x^3-18x^2+52x-48\)\( starting from \)x0=5$
coeffs = {0:-48,1:52,2:-18,3:2}
x0 = 5
newton_horners_method(coeffs,x0,pow(10,-8),100)
1 5.0000000000 4.4545454545 12.0000000000 0.5454545455
2 4.4545454545 4.1510467894 3.2456799399 0.3034986652
3 4.1510467894 4.0253259290 0.7479702580 0.1257208604
4 4.0253259290 4.0009084519 0.1051846202 0.0244174770
5 4.0009084519 4.0000012353 0.0036387610 0.0009072166
6 4.0000012353 4.0000000000 0.0000049412 0.0000012353
7 4.0000000000 4.0000000000 0.0000000000 0.0000000000
Approximate solutions closer than epsilon
Function value smaller than epsilon
3.9999999999999964