TheAlgorithms/Python
1
0
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases Wiki Activity

Delete arithmetic_analysis/ directory and relocate its contents (#10824)

Browse Source 
* Remove eval from arithmetic_analysis/newton_raphson.py

* Relocate contents of arithmetic_analysis/

Delete the arithmetic_analysis/ directory and relocate its files because
the purpose of the directory was always ill-defined. "Arithmetic
analysis" isn't a field of math, and the directory's files contained
algorithms for linear algebra, numerical analysis, and physics.

Relocated the directory's linear algebra algorithms to linear_algebra/,
its numerical analysis algorithms to a new subdirectory called
maths/numerical_analysis/, and its single physics algorithm to physics/.

* updating DIRECTORY.md

---------

Co-authored-by: github-actions <${GITHUB_ACTOR}@users.noreply.github.com>
This commit is contained in:
Tianyi Zheng
2023-10-23 03:31:30 -04:00
committed by GitHub
parent a9cee1d933
commit a8b6bda993
26 changed files with 335 additions and 344 deletions
Download Patch File Download Diff File Expand all files Collapse all files

55
maths/numerical_analysis/bisection.py Normal file
View File

@@ -0,0 +1,55 @@
from collections.abc import Callable
def bisection(function: Callable[[float], float], a: float, b: float) -> float:
"""
finds where function becomes 0 in [a,b] using bolzano
>>> bisection(lambda x: x ** 3 - 1, -5, 5)
1.0000000149011612
>>> bisection(lambda x: x ** 3 - 1, 2, 1000)
Traceback (most recent call last):
...
ValueError: could not find root in given interval.
>>> bisection(lambda x: x ** 2 - 4 * x + 3, 0, 2)
1.0
>>> bisection(lambda x: x ** 2 - 4 * x + 3, 2, 4)
3.0
>>> bisection(lambda x: x ** 2 - 4 * x + 3, 4, 1000)
Traceback (most recent call last):
...
ValueError: could not find root in given interval.
"""
start: float = a
end: float = b
if function(a) == 0: # one of the a or b is a root for the function
return a
elif function(b) == 0:
return b
elif (
function(a) * function(b) > 0
): # if none of these are root and they are both positive or negative,
# then this algorithm can't find the root
raise ValueError("could not find root in given interval.")
else:
mid: float = start + (end - start) / 2.0
while abs(start - mid) > 10**-7: # until precisely equals to 10^-7
if function(mid) == 0:
return mid
elif function(mid) * function(start) < 0:
end = mid
else:
start = mid
mid = start + (end - start) / 2.0
return mid
def f(x: float) -> float:
return x**3 - 2 * x - 5
if __name__ == "__main__":
print(bisection(f, 1, 1000))
import doctest
doctest.testmod()

63
maths/numerical_analysis/bisection_2.py Normal file
View File

@@ -0,0 +1,63 @@
"""
Given a function on floating number f(x) and two floating numbers ‘a’ and ‘b’ such that
f(a) * f(b) < 0 and f(x) is continuous in [a, b].
Here f(x) represents algebraic or transcendental equation.
Find root of function in interval [a, b] (Or find a value of x such that f(x) is 0)
https://en.wikipedia.org/wiki/Bisection_method
"""
def equation(x: float) -> float:
"""
>>> equation(5)
-15
>>> equation(0)
10
>>> equation(-5)
-15
>>> equation(0.1)
9.99
>>> equation(-0.1)
9.99
"""
return 10 - x * x
def bisection(a: float, b: float) -> float:
"""
>>> bisection(-2, 5)
3.1611328125
>>> bisection(0, 6)
3.158203125
>>> bisection(2, 3)
Traceback (most recent call last):
...
ValueError: Wrong space!
"""
# Bolzano theory in order to find if there is a root between a and b
if equation(a) * equation(b) >= 0:
raise ValueError("Wrong space!")
c = a
while (b - a) >= 0.01:
# Find middle point
c = (a + b) / 2
# Check if middle point is root
if equation(c) == 0.0:
break
# Decide the side to repeat the steps
if equation(c) * equation(a) < 0:
b = c
else:
a = c
return c
if __name__ == "__main__":
import doctest
doctest.testmod()
print(bisection(-2, 5))
print(bisection(0, 6))

121
maths/numerical_analysis/integration_by_simpson_approx.py Normal file
View File

@@ -0,0 +1,121 @@
"""
Author : Syed Faizan ( 3rd Year IIIT Pune )
Github : faizan2700
Purpose : You have one function f(x) which takes float integer and returns
float you have to integrate the function in limits a to b.
The approximation proposed by Thomas Simpsons in 1743 is one way to calculate
integration.
( read article : https://cp-algorithms.com/num_methods/simpson-integration.html )
simpson_integration() takes function,lower_limit=a,upper_limit=b,precision and
returns the integration of function in given limit.
"""
# constants
# the more the number of steps the more accurate
N_STEPS = 1000
def f(x: float) -> float:
return x * x
"""
Summary of Simpson Approximation :
By simpsons integration :
1. integration of fxdx with limit a to b is =
f(x0) + 4 * f(x1) + 2 * f(x2) + 4 * f(x3) + 2 * f(x4)..... + f(xn)
where x0 = a
xi = a + i * h
xn = b
"""
def simpson_integration(function, a: float, b: float, precision: int = 4) -> float:
"""
Args:
function : the function which's integration is desired
a : the lower limit of integration
b : upper limit of integration
precision : precision of the result,error required default is 4
Returns:
result : the value of the approximated integration of function in range a to b
Raises:
AssertionError: function is not callable
AssertionError: a is not float or integer
AssertionError: function should return float or integer
AssertionError: b is not float or integer
AssertionError: precision is not positive integer
>>> simpson_integration(lambda x : x*x,1,2,3)
2.333
>>> simpson_integration(lambda x : x*x,'wrong_input',2,3)
Traceback (most recent call last):
...
AssertionError: a should be float or integer your input : wrong_input
>>> simpson_integration(lambda x : x*x,1,'wrong_input',3)
Traceback (most recent call last):
...
AssertionError: b should be float or integer your input : wrong_input
>>> simpson_integration(lambda x : x*x,1,2,'wrong_input')
Traceback (most recent call last):
...
AssertionError: precision should be positive integer your input : wrong_input
>>> simpson_integration('wrong_input',2,3,4)
Traceback (most recent call last):
...
AssertionError: the function(object) passed should be callable your input : ...
>>> simpson_integration(lambda x : x*x,3.45,3.2,1)
-2.8
>>> simpson_integration(lambda x : x*x,3.45,3.2,0)
Traceback (most recent call last):
...
AssertionError: precision should be positive integer your input : 0
>>> simpson_integration(lambda x : x*x,3.45,3.2,-1)
Traceback (most recent call last):
...
AssertionError: precision should be positive integer your input : -1
"""
assert callable(
function
), f"the function(object) passed should be callable your input : {function}"
assert isinstance(a, (float, int)), f"a should be float or integer your input : {a}"
assert isinstance(function(a), (float, int)), (
"the function should return integer or float return type of your function, "
f"{type(a)}"
)
assert isinstance(b, (float, int)), f"b should be float or integer your input : {b}"
assert (
isinstance(precision, int) and precision > 0
), f"precision should be positive integer your input : {precision}"
# just applying the formula of simpson for approximate integration written in
# mentioned article in first comment of this file and above this function
h = (b - a) / N_STEPS
result = function(a) + function(b)
for i in range(1, N_STEPS):
a1 = a + h * i
result += function(a1) * (4 if i % 2 else 2)
result *= h / 3
return round(result, precision)
if __name__ == "__main__":
import doctest
doctest.testmod()

49
maths/numerical_analysis/intersection.py Normal file
View File

@@ -0,0 +1,49 @@
import math
from collections.abc import Callable
def intersection(function: Callable[[float], float], x0: float, x1: float) -> float:
"""
function is the f we want to find its root
x0 and x1 are two random starting points
>>> intersection(lambda x: x ** 3 - 1, -5, 5)
0.9999999999954654
>>> intersection(lambda x: x ** 3 - 1, 5, 5)
Traceback (most recent call last):
...
ZeroDivisionError: float division by zero, could not find root
>>> intersection(lambda x: x ** 3 - 1, 100, 200)
1.0000000000003888
>>> intersection(lambda x: x ** 2 - 4 * x + 3, 0, 2)
0.9999999998088019
>>> intersection(lambda x: x ** 2 - 4 * x + 3, 2, 4)
2.9999999998088023
>>> intersection(lambda x: x ** 2 - 4 * x + 3, 4, 1000)
3.0000000001786042
>>> intersection(math.sin, -math.pi, math.pi)
0.0
>>> intersection(math.cos, -math.pi, math.pi)
Traceback (most recent call last):
...
ZeroDivisionError: float division by zero, could not find root
"""
x_n: float = x0
x_n1: float = x1
while True:
if x_n == x_n1 or function(x_n1) == function(x_n):
raise ZeroDivisionError("float division by zero, could not find root")
x_n2: float = x_n1 - (
function(x_n1) / ((function(x_n1) - function(x_n)) / (x_n1 - x_n))
)
if abs(x_n2 - x_n1) < 10**-5:
return x_n2
x_n = x_n1
x_n1 = x_n2
def f(x: float) -> float:
return math.pow(x, 3) - (2 * x) - 5
if __name__ == "__main__":
print(intersection(f, 3, 3.5))

55
maths/numerical_analysis/nevilles_method.py Normal file
View File

@@ -0,0 +1,55 @@
"""
Python program to show how to interpolate and evaluate a polynomial
using Neville's method.
Neville’s method evaluates a polynomial that passes through a
given set of x and y points for a particular x value (x0) using the
Newton polynomial form.
Reference:
https://rpubs.com/aaronsc32/nevilles-method-polynomial-interpolation
"""
def neville_interpolate(x_points: list, y_points: list, x0: int) -> list:
"""
Interpolate and evaluate a polynomial using Neville's method.
Arguments:
x_points, y_points: Iterables of x and corresponding y points through
which the polynomial passes.
x0: The value of x to evaluate the polynomial for.
Return Value: A list of the approximated value and the Neville iterations
table respectively.
>>> import pprint
>>> neville_interpolate((1,2,3,4,6), (6,7,8,9,11), 5)[0]
10.0
>>> pprint.pprint(neville_interpolate((1,2,3,4,6), (6,7,8,9,11), 99)[1])
[[0, 6, 0, 0, 0],
[0, 7, 0, 0, 0],
[0, 8, 104.0, 0, 0],
[0, 9, 104.0, 104.0, 0],
[0, 11, 104.0, 104.0, 104.0]]
>>> neville_interpolate((1,2,3,4,6), (6,7,8,9,11), 99)[0]
104.0
>>> neville_interpolate((1,2,3,4,6), (6,7,8,9,11), '')
Traceback (most recent call last):
...
TypeError: unsupported operand type(s) for -: 'str' and 'int'
"""
n = len(x_points)
q = [[0] * n for i in range(n)]
for i in range(n):
q[i][1] = y_points[i]
for i in range(2, n):
for j in range(i, n):
q[j][i] = (
(x0 - x_points[j - i + 1]) * q[j][i - 1]
- (x0 - x_points[j]) * q[j - 1][i - 1]
) / (x_points[j] - x_points[j - i + 1])
return [q[n - 1][n - 1], q]
if __name__ == "__main__":
import doctest
doctest.testmod()

57
maths/numerical_analysis/newton_forward_interpolation.py Normal file
View File

@@ -0,0 +1,57 @@
# https://www.geeksforgeeks.org/newton-forward-backward-interpolation/
from __future__ import annotations
import math
# for calculating u value
def ucal(u: float, p: int) -> float:
"""
>>> ucal(1, 2)
0
>>> ucal(1.1, 2)
0.11000000000000011
>>> ucal(1.2, 2)
0.23999999999999994
"""
temp = u
for i in range(1, p):
temp = temp * (u - i)
return temp
def main() -> None:
n = int(input("enter the numbers of values: "))
y: list[list[float]] = []
for _ in range(n):
y.append([])
for i in range(n):
for j in range(n):
y[i].append(j)
y[i][j] = 0
print("enter the values of parameters in a list: ")
x = list(map(int, input().split()))
print("enter the values of corresponding parameters: ")
for i in range(n):
y[i][0] = float(input())
value = int(input("enter the value to interpolate: "))
u = (value - x[0]) / (x[1] - x[0])
# for calculating forward difference table
for i in range(1, n):
for j in range(n - i):
y[j][i] = y[j + 1][i - 1] - y[j][i - 1]
summ = y[0][0]
for i in range(1, n):
summ += (ucal(u, i) * y[0][i]) / math.factorial(i)
print(f"the value at {value} is {summ}")
if __name__ == "__main__":
main()

54
maths/numerical_analysis/newton_method.py Normal file
View File

@@ -0,0 +1,54 @@
"""Newton's Method."""
# Newton's Method - https://en.wikipedia.org/wiki/Newton%27s_method
from collections.abc import Callable
RealFunc = Callable[[float], float] # type alias for a real -> real function
# function is the f(x) and derivative is the f'(x)
def newton(
function: RealFunc,
derivative: RealFunc,
starting_int: int,
) -> float:
"""
>>> newton(lambda x: x ** 3 - 2 * x - 5, lambda x: 3 * x ** 2 - 2, 3)
2.0945514815423474
>>> newton(lambda x: x ** 3 - 1, lambda x: 3 * x ** 2, -2)
1.0
>>> newton(lambda x: x ** 3 - 1, lambda x: 3 * x ** 2, -4)
1.0000000000000102
>>> import math
>>> newton(math.sin, math.cos, 1)
0.0
>>> newton(math.sin, math.cos, 2)
3.141592653589793
>>> newton(math.cos, lambda x: -math.sin(x), 2)
1.5707963267948966
>>> newton(math.cos, lambda x: -math.sin(x), 0)
Traceback (most recent call last):
...
ZeroDivisionError: Could not find root
"""
prev_guess = float(starting_int)
while True:
try:
next_guess = prev_guess - function(prev_guess) / derivative(prev_guess)
except ZeroDivisionError:
raise ZeroDivisionError("Could not find root") from None
if abs(prev_guess - next_guess) < 10**-5:
return next_guess
prev_guess = next_guess
def f(x: float) -> float:
return (x**3) - (2 * x) - 5
def f1(x: float) -> float:
return 3 * (x**2) - 2
if __name__ == "__main__":
print(newton(f, f1, 3))

45
maths/numerical_analysis/newton_raphson.py Normal file
View File

@@ -0,0 +1,45 @@
# 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
"""
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)
if __name__ == "__main__":
import doctest
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)}")

64
maths/numerical_analysis/newton_raphson_2.py Normal file
View File

@@ -0,0 +1,64 @@
"""
Author: P Shreyas Shetty
Implementation of Newton-Raphson method for solving equations of kind
f(x) = 0. It is an iterative method where solution is found by the expression
x[n+1] = x[n] + f(x[n])/f'(x[n])
If no solution exists, then either the solution will not be found when iteration
limit is reached or the gradient f'(x[n]) approaches zero. In both cases, exception
is raised. If iteration limit is reached, try increasing maxiter.
"""
import math as m
from collections.abc import Callable
DerivativeFunc = Callable[[float], float]
def calc_derivative(f: DerivativeFunc, a: float, h: float = 0.001) -> float:
"""
Calculates derivative at point a for function f using finite difference
method
"""
return (f(a + h) - f(a - h)) / (2 * h)
def newton_raphson(
f: DerivativeFunc,
x0: float = 0,
maxiter: int = 100,
step: float = 0.0001,
maxerror: float = 1e-6,
logsteps: bool = False,
) -> tuple[float, float, list[float]]:
a = x0 # set the initial guess
steps = [a]
error = abs(f(a))
f1 = lambda x: calc_derivative(f, x, h=step) # noqa: E731 Derivative of f(x)
for _ in range(maxiter):
if f1(a) == 0:
raise ValueError("No converging solution found")
a = a - f(a) / f1(a) # Calculate the next estimate
if logsteps:
steps.append(a)
if error < maxerror:
break
else:
raise ValueError("Iteration limit reached, no converging solution found")
if logsteps:
# If logstep is true, then log intermediate steps
return a, error, steps
return a, error, []
if __name__ == "__main__":
from matplotlib import pyplot as plt
f = lambda x: m.tanh(x) ** 2 - m.exp(3 * x) # noqa: E731
solution, error, steps = newton_raphson(
f, x0=10, maxiter=1000, step=1e-6, logsteps=True
)
plt.plot([abs(f(x)) for x in steps])
plt.xlabel("step")
plt.ylabel("error")
plt.show()
print(f"solution = {{{solution:f}}}, error = {{{error:f}}}")

83
maths/numerical_analysis/newton_raphson_new.py Normal file
View File

@@ -0,0 +1,83 @@
# Implementing Newton Raphson method in Python
# Author: Saksham Gupta
#
# The Newton-Raphson method (also known as Newton's method) is a way to
# quickly find a good approximation for the root of a functreal-valued ion
# The method can also be extended to complex functions
#
# Newton's Method - https://en.wikipedia.org/wiki/Newton's_method
from sympy import diff, lambdify, symbols
from sympy.functions import * # noqa: F403
def newton_raphson(
function: str,
starting_point: complex,
variable: str = "x",
precision: float = 10**-10,
multiplicity: int = 1,
) -> complex:
"""Finds root from the 'starting_point' onwards by Newton-Raphson method
Refer to https://docs.sympy.org/latest/modules/functions/index.html
for usable mathematical functions
>>> newton_raphson("sin(x)", 2)
3.141592653589793
>>> newton_raphson("x**4 -5", 0.4 + 5j)
(-7.52316384526264e-37+1.4953487812212207j)
>>> newton_raphson('log(y) - 1', 2, variable='y')
2.7182818284590455
>>> newton_raphson('exp(x) - 1', 10, precision=0.005)
1.2186556186174883e-10
>>> newton_raphson('cos(x)', 0)
Traceback (most recent call last):
...
ZeroDivisionError: Could not find root
"""
x = symbols(variable)
func = lambdify(x, function)
diff_function = lambdify(x, diff(function, x))
prev_guess = starting_point
while True:
if diff_function(prev_guess) != 0:
next_guess = prev_guess - multiplicity * func(prev_guess) / diff_function(
prev_guess
)
else:
raise ZeroDivisionError("Could not find root") from None
# Precision is checked by comparing the difference of consecutive guesses
if abs(next_guess - prev_guess) < precision:
return next_guess
prev_guess = next_guess
# Let's Execute
if __name__ == "__main__":
# Find root of trigonometric function
# Find value of pi
print(f"The root of sin(x) = 0 is {newton_raphson('sin(x)', 2)}")
# Find root of polynomial
# Find fourth Root of 5
print(f"The root of x**4 - 5 = 0 is {newton_raphson('x**4 -5', 0.4 +5j)}")
# Find value of e
print(
"The root of log(y) - 1 = 0 is ",
f"{newton_raphson('log(y) - 1', 2, variable='y')}",
)
# Exponential Roots
print(
"The root of exp(x) - 1 = 0 is",
f"{newton_raphson('exp(x) - 1', 10, precision=0.005)}",
)
# Find root of cos(x)
print(f"The root of cos(x) = 0 is {newton_raphson('cos(x)', 0)}")

65
maths/numerical_analysis/numerical_integration.py Normal file
View File

@@ -0,0 +1,65 @@
"""
Approximates the area under the curve using the trapezoidal rule
"""
from __future__ import annotations
from collections.abc import Callable
def trapezoidal_area(
fnc: Callable[[float], float],
x_start: float,
x_end: float,
steps: int = 100,
) -> float:
"""
Treats curve as a collection of linear lines and sums the area of the
trapezium shape they form
:param fnc: a function which defines a curve
:param x_start: left end point to indicate the start of line segment
:param x_end: right end point to indicate end of line segment
:param steps: an accuracy gauge; more steps increases the accuracy
:return: a float representing the length of the curve
>>> def f(x):
... return 5
>>> '%.3f' % trapezoidal_area(f, 12.0, 14.0, 1000)
'10.000'
>>> def f(x):
... return 9*x**2
>>> '%.4f' % trapezoidal_area(f, -4.0, 0, 10000)
'192.0000'
>>> '%.4f' % trapezoidal_area(f, -4.0, 4.0, 10000)
'384.0000'
"""
x1 = x_start
fx1 = fnc(x_start)
area = 0.0
for _ in range(steps):
# Approximates small segments of curve as linear and solve
# for trapezoidal area
x2 = (x_end - x_start) / steps + x1
fx2 = fnc(x2)
area += abs(fx2 + fx1) * (x2 - x1) / 2
# Increment step
x1 = x2
fx1 = fx2
return area
if __name__ == "__main__":
def f(x):
return x**3
print("f(x) = x^3")
print("The area between the curve, x = -10, x = 10 and the x axis is:")
i = 10
while i <= 100000:
area = trapezoidal_area(f, -5, 5, i)
print(f"with {i} steps: {area}")
i *= 10

44
maths/numerical_analysis/runge_kutta.py Normal file
View File

@@ -0,0 +1,44 @@
import numpy as np
def runge_kutta(f, y0, x0, h, x_end):
"""
Calculate the numeric solution at each step to the ODE f(x, y) using RK4
https://en.wikipedia.org/wiki/Runge-Kutta_methods
Arguments:
f -- The ode as a function of x and y
y0 -- the initial value for y
x0 -- the initial value for x
h -- the stepsize
x_end -- the end value for x
>>> # the exact solution is math.exp(x)
>>> def f(x, y):
... return y
>>> y0 = 1
>>> y = runge_kutta(f, y0, 0.0, 0.01, 5)
>>> y[-1]
148.41315904125113
"""
n = int(np.ceil((x_end - x0) / h))
y = np.zeros((n + 1,))
y[0] = y0
x = x0
for k in range(n):
k1 = f(x, y[k])
k2 = f(x + 0.5 * h, y[k] + 0.5 * h * k1)
k3 = f(x + 0.5 * h, y[k] + 0.5 * h * k2)
k4 = f(x + h, y[k] + h * k3)
y[k + 1] = y[k] + (1 / 6) * h * (k1 + 2 * k2 + 2 * k3 + k4)
x += h
return y
if __name__ == "__main__":
import doctest
doctest.testmod()

114
maths/numerical_analysis/runge_kutta_fehlberg_45.py Normal file
View File

@@ -0,0 +1,114 @@
"""
Use the Runge-Kutta-Fehlberg method to solve Ordinary Differential Equations.
"""
from collections.abc import Callable
import numpy as np
def runge_kutta_fehlberg_45(
func: Callable,
x_initial: float,
y_initial: float,
step_size: float,
x_final: float,
) -> np.ndarray:
"""
Solve an Ordinary Differential Equations using Runge-Kutta-Fehlberg Method (rkf45)
of order 5.
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
args:
func: An ordinary differential equation (ODE) as function of x and y.
x_initial: The initial value of x.
y_initial: The initial value of y.
step_size: The increment value of x.
x_final: The final value of x.
Returns:
Solution of y at each nodal point
# exact value of y[1] is tan(0.2) = 0.2027100937470787
>>> def f(x, y):
... return 1 + y**2
>>> y = runge_kutta_fehlberg_45(f, 0, 0, 0.2, 1)
>>> y[1]
0.2027100937470787
>>> def f(x,y):
... return x
>>> y = runge_kutta_fehlberg_45(f, -1, 0, 0.2, 0)
>>> y[1]
-0.18000000000000002
>>> y = runge_kutta_fehlberg_45(5, 0, 0, 0.1, 1)
Traceback (most recent call last):
...
TypeError: 'int' object is not callable
>>> def f(x, y):
... return x + y
>>> y = runge_kutta_fehlberg_45(f, 0, 0, 0.2, -1)
Traceback (most recent call last):
...
ValueError: The final value of x must be greater than initial value of x.
>>> def f(x, y):
... return x
>>> y = runge_kutta_fehlberg_45(f, -1, 0, -0.2, 0)
Traceback (most recent call last):
...
ValueError: Step size must be positive.
"""
if x_initial >= x_final:
raise ValueError(
"The final value of x must be greater than initial value of x."
)
if step_size <= 0:
raise ValueError("Step size must be positive.")
n = int((x_final - x_initial) / step_size)
y = np.zeros(
(n + 1),
)
x = np.zeros(n + 1)
y[0] = y_initial
x[0] = x_initial
for i in range(n):
k1 = step_size * func(x[i], y[i])
k2 = step_size * func(x[i] + step_size / 4, y[i] + k1 / 4)
k3 = step_size * func(
x[i] + (3 / 8) * step_size, y[i] + (3 / 32) * k1 + (9 / 32) * k2
)
k4 = step_size * func(
x[i] + (12 / 13) * step_size,
y[i] + (1932 / 2197) * k1 - (7200 / 2197) * k2 + (7296 / 2197) * k3,
)
k5 = step_size * func(
x[i] + step_size,
y[i] + (439 / 216) * k1 - 8 * k2 + (3680 / 513) * k3 - (845 / 4104) * k4,
)
k6 = step_size * func(
x[i] + step_size / 2,
y[i]
- (8 / 27) * k1
+ 2 * k2
- (3544 / 2565) * k3
+ (1859 / 4104) * k4
- (11 / 40) * k5,
)
y[i + 1] = (
y[i]
+ (16 / 135) * k1
+ (6656 / 12825) * k3
+ (28561 / 56430) * k4
- (9 / 50) * k5
+ (2 / 55) * k6
)
x[i + 1] = step_size + x[i]
return y
if __name__ == "__main__":
import doctest
doctest.testmod()

29
maths/numerical_analysis/secant_method.py Normal file
View File

@@ -0,0 +1,29 @@
"""
Implementing Secant method in Python
Author: dimgrichr
"""
from math import exp
def f(x: float) -> float:
"""
>>> f(5)
39.98652410600183
"""
return 8 * x - 2 * exp(-x)
def secant_method(lower_bound: float, upper_bound: float, repeats: int) -> float:
"""
>>> secant_method(1, 3, 2)
0.2139409276214589
"""
x0 = lower_bound
x1 = upper_bound
for _ in range(repeats):
x0, x1 = x1, x1 - (f(x1) * (x1 - x0)) / (f(x1) - f(x0))
return x1
if __name__ == "__main__":
print(f"Example: {secant_method(1, 3, 2)}")

86
maths/numerical_analysis/simpson_rule.py Normal file
View File

@@ -0,0 +1,86 @@
"""
Numerical integration or quadrature for a smooth function f with known values at x_i
This method is the classical approach of summing 'Equally Spaced Abscissas'
method 2:
"Simpson Rule"
"""
def method_2(boundary: list[int], steps: int) -> float:
# "Simpson Rule"
# int(f) = delta_x/2 * (b-a)/3*(f1 + 4f2 + 2f_3 + ... + fn)
"""
Calculate the definite integral of a function using Simpson's Rule.
:param boundary: A list containing the lower and upper bounds of integration.
:param steps: The number of steps or resolution for the integration.
:return: The approximate integral value.
>>> round(method_2([0, 2, 4], 10), 10)
2.6666666667
>>> round(method_2([2, 0], 10), 10)
-0.2666666667
>>> round(method_2([-2, -1], 10), 10)
2.172
>>> round(method_2([0, 1], 10), 10)
0.3333333333
>>> round(method_2([0, 2], 10), 10)
2.6666666667
>>> round(method_2([0, 2], 100), 10)
2.5621226667
>>> round(method_2([0, 1], 1000), 10)
0.3320026653
>>> round(method_2([0, 2], 0), 10)
Traceback (most recent call last):
...
ZeroDivisionError: Number of steps must be greater than zero
>>> round(method_2([0, 2], -10), 10)
Traceback (most recent call last):
...
ZeroDivisionError: Number of steps must be greater than zero
"""
if steps <= 0:
raise ZeroDivisionError("Number of steps must be greater than zero")
h = (boundary[1] - boundary[0]) / steps
a = boundary[0]
b = boundary[1]
x_i = make_points(a, b, h)
y = 0.0
y += (h / 3.0) * f(a)
cnt = 2
for i in x_i:
y += (h / 3) * (4 - 2 * (cnt % 2)) * f(i)
cnt += 1
y += (h / 3.0) * f(b)
return y
def make_points(a, b, h):
x = a + h
while x < (b - h):
yield x
x = x + h
def f(x): # enter your function here
y = (x - 0) * (x - 0)
return y
def main():
a = 0.0 # Lower bound of integration
b = 1.0 # Upper bound of integration
steps = 10.0 # number of steps or resolution
boundary = [a, b] # boundary of integration
y = method_2(boundary, steps)
print(f"y = {y}")
if __name__ == "__main__":
import doctest
doctest.testmod()
main()

63
maths/numerical_analysis/square_root.py Normal file
View File

@@ -0,0 +1,63 @@
import math
def fx(x: float, a: float) -> float:
return math.pow(x, 2) - a
def fx_derivative(x: float) -> float:
return 2 * x
def get_initial_point(a: float) -> float:
start = 2.0
while start <= a:
start = math.pow(start, 2)
return start
def square_root_iterative(
a: float, max_iter: int = 9999, tolerance: float = 1e-14
) -> float:
"""
Square root approximated using Newton's method.
https://en.wikipedia.org/wiki/Newton%27s_method
>>> all(abs(square_root_iterative(i) - math.sqrt(i)) <= 1e-14 for i in range(500))
True
>>> square_root_iterative(-1)
Traceback (most recent call last):
...
ValueError: math domain error
>>> square_root_iterative(4)
2.0
>>> square_root_iterative(3.2)
1.788854381999832
>>> square_root_iterative(140)
11.832159566199232
"""
if a < 0:
raise ValueError("math domain error")
value = get_initial_point(a)
for _ in range(max_iter):
prev_value = value
value = value - fx(value, a) / fx_derivative(value)
if abs(prev_value - value) < tolerance:
return value
return value
if __name__ == "__main__":
from doctest import testmod
testmod()