-
-
Notifications
You must be signed in to change notification settings - Fork 44k
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Consolidate Newton-Raphson implementations (#10859)
* updating DIRECTORY.md * updating DIRECTORY.md * Consolidate Newton-Raphson duplicates * Rename consolidated Newton-Raphson file * updating DIRECTORY.md * updating DIRECTORY.md * Fix doctest precision * Fix doctest precision again --------- Co-authored-by: github-actions <${GITHUB_ACTOR}@users.noreply.github.com>
- Loading branch information
1 parent
e5d6969
commit b98312c
Showing
5 changed files
with
105 additions
and
241 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file was deleted.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,45 +1,113 @@ | ||
# Implementing Newton Raphson method in Python | ||
# Author: Syed Haseeb Shah (github.com/QuantumNovice) | ||
# The Newton-Raphson method (also known as Newton's method) is a way to | ||
# quickly find a good approximation for the root of a real-valued function | ||
from __future__ import annotations | ||
|
||
from decimal import Decimal | ||
|
||
from sympy import diff, lambdify, symbols | ||
|
||
|
||
def newton_raphson(func: str, a: float | Decimal, precision: float = 1e-10) -> float: | ||
"""Finds root from the point 'a' onwards by Newton-Raphson method | ||
>>> newton_raphson("sin(x)", 2) | ||
3.1415926536808043 | ||
>>> newton_raphson("x**2 - 5*x + 2", 0.4) | ||
0.4384471871911695 | ||
>>> newton_raphson("x**2 - 5", 0.1) | ||
2.23606797749979 | ||
>>> newton_raphson("log(x) - 1", 2) | ||
2.718281828458938 | ||
""" | ||
The Newton-Raphson method (aka the Newton method) is a root-finding algorithm that | ||
approximates a root of a given real-valued function f(x). It is an iterative method | ||
given by the formula | ||
x_{n + 1} = x_n + f(x_n) / f'(x_n) | ||
with the precision of the approximation increasing as the number of iterations increase. | ||
Reference: https://en.wikipedia.org/wiki/Newton%27s_method | ||
""" | ||
from collections.abc import Callable | ||
|
||
RealFunc = Callable[[float], float] | ||
|
||
|
||
def calc_derivative(f: RealFunc, x: float, delta_x: float = 1e-3) -> float: | ||
""" | ||
Approximate the derivative of a function f(x) at a point x using the finite | ||
difference method | ||
>>> import math | ||
>>> tolerance = 1e-5 | ||
>>> derivative = calc_derivative(lambda x: x**2, 2) | ||
>>> math.isclose(derivative, 4, abs_tol=tolerance) | ||
True | ||
>>> derivative = calc_derivative(math.sin, 0) | ||
>>> math.isclose(derivative, 1, abs_tol=tolerance) | ||
True | ||
""" | ||
return (f(x + delta_x / 2) - f(x - delta_x / 2)) / delta_x | ||
|
||
|
||
def newton_raphson( | ||
f: RealFunc, | ||
x0: float = 0, | ||
max_iter: int = 100, | ||
step: float = 1e-6, | ||
max_error: float = 1e-6, | ||
log_steps: bool = False, | ||
) -> tuple[float, float, list[float]]: | ||
""" | ||
x = symbols("x") | ||
f = lambdify(x, func, "math") | ||
f_derivative = lambdify(x, diff(func), "math") | ||
x_curr = a | ||
while True: | ||
x_curr = Decimal(x_curr) - Decimal(f(x_curr)) / Decimal(f_derivative(x_curr)) | ||
if abs(f(x_curr)) < precision: | ||
return float(x_curr) | ||
Find a root of the given function f using the Newton-Raphson method. | ||
:param f: A real-valued single-variable function | ||
:param x0: Initial guess | ||
:param max_iter: Maximum number of iterations | ||
:param step: Step size of x, used to approximate f'(x) | ||
:param max_error: Maximum approximation error | ||
:param log_steps: bool denoting whether to log intermediate steps | ||
:return: A tuple containing the approximation, the error, and the intermediate | ||
steps. If log_steps is False, then an empty list is returned for the third | ||
element of the tuple. | ||
:raises ZeroDivisionError: The derivative approaches 0. | ||
:raises ArithmeticError: No solution exists, or the solution isn't found before the | ||
iteration limit is reached. | ||
>>> import math | ||
>>> tolerance = 1e-15 | ||
>>> root, *_ = newton_raphson(lambda x: x**2 - 5*x + 2, 0.4, max_error=tolerance) | ||
>>> math.isclose(root, (5 - math.sqrt(17)) / 2, abs_tol=tolerance) | ||
True | ||
>>> root, *_ = newton_raphson(lambda x: math.log(x) - 1, 2, max_error=tolerance) | ||
>>> math.isclose(root, math.e, abs_tol=tolerance) | ||
True | ||
>>> root, *_ = newton_raphson(math.sin, 1, max_error=tolerance) | ||
>>> math.isclose(root, 0, abs_tol=tolerance) | ||
True | ||
>>> newton_raphson(math.cos, 0) | ||
Traceback (most recent call last): | ||
... | ||
ZeroDivisionError: No converging solution found, zero derivative | ||
>>> newton_raphson(lambda x: x**2 + 1, 2) | ||
Traceback (most recent call last): | ||
... | ||
ArithmeticError: No converging solution found, iteration limit reached | ||
""" | ||
|
||
def f_derivative(x: float) -> float: | ||
return calc_derivative(f, x, step) | ||
|
||
a = x0 # Set initial guess | ||
steps = [] | ||
for _ in range(max_iter): | ||
if log_steps: # Log intermediate steps | ||
steps.append(a) | ||
|
||
error = abs(f(a)) | ||
if error < max_error: | ||
return a, error, steps | ||
|
||
if f_derivative(a) == 0: | ||
raise ZeroDivisionError("No converging solution found, zero derivative") | ||
a -= f(a) / f_derivative(a) # Calculate next estimate | ||
raise ArithmeticError("No converging solution found, iteration limit reached") | ||
|
||
|
||
if __name__ == "__main__": | ||
import doctest | ||
from math import exp, tanh | ||
|
||
doctest.testmod() | ||
|
||
# Find value of pi | ||
print(f"The root of sin(x) = 0 is {newton_raphson('sin(x)', 2)}") | ||
# Find root of polynomial | ||
print(f"The root of x**2 - 5*x + 2 = 0 is {newton_raphson('x**2 - 5*x + 2', 0.4)}") | ||
# Find value of e | ||
print(f"The root of log(x) - 1 = 0 is {newton_raphson('log(x) - 1', 2)}") | ||
# Find root of exponential function | ||
print(f"The root of exp(x) - 1 = 0 is {newton_raphson('exp(x) - 1', 0)}") | ||
def func(x: float) -> float: | ||
return tanh(x) ** 2 - exp(3 * x) | ||
|
||
solution, err, steps = newton_raphson( | ||
func, x0=10, max_iter=100, step=1e-6, log_steps=True | ||
) | ||
print(f"{solution=}, {err=}") | ||
print("\n".join(str(x) for x in steps)) |
This file was deleted.
Oops, something went wrong.
This file was deleted.
Oops, something went wrong.