|
1 | | -# Implementing Newton Raphson method in Python |
2 | | -# Author: Syed Haseeb Shah (github.com/QuantumNovice) |
3 | | -# The Newton-Raphson method (also known as Newton's method) is a way to |
4 | | -# quickly find a good approximation for the root of a real-valued function |
5 | | -from __future__ import annotations |
6 | | - |
7 | | -from decimal import Decimal |
8 | | - |
9 | | -from sympy import diff, lambdify, symbols |
10 | | - |
11 | | - |
12 | | -def newton_raphson(func: str, a: float | Decimal, precision: float = 1e-10) -> float: |
13 | | - """Finds root from the point 'a' onwards by Newton-Raphson method |
14 | | - >>> newton_raphson("sin(x)", 2) |
15 | | - 3.1415926536808043 |
16 | | - >>> newton_raphson("x**2 - 5*x + 2", 0.4) |
17 | | - 0.4384471871911695 |
18 | | - >>> newton_raphson("x**2 - 5", 0.1) |
19 | | - 2.23606797749979 |
20 | | - >>> newton_raphson("log(x) - 1", 2) |
21 | | - 2.718281828458938 |
| 1 | +""" |
| 2 | +The Newton-Raphson method (aka the Newton method) is a root-finding algorithm that |
| 3 | +approximates a root of a given real-valued function f(x). It is an iterative method |
| 4 | +given by the formula |
| 5 | + |
| 6 | +x_{n + 1} = x_n + f(x_n) / f'(x_n) |
| 7 | + |
| 8 | +with the precision of the approximation increasing as the number of iterations increase. |
| 9 | + |
| 10 | +Reference: https://en.wikipedia.org/wiki/Newton%27s_method |
| 11 | +""" |
| 12 | +from collections.abc import Callable |
| 13 | + |
| 14 | +RealFunc = Callable[[float], float] |
| 15 | + |
| 16 | + |
| 17 | +def calc_derivative(f: RealFunc, x: float, delta_x: float = 1e-3) -> float: |
| 18 | + """ |
| 19 | + Approximate the derivative of a function f(x) at a point x using the finite |
| 20 | + difference method |
| 21 | + |
| 22 | + >>> import math |
| 23 | + >>> tolerance = 1e-5 |
| 24 | + >>> derivative = calc_derivative(lambda x: x**2, 2) |
| 25 | + >>> math.isclose(derivative, 4, abs_tol=tolerance) |
| 26 | + True |
| 27 | + >>> derivative = calc_derivative(math.sin, 0) |
| 28 | + >>> math.isclose(derivative, 1, abs_tol=tolerance) |
| 29 | + True |
| 30 | + """ |
| 31 | + return (f(x + delta_x / 2) - f(x - delta_x / 2)) / delta_x |
| 32 | + |
| 33 | + |
| 34 | +def newton_raphson( |
| 35 | + f: RealFunc, |
| 36 | + x0: float = 0, |
| 37 | + max_iter: int = 100, |
| 38 | + step: float = 1e-6, |
| 39 | + max_error: float = 1e-6, |
| 40 | + log_steps: bool = False, |
| 41 | +) -> tuple[float, float, list[float]]: |
22 | 42 | """ |
23 | | - x = symbols("x") |
24 | | - f = lambdify(x, func, "math") |
25 | | - f_derivative = lambdify(x, diff(func), "math") |
26 | | - x_curr = a |
27 | | - while True: |
28 | | - x_curr = Decimal(x_curr) - Decimal(f(x_curr)) / Decimal(f_derivative(x_curr)) |
29 | | - if abs(f(x_curr)) < precision: |
30 | | - return float(x_curr) |
| 43 | + Find a root of the given function f using the Newton-Raphson method. |
| 44 | + |
| 45 | + :param f: A real-valued single-variable function |
| 46 | + :param x0: Initial guess |
| 47 | + :param max_iter: Maximum number of iterations |
| 48 | + :param step: Step size of x, used to approximate f'(x) |
| 49 | + :param max_error: Maximum approximation error |
| 50 | + :param log_steps: bool denoting whether to log intermediate steps |
| 51 | + |
| 52 | + :return: A tuple containing the approximation, the error, and the intermediate |
| 53 | + steps. If log_steps is False, then an empty list is returned for the third |
| 54 | + element of the tuple. |
| 55 | + |
| 56 | + :raises ZeroDivisionError: The derivative approaches 0. |
| 57 | + :raises ArithmeticError: No solution exists, or the solution isn't found before the |
| 58 | + iteration limit is reached. |
| 59 | + |
| 60 | + >>> import math |
| 61 | + >>> tolerance = 1e-15 |
| 62 | + >>> root, *_ = newton_raphson(lambda x: x**2 - 5*x + 2, 0.4, max_error=tolerance) |
| 63 | + >>> math.isclose(root, (5 - math.sqrt(17)) / 2, abs_tol=tolerance) |
| 64 | + True |
| 65 | + >>> root, *_ = newton_raphson(lambda x: math.log(x) - 1, 2, max_error=tolerance) |
| 66 | + >>> math.isclose(root, math.e, abs_tol=tolerance) |
| 67 | + True |
| 68 | + >>> root, *_ = newton_raphson(math.sin, 1, max_error=tolerance) |
| 69 | + >>> math.isclose(root, 0, abs_tol=tolerance) |
| 70 | + True |
| 71 | + >>> newton_raphson(math.cos, 0) |
| 72 | + Traceback (most recent call last): |
| 73 | + ... |
| 74 | + ZeroDivisionError: No converging solution found, zero derivative |
| 75 | + >>> newton_raphson(lambda x: x**2 + 1, 2) |
| 76 | + Traceback (most recent call last): |
| 77 | + ... |
| 78 | + ArithmeticError: No converging solution found, iteration limit reached |
| 79 | + """ |
| 80 | + |
| 81 | + def f_derivative(x: float) -> float: |
| 82 | + return calc_derivative(f, x, step) |
| 83 | + |
| 84 | + a = x0 # Set initial guess |
| 85 | + steps = [] |
| 86 | + for _ in range(max_iter): |
| 87 | + if log_steps: # Log intermediate steps |
| 88 | + steps.append(a) |
| 89 | + |
| 90 | + error = abs(f(a)) |
| 91 | + if error < max_error: |
| 92 | + return a, error, steps |
| 93 | + |
| 94 | + if f_derivative(a) == 0: |
| 95 | + raise ZeroDivisionError("No converging solution found, zero derivative") |
| 96 | + a -= f(a) / f_derivative(a) # Calculate next estimate |
| 97 | + raise ArithmeticError("No converging solution found, iteration limit reached") |
31 | 98 |
|
32 | 99 |
|
33 | 100 | if __name__ == "__main__": |
34 | 101 | import doctest |
| 102 | + from math import exp, tanh |
35 | 103 |
|
36 | 104 | doctest.testmod() |
37 | 105 |
|
38 | | - # Find value of pi |
39 | | - print(f"The root of sin(x) = 0 is {newton_raphson('sin(x)', 2)}") |
40 | | -# Find root of polynomial |
41 | | - print(f"The root of x**2 - 5*x + 2 = 0 is {newton_raphson('x**2 - 5*x + 2', 0.4)}") |
42 | | - # Find value of e |
43 | | - print(f"The root of log(x) - 1 = 0 is {newton_raphson('log(x) - 1', 2)}") |
44 | | - # Find root of exponential function |
45 | | - print(f"The root of exp(x) - 1 = 0 is {newton_raphson('exp(x) - 1', 0)}") |
| 106 | + deffunc(x: float) ->float: |
| 107 | + returntanh(x) **2-exp(3*x) |
| 108 | + |
| 109 | + solution, err, steps=newton_raphson( |
| 110 | + func, x0=10, max_iter=100, step=1e-6, log_steps=True |
| 111 | + ) |
| 112 | + print(f"{solution=}, {err=}") |
| 113 | + print("\n".join(str(x) forxinsteps)) |
0 commit comments