Module root_finding¶
Module generated from eponymous notebook.
Author: Brenton LeMesurier, brenton.lemesurier@unco.edu and lemesurierb@cofc.edu
Last Revised: February 25, 2021
A collection of functions for rootfinding, as used for example in MATH375 Assignments 2 and 3.
import numpy as np
def bisection(f, a, b, errorTolerance, maxIterations=100, demoMode=False):
"""Solve f(x) = 0 to within absolute error errorTolerance by the Bisection Method,
but avoiding infinite loops by giving up after maxIterations iterations.
Note: There is an optional input parameter "demoMode" which controls whether to
- print intermediate results (for "study" purposes), or to
- work silently (for "production" use).
"""
functionEvaluations = 0
# Whenever a (resp. c) is updated, also update fa (resp. fc).
fa = f(a)
functionEvaluations += 1
for iteration in range(1, maxIterations):
c = (a + b)/2
fc = f(c)
functionEvaluations += 1
if fa * fc < 0:
b = c
else:
a = c
fa = fc # Copy rather than redundantly evaluating f(a) for the new a, which is the old c
errorBound = (b-a)/2
if demoMode: print(f"{iteration=}, approximate root {c}, {errorBound=:0.3}")
if errorBound <= errorTolerance: break
root = c
return (root, errorBound, iteration, functionEvaluations)
def newton(f, Df, x0, errorTolerance, maxIterations=100, demoMode=False):
"""
Solve f(x) = 0 to within absolute error errorTolerance by Newton's Method
but avoiding infinite loops by giving up after maxIterations iterations.
Note: There is an optional input parameter "demoMode" which controls whether to
- print intermediate results (for "study" purposes), or to
- work silently (for "production" use).
The default is silence.
"""
functionEvaluations = 0
x = x0
for iteration in range(1, maxIterations + 1):
fx = f(x)
Dfx = Df(x)
functionEvaluations += 2
# Note: a more careful, robust code would check for the possibility of division by zero here.
dx = fx/Dfx
x -= dx # Aside: this is shorthand for "x = x - dx"
errorEstimate = abs(dx)
if demoMode: print(f"{iteration=}, {x=}, {errorEstimate=:0.3}")
if errorEstimate <= errorTolerance: break
return (x, errorEstimate, iteration, functionEvaluations)
def fpi(g, x_0, errorTolerance=1e-6, maxIterations=20, demoMode=False):
"""Fixed point iteration for approximately solving x = f(x),
x_0: the initial value
"""
x = x_0
for iteration in range(maxIterations):
x_new = g(x)
errorEstimate = np.abs(x_new - x)
x = x_new
if demoMode: print(f"x_{iteration} = {x}, {errorEstimate=:0.3}")
if errorEstimate <= errorTolerance: break
return (x, errorEstimate)
def false_position(f, a, b, errorTolerance=1e-15, maxIterations=15, demoMode=False):
"""Solve f(x)=0 in the interval [a, b] by the Method of False Position.
This code also illustrates a few ideas that I encourage, such as:
- Avoiding infinite loops, by using for loops sand break
- Avoiding repeated evaluation of the same quantity
- Use of descriptive variable names
- Use of "camelCase" to turn descriptive phrases into valid Python variable names
- An optional "demonstration mode" to display intermediate results.
"""
fa = f(a)
fb = f(b)
functionEvaluations = 2
for iteration in range(1, maxIterations+1):
if demoMode: print(f"\nIteration {iteration}:")
c = (a * fb - fa * b)/(fb - fa)
fc = f(c)
functionEvaluations += 1
if fa * fc < 0:
b = c
fb = fc # N.B. When b is updated, so must be fb = f(b)
else:
a = c
fa = fc
errorBound = b - a
if demoMode:
print(f"The root is in interval [{a}, {b}]")
print(f"The new approximation is {c}, with error bound {errorBound:0.3}, backward error {abs(fc):0.3}")
if errorBound < errorTolerance: break
# Whether we got here due to accuracy of running out of iterations,
# return the information we have, including an error bound:
root = c # the newest value is probably the most accurate
iterations = iteration
return (root, errorBound, iterations, functionEvaluations)
def secant(f, a, b, errorTolerance=1e-15, maxIterations=15, demoMode=False):
"""Solve f(x)=0 in the interval [a, b] by the Secant Method."""
# Some more descriptive names
x_older = a
x_more_recent = b
f_x_older = f(x_older)
f_x_more_recent = f(x_more_recent)
functionEvaluations = 2
for iteration in range(1, maxIterations+1):
x_new = (x_older * f_x_more_recent - f_x_older * x_more_recent)/(f_x_more_recent - f_x_older)
f_x_new = f(x_new)
functionEvaluations += 1
(x_older, x_more_recent) = (x_more_recent, x_new)
(f_x_older, f_x_more_recent) = (f_x_more_recent, f_x_new)
errorEstimate = abs(x_older - x_more_recent)
if demoMode: print(f"{iteration=}, x={x_more_recent}, {errorEstimate=:0.3}")
if errorEstimate < errorTolerance: break
# Whether we got here due to accuracy of running out of iterations,
# return the information we have, including an error estimate:
root = x_new
iterations = iteration
return (root, errorEstimate, iterations, functionEvaluations)
This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International