Non-Linear Solvers
Contents:
1 Newton’s Method
We are going to start with my favorite solver, Newton’s Method. Out of the four methods we will discuss, this one has the shortest pseudocode.
Compute $x_{n+1} = x_n - \frac{f(x_n)}{df(x_n)}$ until $|x_{n+1} - x_n| \leq eps$
If we look at the pieces:
- $x_{n+1}$ - next guess for a solution
- $x_n$ - current guess for a solution
- $f$ - function for which a solution is desired
- $df$ or $f’$ - derivative of $f$
- $eps$ - tolerance
Notice the requirements?
- We must have both $f$ and its derivative $df$.
- $df(x_n)$ must be non-zero.
1.1 Refining the Pseudocode
Our starting pseudocode was:
Compute $x_{n+1} = x_n - \frac{f(x_n)}{df(x_n)}$ until $|x_{n+1} - x_n| \leq eps$
While this code is complete, it is not something we can immediately start coding. Let us start by a slight change in notation.
while$|x_{n+1} - x_n| > eps$:$x_{n+1} = x_n - \frac{f(x_n)}{df(x_n)}$
This needs to be something we can write in code. Mathematical notation is usually difficult to preserve in formal code. Let us preemptively adjust our notation.
while abs(x_np1 - x_n) > eps:
x_np1 = x_n - (f(x_n) / df(x_n))
Let us clarify our inputs by making this a proper function.
def newtons_method(f, df, eps):
while abs(x_np1 - x_n) > eps:
x_np1 = x_n - (f(x_n) / df(x_n))
It is now clear that $f$, $df$, and $eps$ are arguments passed in when newtons_method is invoked. Both x_n and x_np1 are function local variables. This also brings up an issue. Where is our initial guess for $x_n$ ?
def newtons_method(f, df, x_n, eps):
while abs(x_np1 - x_n) > eps:
x_np1 = x_n - (f(x_n) / df(x_n))
That is much better.
1.2 Where is the Code?
You are probably definitely surprised that I have not started writing C++, Python, or Rust code yet. Why am I not ready, yet?
x_np1needs a better name.- What if $|x_n - x_{pn1}| > eps$ is never false?
- What if $df(x_n)$ is small?
- $x_n$ needs to be updated.
These can all be addressed quickly. Let us start by addressing 1 and 4:
def newtons_method(f, df, x_n, eps):
while abs(x_np1 - x_n) > eps:
next_x_n = x_n - (f(x_n) / df(x_n))
x_n = next_x_n
We are now left with two problems:
x_np1needs a better name.- What if $|x_n - x_{pn1}| > eps$ is never false?
- What if $df(x_n)$ is small?
- $x_n$
needs to be updated.
These two problems are both solved with a loop counter.
MAX_ITERATIONS = 100
def newtons_method(f, df, x_n, eps):
n = 0
while abs(x_np1 - x_n) > eps:
next_x_n = x_n - (f(x_n) / df(x_n))
x_n = next_x_n
n = n + 1
if n >= MAX_ITERATIONS:
break
return x_n
Yes, I did silently add the missing return statement. Did you notice it was missing?
1.3 Now We are Ready!
I am going to implement this code in Python 3 (3.7 to be exact). Let us make a quick draft.
newtons_method_1.py
#! /usr/bin/env python3
import sys
from typing import (Callable)
import typing
from fractions import (Fraction)
EPSILON = 10e-6
MAX_ITERATIONS = 100
def newtons_method(f, df, x_n, eps=EPSILON):
n = 0
next_x_n = x_n + 1000 * eps
while abs(x_n - next_x_n) > eps:
x_n = next_x_n
next_x_n = x_n - (f(x_n) / df(x_n))
if n >= MAX_ITERATIONS:
break
n += 1
return x_n
def main():
try:
initial_guess = float(sys.argv[1])
except IndexError as error:
print("Usage: {0} initial_guess".format(*sys.argv))
sys.exit(1)
except ValueError as error:
print("ERROR: {1} is not a valid number".format(*sys.argv))
print(" " + str(error))
sys.exit(2)
# Function (f) and its derivative (dx)
def f(x):
return (x ** 2) - 1
def df(x):
return 2 * x
try:
solution_newton = newtons_method(f, df, initial_guess)
fx_newton = f(solution_newton)
print("x = {:.4f} | f(x) = {:.4f}".format(solution_newton, fx_newton))
except ZeroDivisionError as error:
print(str(error))
if __name__ == "__main__":
main()
Notice how I used a direct translation of the pseudocode newtons_method function? Sometimes we can get away with a transliteration.
The tolerance, eps is now an argument to newtons_method (albeit one with a default argument value). This is something I usually define and use throughout a codebase.
Let us make some quick improvements.
newtons_method_2.py
#! /usr/bin/env python3
import sys
from typing import (Callable)
import typing
from fractions import (Fraction)
EPSILON = 10e-6
MAX_ITERATIONS = 100
def newtons_method(f, df, x_n, eps=EPSILON):
n = 0
next_x_n = x_n + 1000 * eps
while abs(x_n - next_x_n) > eps and n < MAX_ITERATIONS:
n += 1
x_n = next_x_n
next_x_n = x_n - (f(x_n) / df(x_n))
return n, x_n
def main():
try:
initial_guess = float(sys.argv[1])
except IndexError as error:
print("Usage: {0} initial_guess".format(*sys.argv))
sys.exit(1)
except ValueError as error:
print("ERROR: {1} is not a valid number".format(*sys.argv))
print(" " + str(error))
sys.exit(2)
# Function (f) and its derivative (dx)
def f(x):
return (x ** 2) - 1
def df(x):
return 2 * x
try:
num_iterations, solution_newton = newtons_method(f, df, initial_guess)
fx_newton = f(solution_newton)
output_str = "x = {:.4f} | f(x) = {:.4f} | {} iterations"
print(output_str.format(solution_newton, fx_newton, num_iterations))
except ZeroDivisionError as error:
print(str(error))
if __name__ == "__main__":
main()
I added the number of iterations as a return value. Python allows us to return multiple values as a tuple
return n, x_n
and unpack the result
num_iterations, solution_newton = newtons_method(f, df, initial_guess)
Now if I run the code with
./newtons_method_2.py 2
the output includes the number of iterations.
x = 1.0000 | f(x) = 0.0000 | 5 iterations
This is a reasonable implementation.
1.4 Cleaning Up The Implementation
It is now time to add type hints.
newtons_method_3.py
#! /usr/bin/env python3
import sys
from typing import (Callable, Tuple)
import typing
from fractions import (Fraction)
EPSILON = 10e-6
MAX_ITERATIONS = 100
def newtons_method(f: Callable[[float], float],
df: Callable[[float], float],
x_n: float,
eps: float = EPSILON) -> Tuple[int, float]:
n = 0
next_x_n = x_n + 1000 * eps
while abs(x_n - next_x_n) > eps and n < MAX_ITERATIONS:
n += 1
x_n = next_x_n
next_x_n = x_n - (f(x_n) / df(x_n))
return n, x_n
def main():
try:
initial_guess = float(sys.argv[1])
except IndexError as error:
print("Usage: {0} initial_guess".format(*sys.argv))
sys.exit(1)
except ValueError as error:
print("ERROR: {1} is not a valid number".format(*sys.argv))
print(" " + str(error))
sys.exit(2)
# Function (f) and its derivative (dx)
def f(x):
return (x ** 2) - 1
def df(x):
return 2 * x
try:
num_iterations, solution_newton = newtons_method(f, df, initial_guess)
fx_newton = f(solution_newton)
output_str = "x = {:.4f} | f(x) = {:.4f} | {} iterations"
print(output_str.format(solution_newton, fx_newton, num_iterations))
except ZeroDivisionError as error:
print(str(error))
if __name__ == "__main__":
main()
Let us add some docstrings in the Google/Sphinx style.
newtons_method_4.py
#! /usr/bin/env python3
import sys
from typing import (Callable, Tuple)
import typing
from fractions import (Fraction)
EPSILON = 10e-6
MAX_ITERATIONS = 100
def newtons_method(f: Callable[[float], float],
df: Callable[[float], float],
x_n: float,
eps: float = EPSILON) -> Tuple[int, float]:
"""
Compute an approximate solution to f(x) = 0 using Newton's Method.
Args:
f: mathematical function
df: derivative of f
x_n: an initial guess (i.e., x_0)
eps: tolerance within which two guesses will be considered equal.
Returns:
A Tuple with the number of iterations and approximate solution.
"""
n = 0
next_x_n = x_n + 1000 * eps
while abs(x_n - next_x_n) > eps and n < MAX_ITERATIONS:
n += 1
x_n = next_x_n
next_x_n = x_n - (f(x_n) / df(x_n))
return n, x_n
def main():
try:
initial_guess = float(sys.argv[1])
except IndexError as error:
print("Usage: {0} initial_guess".format(*sys.argv))
sys.exit(1)
except ValueError as error:
print("ERROR: {1} is not a valid number".format(*sys.argv))
print(" " + str(error))
sys.exit(2)
# Function (f) and its derivative (dx)
def f(x):
return (x ** 2) - 1
def df(x):
return 2 * x
try:
num_iterations, solution_newton = newtons_method(f, df, initial_guess)
fx_newton = f(solution_newton)
output_str = "x = {:.4f} | f(x) = {:.4f} | {} iterations"
print(output_str.format(solution_newton, fx_newton, num_iterations))
except ZeroDivisionError as error:
print(str(error))
if __name__ == "__main__":
main()
If we run pylint using python -m pylint newtons_method_4.py, we receive a few warnings:
************* Module newtons_method_4
newtons_method_4.py:1:0: C0111: Missing module docstring (missing-docstring)
newtons_method_4.py:13:0: C0103: Argument name "f" doesn't conform to snake_case naming style (invalid-name)
newtons_method_4.py:13:0: C0103: Argument name "df" doesn't conform to snake_case naming style (invalid-name)
newtons_method_4.py:31:4: C0103: Variable name "n" doesn't conform to snake_case naming style (invalid-name)
newtons_method_4.py:35:8: C0103: Variable name "n" doesn't conform to snake_case naming style (invalid-name)
newtons_method_4.py:43:0: C0111: Missing function docstring (missing-docstring)
newtons_method_4.py:57:4: C0103: Function name "f" doesn't conform to snake_case naming style (invalid-name)
newtons_method_4.py:57:4: C0103: Argument name "x" doesn't conform to snake_case naming style (invalid-name)
newtons_method_4.py:60:4: C0103: Function name "df" doesn't conform to snake_case naming style (invalid-name)
newtons_method_4.py:60:4: C0103: Argument name "x" doesn't conform to snake_case naming style (invalid-name)
newtons_method_4.py:6:0: W0611: Unused import typing (unused-import)
newtons_method_4.py:7:0: W0611: Unused Fraction imported from fractions (unused-import)
------------------------------------------------------------------
Your code has been rated at 6.76/10 (previous run: 6.76/10, +0.00)
In this case we can ignore the warnings dealing with x, n, f, and df. Why? We are sticking with the mathematical notation from the original pseudocode. That leaves us with:
newtons_method_4.py:1:0: C0111: Missing module docstring (missing-docstring)
newtons_method_4.py:43:0: C0111: Missing function docstring (missing-docstring)
newtons_method_4.py:6:0: W0611: Unused import typing (unused-import)
newtons_method_4.py:7:0: W0611: Unused Fraction imported from fractions (unused-import)
Since this is a quick example, we can forgo a module docstring. That leaves us with a docstring for main and a few unused imports.
newtons_method_5.py
#! /usr/bin/env python3
import sys
from typing import (Callable, Tuple)
from fractions import (Fraction)
EPSILON = 10e-6
MAX_ITERATIONS = 100
def newtons_method(f: Callable[[float], float],
df: Callable[[float], float],
x_n: float,
eps: float = EPSILON) -> Tuple[int, float]:
"""
Compute an approximate solution to f(x) = 0 using Newton's Method.
Args:
f: mathematical function
df: derivative of f
x_n: an intitial guess (i.e., x_0)
eps: tolerance within which two guesses will be considered equal.
Returns:
A Tuple with the number of iterations and approximate solution.
"""
n = 0
next_x_n = x_n + 1000 * eps
while abs(x_n - next_x_n) > eps and n < MAX_ITERATIONS:
n += 1
x_n = next_x_n
next_x_n = x_n - (f(x_n) / df(x_n))
return n, x_n
def main():
"""
This main function serves as the driver for the demo. Such functions
are not required in Python. However, we want to prevent unnecessary module
level (i.e., global) variables.
"""
try:
initial_guess = float(sys.argv[1])
except IndexError as error:
print("Usage: {0} initial_guess".format(*sys.argv))
sys.exit(1)
except ValueError as error:
print("ERROR: {1} is not a valid number".format(*sys.argv))
print(" " + str(error))
sys.exit(2)
# Function (f) and its derivative (dx)
def f(x):
return (x ** 2) - 1
def df(x):
return 2 * x
try:
num_iterations, solution_newton = newtons_method(f, df, initial_guess)
fx_newton = f(solution_newton)
output_str = "x = {:.4f} | f(x) = {:.4f} | {} iterations"
print(output_str.format(solution_newton, fx_newton, num_iterations))
except ZeroDivisionError as error:
print(str(error))
if __name__ == "__main__":
main()
We are left with one more warning:
newtons_method_5.py:6:0: W0611: Unused Fraction imported from fractions (unused-import)
So far we have been using the float data type (i.e., finite precision). Let us switch to the Python fractions library.
newtons_method_6.py
#! /usr/bin/env python3
import sys
from typing import (Callable, Tuple)
from fractions import (Fraction)
EPSILON = Fraction(1, 10**6)
MAX_ITERATIONS = 100
def newtons_method(f: Callable[[Fraction], Fraction],
df: Callable[[Fraction], Fraction],
x_n: Fraction,
eps: Fraction = EPSILON) -> Tuple[int, Fraction]:
"""
Compute an approximate solution to f(x) = 0 using Newton's Method.
Args:
f: mathematical function
df: derivative of f
x_n: an intitial guess (i.e., x_0)
eps: tolerance within which two guesses will be considered equal.
Returns:
A Tuple with the number of iterations and approximate solution.
"""
n = 0
next_x_n = x_n + 1000 * eps
while abs(x_n - next_x_n) > eps and n < MAX_ITERATIONS:
n += 1
x_n = next_x_n
next_x_n = x_n - (f(x_n) / df(x_n))
return n, x_n
def main():
"""
This main function serves as the driver for the demo. Such functions
are not required in Python. However, we want to prevent unnecessary module
level (i.e., global) variables.
"""
try:
initial_guess = float(sys.argv[1])
except IndexError as error:
print("Usage: {0} initial_guess".format(*sys.argv))
sys.exit(1)
except ValueError as error:
print("ERROR: {1} is not a valid number".format(*sys.argv))
print(" " + str(error))
sys.exit(2)
# Function (f) and its derivative (dx)
def f(x):
return (x ** 2) - 1
def df(x):
return 2 * x
try:
num_iterations, solution_newton = newtons_method(f, df, initial_guess)
fx_newton = f(solution_newton)
output_str = "x = {} | f(x) = {} | {} iterations"
print(output_str.format(solution_newton,
str(fx_newton),
str(num_iterations)))
except ZeroDivisionError as error:
print(str(error))
if __name__ == "__main__":
main()
That leaves us with
x = 1.0000000469590518 | f(x) = 9.391810573688986e-08 | 5 iterations
That is worse than before! Since are using fractions, rounding error propagates differently. This can be solved by tweaking
EPSILON = Fraction(1, 10**6)
However, 9.391810573688986e-08 is probably close enough to zero.
The generator expression and S.O.L.I.D discussion code can be found below.
newtons_method_7.py
#! /usr/bin/env python3
import sys
from typing import (Callable, Tuple)
from fractions import (Fraction)
EPSILON = Fraction(1, 10**6)
MAX_ITERATIONS = 100
def newtons_method(f: Callable[[Fraction], Fraction],
df: Callable[[Fraction], Fraction],
x_0: Fraction,
eps: Fraction = EPSILON) -> Tuple[int, Fraction]:
"""
Compute an approximate solution to f(x) = 0 using Newton's Method.
Args:
f: mathematical function
df: derivative of f
x_0: an initial guess (i.e., x_0)
eps: tolerance within which two guesses will be considered equal.
Returns:
A Tuple with the number of iterations and approximate solution.
"""
for n, numerical_solution in enumerate(__inner_newtons_method(f, df, x_0, eps)):
if n >= MAX_ITERATIONS:
break
return n, numerical_solution
def __inner_newtons_method(f: Callable[[Fraction], Fraction],
df: Callable[[Fraction], Fraction],
x_n: Fraction,
eps: Fraction = EPSILON) -> Tuple[int, Fraction]:
next_x_n = x_n + 1000 * eps
while abs(x_n - next_x_n) > eps:
x_n = next_x_n
next_x_n = x_n - (f(x_n) / df(x_n))
yield x_n
return x_n
def main():
"""
This main function serves as the driver for the demo. Such functions
are not required in Python. However, we want to prevent unnecessary module
level (i.e., global) variables.
"""
try:
initial_guess = float(sys.argv[1])
except IndexError as error:
print("Usage: {0} initial_guess".format(*sys.argv))
sys.exit(1)
except ValueError as error:
print("ERROR: {1} is not a valid number".format(*sys.argv))
print(" " + str(error))
sys.exit(2)
# Function (f) and its derivative (dx)
def f(x):
return (x ** 2) - 1
def df(x):
return 2 * x
try:
num_iterations, solution_newton = newtons_method(f, df, initial_guess)
fx_newton = f(solution_newton)
output_str = "x = {} | f(x) = {} | {} iterations"
print(output_str.format(solution_newton,
str(fx_newton),
str(num_iterations)))
except ZeroDivisionError as error:
print(str(error))
if __name__ == "__main__":
main()