tessl/pypi-mpmath
Python library for arbitrary-precision floating-point arithmetic
—
Pending
numerical-calculus.mddocs/
Numerical Calculus
Numerical integration, differentiation, root finding, optimization, series summation, and differential equation solving with high precision.
Capabilities
Numerical Integration
Adaptive quadrature methods for computing definite integrals.
def quad(f, a, b, **kwargs):
"""
Adaptive quadrature integration.
Args:
f: Function to integrate
a, b: Integration limits (can be infinite)
**kwargs: Additional options (maxdegree, method, etc.)
Returns:
Definite integral ∫_a^b f(x) dx
"""
def quadgl(f, a, b, n):
"""
Gauss-Legendre quadrature.
Args:
f: Function to integrate
a, b: Integration limits
n: Number of quadrature points
Returns:
Approximation to ∫_a^b f(x) dx
"""
def quadts(f, a, b, **kwargs):
"""
Tanh-sinh (doubly exponential) quadrature.
Args:
f: Function to integrate
a, b: Integration limits
**kwargs: Additional options
Returns:
Integral using tanh-sinh transformation
"""
def quadosc(f, a, b, omega, **kwargs):
"""
Integration of oscillatory functions.
Args:
f: Function to integrate
a, b: Integration limits
omega: Oscillation frequency
**kwargs: Additional options
Returns:
Integral of f(x)*cos(omega*x) or f(x)*sin(omega*x)
"""
def quadsubdiv(f, a, b, **kwargs):
"""
Subdivision-based quadrature.
Args:
f: Function to integrate
a, b: Integration limits
**kwargs: Additional options
Returns:
Integral using adaptive subdivision
"""
def gauss_quadrature(n, a=-1, b=1):
"""
Generate Gauss-Legendre quadrature rule.
Args:
n: Number of quadrature points
a: Lower limit (default: -1)
b: Upper limit (default: 1)
Returns:
tuple: (nodes, weights) for quadrature
"""Inverse Laplace Transform
Methods for computing inverse Laplace transforms.
def invertlaplace(f, t, **kwargs):
"""
Numerical inverse Laplace transform.
Args:
f: Function in s-domain
t: Time value or array
**kwargs: Algorithm options
Returns:
f^(-1)(t) - inverse Laplace transform
"""
def invlaptalbot(f, t, **kwargs):
"""
Inverse Laplace transform using Talbot's method.
Args:
f: Function in s-domain
t: Time value
**kwargs: Algorithm parameters
Returns:
Inverse transform using Talbot contour
"""
def invlapstehfest(f, t, **kwargs):
"""
Inverse Laplace transform using Stehfest's method.
Args:
f: Function in s-domain
t: Time value
**kwargs: Algorithm parameters
Returns:
Inverse transform using Stehfest algorithm
"""
def invlapdehoog(f, t, **kwargs):
"""
Inverse Laplace transform using de Hoog's method.
Args:
f: Function in s-domain
t: Time value
**kwargs: Algorithm parameters
Returns:
Inverse transform using de Hoog method
"""Numerical Differentiation
Compute derivatives numerically with high accuracy.
def diff(f, x, n=1, **kwargs):
"""
Numerical differentiation.
Args:
f: Function to differentiate
x: Point at which to evaluate derivative
n: Order of derivative (default: 1)
**kwargs: Algorithm options (h, method, etc.)
Returns:
nth derivative f^(n)(x)
"""
def diffs(f, x, n):
"""
Compute multiple derivatives simultaneously.
Args:
f: Function to differentiate
x: Evaluation point
n: Maximum derivative order
Returns:
List [f(x), f'(x), f''(x), ..., f^(n)(x)]
"""
def diffs_prod(f, g, x, n):
"""
Derivatives of product using Leibniz rule.
Args:
f, g: Functions
x: Evaluation point
n: Maximum derivative order
Returns:
Derivatives of f(x)*g(x)
"""
def diffs_exp(f, x, n):
"""
Derivatives of exponential composition.
Args:
f: Function in exponent
x: Evaluation point
n: Maximum derivative order
Returns:
Derivatives of exp(f(x))
"""
def diffun(f, n=1):
"""
Create derivative function.
Args:
f: Function to differentiate
n: Derivative order
Returns:
Function computing nth derivative
"""
def differint(f, x, n):
"""
Fractional derivative/integral.
Args:
f: Function
x: Evaluation point
n: Fractional order (positive for derivative, negative for integral)
Returns:
Fractional derivative/integral of order n
"""
def difference(f, n=1):
"""
Finite difference operator.
Args:
f: Sequence or function
n: Order of difference
Returns:
nth finite difference
"""Series Expansions and Approximations
Taylor series, Padé approximants, and other expansions.
def taylor(f, x, n, **kwargs):
"""
Taylor series expansion.
Args:
f: Function to expand
x: Expansion point
n: Order of expansion
**kwargs: Additional options
Returns:
Coefficients of Taylor series
"""
def pade(p, q, **kwargs):
"""
Padé approximant from series coefficients.
Args:
p: Numerator degree
q: Denominator degree
**kwargs: Series coefficients or function
Returns:
Padé approximant coefficients
"""
def fourier(f, a, b, n):
"""
Fourier series coefficients.
Args:
f: Function to expand
a, b: Interval [a, b]
n: Number of coefficients
Returns:
Fourier series coefficients
"""
def fourierval(coeffs, x):
"""
Evaluate Fourier series.
Args:
coeffs: Fourier coefficients
x: Evaluation point
Returns:
Value of Fourier series at x
"""
def chebyfit(f, a, b, n):
"""
Chebyshev polynomial approximation.
Args:
f: Function to approximate
a, b: Interval [a, b]
n: Degree of approximation
Returns:
Chebyshev coefficients
"""Polynomial Operations
Operations on polynomials and polynomial root finding.
def polyval(coeffs, x):
"""
Evaluate polynomial at given point.
Args:
coeffs: Polynomial coefficients (highest degree first)
x: Evaluation point
Returns:
Value of polynomial at x
"""
def polyroots(coeffs):
"""
Find roots of polynomial.
Args:
coeffs: Polynomial coefficients
Returns:
List of polynomial roots
"""Series Summation and Extrapolation
Acceleration methods for infinite series and sequences.
def nsum(f, n1, n2=None, **kwargs):
"""
Numerical summation with acceleration.
Args:
f: Function or sequence to sum
n1: Starting index (or ending index if n2 is None)
n2: Ending index (optional, for infinite sums)
**kwargs: Acceleration method options
Returns:
Sum with acceleration/extrapolation
"""
def nprod(f, n1, n2=None, **kwargs):
"""
Numerical product with acceleration.
Args:
f: Function or sequence
n1: Starting index
n2: Ending index (optional)
**kwargs: Acceleration options
Returns:
Product with acceleration
"""
def richardson(sequence, **kwargs):
"""
Richardson extrapolation.
Args:
sequence: Sequence to extrapolate
**kwargs: Extrapolation parameters
Returns:
Extrapolated limit
"""
def shanks(sequence):
"""
Shanks transformation for sequence acceleration.
Args:
sequence: Sequence to accelerate
Returns:
Accelerated sequence
"""
def levin(sequence, **kwargs):
"""
Levin transformation.
Args:
sequence: Sequence to transform
**kwargs: Transformation parameters
Returns:
Transformed sequence
"""
def cohen_alt(sequence):
"""
Cohen alternating series acceleration.
Args:
sequence: Alternating series
Returns:
Accelerated sum
"""
def sumem(f, a, b, **kwargs):
"""
Euler-Maclaurin summation formula.
Args:
f: Function to sum
a, b: Summation range
**kwargs: Method parameters
Returns:
Sum using Euler-Maclaurin formula
"""
def sumap(f, a, b, **kwargs):
"""
Abel-Plana summation formula.
Args:
f: Function to sum
a, b: Summation range
**kwargs: Method parameters
Returns:
Sum using Abel-Plana formula
"""Root Finding and Optimization
Find zeros and extrema of functions.
def findroot(f, x0, **kwargs):
"""
Find root of equation f(x) = 0.
Args:
f: Function or system of equations
x0: Initial guess
**kwargs: Solver options (method, tol, etc.)
Returns:
Root of equation
"""
def multiplicity(f, x0):
"""
Determine multiplicity of root.
Args:
f: Function
x0: Root location
Returns:
Multiplicity of root at x0
"""
def jacobian(f, x, **kwargs):
"""
Compute Jacobian matrix of vector function.
Args:
f: Vector function
x: Evaluation point
**kwargs: Differentiation options
Returns:
Jacobian matrix
"""Differential Equations
Solve ordinary differential equations.
def odefun(f, x0, y0, **kwargs):
"""
Solve ODE system y' = f(x, y).
Args:
f: Right-hand side function f(x, y)
x0: Initial x value
y0: Initial condition y(x0)
**kwargs: Solver options
Returns:
ODE solution function
"""Limits
Compute limits of sequences and functions.
def limit(f, x, direction='+', **kwargs):
"""
Numerical limit computation.
Args:
f: Function
x: Limit point
direction: Approach direction ('+', '-', or '+-')
**kwargs: Algorithm options
Returns:
lim_{t→x} f(t)
"""Usage Examples
import mpmath
from mpmath import mp
# Set precision
mp.dps = 25
# Numerical integration
def f(x):
return mp.exp(-x**2)
# Gaussian integral from 0 to infinity
integral = mp.quad(f, [0, mp.inf])
print(f"∫₀^∞ e^(-x²) dx = {integral}")
print(f"√π/2 = {mp.sqrt(mp.pi)/2}")
# Different quadrature methods
integral_gl = mp.quadgl(f, 0, 2, 20) # Gauss-Legendre
integral_ts = mp.quadts(f, 0, 2) # Tanh-sinh
print(f"Gauss-Legendre: {integral_gl}")
print(f"Tanh-sinh: {integral_ts}")
# Numerical differentiation
def g(x):
return mp.sin(x**2)
derivative = mp.diff(g, 1) # g'(1)
print(f"d/dx sin(x²) at x=1: {derivative}")
print(f"2*cos(1)*1 = {2*mp.cos(1)*1}") # Analytical result
# Multiple derivatives
derivs = mp.diffs(mp.exp, 0, 5) # exp and its derivatives at x=0
print(f"exp and derivatives at 0: {derivs}")
# Root finding
def h(x):
return x**3 - 2*x - 5
root = mp.findroot(h, 2) # Initial guess x=2
print(f"Root of x³ - 2x - 5 = 0: {root}")
print(f"Verification h({root}) = {h(root)}")
# Series summation with acceleration
def term(n):
return (-1)**(n-1) / n # Alternating harmonic series
# Sum first 100 terms with acceleration
partial_sum = mp.nsum(term, [1, 100])
print(f"Accelerated sum: {partial_sum}")
print(f"ln(2) = {mp.ln(2)}") # Should converge to ln(2)
# Taylor series
taylor_coeffs = mp.taylor(mp.exp, 0, 10) # exp(x) around x=0
print(f"Taylor coefficients of exp(x): {taylor_coeffs[:5]}")
# Polynomial evaluation and roots
poly_coeffs = [1, 0, -1] # x² - 1
roots = mp.polyroots(poly_coeffs)
print(f"Roots of x² - 1: {roots}")
# Evaluate polynomial
value = mp.polyval(poly_coeffs, 2) # 2² - 1 = 3
print(f"Value of x² - 1 at x=2: {value}")
# Limits
def limit_func(x):
return mp.sin(x) / x
limit_val = mp.limit(limit_func, 0)
print(f"lim_{x→0} sin(x)/x = {limit_val}")
# ODE solving
def ode_rhs(x, y):
return -y # y' = -y, solution is y = y0 * exp(-x)
ode_solution = mp.odefun(ode_rhs, 0, 1) # y(0) = 1
y_at_1 = ode_solution(1)
print(f"ODE solution at x=1: {y_at_1}")
print(f"Analytical: exp(-1) = {mp.exp(-1)}")
# Inverse Laplace transform
def laplace_func(s):
return 1 / (s + 1) # L^(-1){1/(s+1)} = exp(-t)
time_val = 1
inverse = mp.invertlaplace(laplace_func, time_val)
print(f"Inverse Laplace transform at t=1: {inverse}")
print(f"Expected exp(-1): {mp.exp(-1)}")