Numerical

A clean, elegant Python library for scientific computing and numerical methods.

Features

Numerical provides four main domains of numerical computation:

Root Finding

• Bisection - Robust bracketing method
• False Position - Improved bisection with linear interpolation
• Newton-Raphson - Fast convergence with automatic differentiation

Numerical Integration

• Riemann Sum - Simple rectangle rule
• Trapezoid Rule - Piecewise linear approximation
• Simpson's Rule - Quadratic approximation with even subdivision requirement

Numerical Differentiation

• Finite Differences - Central, forward, and backward differences
• Automatic numerical derivative computation

ODE Solving

• Euler Method - First-order explicit method
• Runge-Kutta 4th Order - Fourth-order explicit method

Design Principles

• Semantic Coherence - Clear, intuitive organization by mathematical domain
• Minimal Abstraction - No unnecessary complexity or overengineering
• Low Coupling - Modules work independently
• Pythonic API - Simple, predictable interfaces
• Type Safety - Complete type hints throughout
• Zero Dependencies - Core library has no external dependencies

Installation

From Source

git clone https://github.com/MiguelMochizukiDev/numerical.git
cd numerical
python -m pip install -e .

Requirements

• Python 3.12 or later

Quick Start

Root Finding

from numerical.roots import Bisection, Newton, FalsePosition

# Define your function
def f(x):
    return x**3 - 2

# Create a solver
solver = Bisection(f, a=0, b=2, tol=1e-6)

# Get the solution
result = solver.solve()
print(f"Root: {result.root}")
print(f"Residual: {result.residual}")
print(f"Iterations: {result.iterations}")

Iterate through steps:

solver = Newton(f, a=0, b=2)

for step in solver.steps():
    print(f"[{step.iteration}] x = {step.approximation:.6f}, "
          f"f(x) = {step.residual:.2e}")

Numerical Integration

from numerical.integration import Simpson, Trapezoid

def g(x):
    return x**2

# Simpson's Rule requires even number of subintervals
integrator = Simpson(g, a=0, b=1, n=100)
result = integrator.compute()
print(f"Integral: {result}")

# Trapezoid rule (no restrictions)
integrator = Trapezoid(g, a=0, b=1, n=100)
result = integrator.compute()
print(f"Integral: {result}")

Numerical Differentiation

from numerical.differentiation import finite_difference

def f(x):
    return x**2

# Get derivative function
df = finite_difference(f, h=1e-5, method="central")

# Evaluate at point
print(f"f'(2) = {df(2)}")

# Available methods: "central", "forward", "backward"

ODE Solving

from numerical.ode import Euler, RungeKutta

def dydt(t, y):
    return -2*y + 1

# Euler method
solver = Euler(dydt, t0=0, y0=1, tf=2, n=100)
t, y = solver.solve()

# Runge-Kutta 4th order
solver = RungeKutta(dydt, t0=0, y0=1, tf=2, n=100)
t, y = solver.solve()

API Reference

Root Finding

RootFinder (Abstract Base Class)

solver = Bisection(f, a, b, tol=1e-6, max_iter=1000)

# Methods
result = solver.solve()  # -> RootSolution
for step in solver.steps():  # -> Iterator[IterationResult]
    ...

Parameters:

• f: ScalarFunction - Function for which to find root
• a, b: float - Initial interval bounds (or initial guess for Newton)
• tol: float - Convergence tolerance (default: 1e-6)
• max_iter: int - Maximum iterations (default: 1000)

Returns:

• RootSolution - Final result with root, residual, iterations, convergence status

Specific Methods

• Bisection: Requires sign change in interval
• False Position: Linear interpolation-based bracketing
• Newton: Requires derivative (computed automatically if not provided)

Integration

Integrator (Abstract Base Class)

integrator = Simpson(f, a, b, n)
result = integrator.compute()  # -> float

Parameters:

• f: ScalarFunction - Function to integrate
• a, b: float - Integration bounds (a < b)
• n: int - Number of subintervals (must be even for Simpson)

Specific Methods

• Riemann: Any positive n
• Trapezoid: Any positive n
• Simpson: Requires even n (raises ValidationError otherwise)

Differentiation

finite_difference()

df = finite_difference(f, h=1e-5, method="central")
result = df(x)  # -> float

Parameters:

• f: ScalarFunction - Function to differentiate
• h: float - Step size (default: 1e-5)
• method: str - One of "central", "forward", "backward" (default: "central")

Returns:

• ScalarFunction - Function representing numerical derivative

Convergence Rates:

• Central: O(h²)
• Forward/Backward: O(h)

ODE Solving

ODESolver (Abstract Base Class)

solver = Euler(dydt, t0, y0, tf, n)
t, y = solver.solve()  # -> tuple[List[float], List[float]]

Parameters:

• dydt: Callable[[float, float], float] - dy/dt = f(t, y)
• t0: float - Initial time
• y0: float - Initial value
• tf: float - Final time
• n: int - Number of steps

Returns:

• Tuple[List[float], List[float]] - Time points and solution values

Type Aliases

from numerical.types import ScalarFunction, Tolerance, MaxIterations

# ScalarFunction: Callable[[float], float]
# Tolerance: float
# MaxIterations: int

Exceptions

from numerical.exceptions import (
    NumericalError,      # Base exception
    ConvergenceError,    # Method did not converge
    BracketingError,     # Could not bracket root
    ValidationError      # Invalid parameters
)

Examples

Finding Multiple Roots

import numpy as np
from numerical.roots import Newton

def f(x):
    return x**3 - 5*x**2 + 8*x - 4

roots = []
for initial_guess in [0.5, 1.5, 2.5, 3.5]:
    solver = Newton(f, a=initial_guess, b=initial_guess+1)
    try:
        result = solver.solve()
        roots.append(result.root)
    except ConvergenceError:
        pass

print(f"Found roots: {roots}")

Numerical vs Analytical Integration

import math
from numerical.integration import Simpson

# Analytical: integral of x^2 from 0 to 1 = 1/3
numerical_result = Simpson(lambda x: x**2, 0, 1, n=100).compute()
analytical_result = 1/3

print(f"Numerical: {numerical_result:.10f}")
print(f"Analytical: {analytical_result:.10f}")
print(f"Error: {abs(numerical_result - analytical_result):.2e}")

Derivative Approximation

from numerical.differentiation import finite_difference

def f(x):
    return math.sin(x)

# Analytical: cos(x)
df = finite_difference(f)

x = math.pi / 4
numerical_derivative = df(x)
analytical_derivative = math.cos(x)

print(f"Numerical: {numerical_derivative:.10f}")
print(f"Analytical: {analytical_derivative:.10f}")

Project Structure

numerical/
├── __init__.py           # Main API exports
├── types.py              # Type aliases (ScalarFunction, etc)
├── exceptions.py         # Custom exceptions
├── result.py             # Result dataclasses
│
├── roots/
│   ├── __init__.py
│   ├── base.py          # RootFinder abstract base
│   ├── bisection.py
│   ├── false_position.py
│   ├── newton.py
│   └── utils.py
│
├── integration/
│   ├── __init__.py
│   ├── base.py          # Integrator abstract base
│   ├── riemann.py
│   ├── trapezoid.py
│   └── simpson.py
│
├── differentiation/
│   ├── __init__.py
│   └── finite_difference.py
│
└── ode/
    ├── __init__.py
    ├── base.py          # ODESolver abstract base
    ├── euler.py
    └── rk4.py

Performance Considerations

• Root Finding: O(n) function evaluations per iteration, where n depends on method
• Integration: Accuracy improves with more subintervals; CPU time grows linearly
• Differentiation: Central differences provide O(h²) accuracy for smooth functions
• ODE Solving: Local truncation error depends on method order and step size

Limitations and Known Issues

• All solvers assume scalar functions (single variable/ODE)
• Vector/system solvers not yet implemented
• Adaptive step size methods not yet implemented
• No automatic error estimation in ODE solvers

Contributing

Contributions welcome! Areas of interest:

• Vector-valued functions support
• Adaptive step size methods
• Error estimation
• Additional ODE solvers (RK45, BDF, etc)
• Performance optimizations
• Comprehensive test suite

License

MIT License - see LICENSE file for details

Author

Miguel Mochizuki Silva

• GitHub: @MiguelMochizukiDev
• Email: miguelmochizukisilva@gmail.com

Citation

If you use this library in academic work, please cite:

@software{numerical2026,
  author = {Silva, Miguel Mochizuki},
  title = {Numerical: Python Scientific Computing Library},
  year = {2026},
  url = {https://github.com/MiguelMochizukiDev/numerical}
}

Changelog

v1.0.0 (2026-05-07)

• Initial release
• Root finding (Bisection, False Position, Newton-Raphson)
• Numerical integration (Riemann, Trapezoid, Simpson)
• Numerical differentiation (Finite differences)
• ODE solving (Euler, Runge-Kutta 4th order)
• Complete type hints and documentation

References

• Burden & Faires - Numerical Analysis (10th ed.)