178
7 Solving Nonlinear Algebraic Equations
Implementation Given some Python implementation f(x) of our mathematical
function, a straightforward implementation of the above algorithm that quits after
finding one root, looks like
x = linspace(0, 4, 10001)
y = f(x)
root = None # Initialization
for i in range(len(x)-1):
if y[i]*y[i+1] < 0:
root = x[i] - (x[i+1] - x[i])/(y[i+1] - y[i])*y[i]
break # Jump out of loop
elif y[i] == 0:
root = x[i]
break # Jump out of loop
if root is None:
print(’Could not find any root in [{:g}, {:g}]’.format(x[0], x[-1]))
else:
print(’Find (the first) root as x={:.17f}’.format(root))
(See the file brute_force_root_finder_flat.py.)
Note the nice use of setting root to None: we can simply test if root is None
to see if we found a root and overwrote the None value, or if we did not find any
root among the tested points.
Running this program with some function, say f (x) = e −x 2 cos(4x) (which has
a solution at x =
π
8 ), gives the root 0.39269910538048097, which has an error of
2.4 × 10 −8 . Increasing the number of points with a factor of ten gives a root with an
error of 2.9 × 10 −10 .
After such a quick “flat” implementation of an algorithm, we should always try
to offer the algorithm as a Python function, applicable to as wide a problem domain
as possible. The function should take f and an associated interval [a, b] as input, as
well as a number of points (n), and return a list of all the roots in [a, b]. Here is our
candidate for a good implementation of the brute force root finding algorithm:
def brute_force_root_finder(f, a, b, n):
from numpy import linspace
x = linspace(a, b, n)
y = f(x)
roots = []
for i in range(n-1):
if y[i]*y[i+1] < 0:
root = x[i] - (x[i+1] - x[i])/(y[i+1] - y[i])*y[i]
roots.append(root)
elif y[i] == 0:
root = x[i]
roots.append(root)
return roots
(See the file brute_force_root_finder_function.py.)
This time we use another elegant technique to indicate if roots were found or not:
roots is an empty list if the root finding was unsuccessful, otherwise it contains all
the roots. Application of the function to the previous example can be coded as
7 Solving Nonlinear Algebraic Equations
Implementation Given some Python implementation f(x) of our mathematical
function, a straightforward implementation of the above algorithm that quits after
finding one root, looks like
x = linspace(0, 4, 10001)
y = f(x)
root = None # Initialization
for i in range(len(x)-1):
if y[i]*y[i+1] < 0:
root = x[i] - (x[i+1] - x[i])/(y[i+1] - y[i])*y[i]
break # Jump out of loop
elif y[i] == 0:
root = x[i]
break # Jump out of loop
if root is None:
print(’Could not find any root in [{:g}, {:g}]’.format(x[0], x[-1]))
else:
print(’Find (the first) root as x={:.17f}’.format(root))
(See the file brute_force_root_finder_flat.py.)
Note the nice use of setting root to None: we can simply test if root is None
to see if we found a root and overwrote the None value, or if we did not find any
root among the tested points.
Running this program with some function, say f (x) = e −x 2 cos(4x) (which has
a solution at x =
π
8 ), gives the root 0.39269910538048097, which has an error of
2.4 × 10 −8 . Increasing the number of points with a factor of ten gives a root with an
error of 2.9 × 10 −10 .
After such a quick “flat” implementation of an algorithm, we should always try
to offer the algorithm as a Python function, applicable to as wide a problem domain
as possible. The function should take f and an associated interval [a, b] as input, as
well as a number of points (n), and return a list of all the roots in [a, b]. Here is our
candidate for a good implementation of the brute force root finding algorithm:
def brute_force_root_finder(f, a, b, n):
from numpy import linspace
x = linspace(a, b, n)
y = f(x)
roots = []
for i in range(n-1):
if y[i]*y[i+1] < 0:
root = x[i] - (x[i+1] - x[i])/(y[i+1] - y[i])*y[i]
roots.append(root)
elif y[i] == 0:
root = x[i]
roots.append(root)
return roots
(See the file brute_force_root_finder_function.py.)
This time we use another elegant technique to indicate if roots were found or not:
roots is an empty list if the root finding was unsuccessful, otherwise it contains all
the roots. Application of the function to the previous example can be coded as