mathematicskit.calculus#

mathematicskit.calculus hub · History · Examples · API reference

Finite-difference derivatives with Richardson extrapolation (hand-rolled); numerical quadrature built on scipy.integrate (trapezoidal, Simpson’s, Gauss-Legendre via fixed_quad, adaptive via quad); forward-mode automatic differentiation via dual numbers and a reverse-mode (backpropagation-style) autodiff engine (hand-rolled – the pedagogical subject); Taylor/Maclaurin series expansion and convergence-radius estimation.

mathematicskit.calculus: numerical differentiation and integration.

Archimedes’ polygon bounds on pi; finite-difference derivatives (forward/backward/central, complex step) with Richardson extrapolation for higher accuracy (hand-rolled – no scipy equivalent); numerical quadrature built on scipy.integrate (composite trapezoidal, composite Simpson’s, Gauss-Legendre via fixed_quad, adaptive quadrature via quad) alongside hand-rolled Riemann sums, Romberg, Clenshaw-Curtis, tanh-sinh, and Euler-Maclaurin-corrected rules; forward-mode automatic differentiation via dual numbers and a small reverse-mode (backpropagation-style) autodiff engine (hand-rolled – autodiff is the pedagogical subject); Taylor/ Maclaurin series expansion and convergence-radius estimation for standard functions.

class mathematicskit.calculus.AdaptiveQuadrature(tol=1e-08, max_depth=50)[source]#

Bases: Quadrature

Adaptive quadrature with error estimation, via scipy.integrate.quad().

quad wraps QUADPACK’s QAGS: adaptive subdivision plus Gauss-Kronrod extrapolation, refining subintervals where the integrand is hardest to approximate until the requested absolute tolerance is met. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.6 (“Adaptive Quadrature Methods”), for the algorithm this generalizes.

Parameters:

  • tol (float) – Target absolute error (passed as epsabs).

  • max_depth (int) – Upper bound on subintervals (passed as limit).

Examples

>>> result = AdaptiveQuadrature(tol=1e-10).integrate(lambda x: x**2, 0.0, 1.0)
>>> round(result.value, 8)
0.33333333
integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.ClenshawCurtisQuadrature(n=16)[source]#

Bases: Quadrature

Clenshaw-Curtis quadrature: integrate the Chebyshev interpolant at \(\cos(k\pi/n)\).

For smooth integrands it converges nearly as fast as Gauss-Legendre quadrature with the same number of points, and its nodes nest when n doubles. See C. W. Clenshaw and A. R. Curtis, “A Method for Numerical Integration on an Automatic Computer,” Numerische Mathematik 2 (1960), 197-205.

Parameters:

n (int) – Number of intervals (n + 1 function evaluations).

Examples

>>> result = ClenshawCurtisQuadrature(n=16).integrate(np.exp, 0.0, 1.0)
>>> abs(result.value - (np.e - 1)) < 1e-14
True
integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.DerivativeResult(value, step=0.0, method='', error_estimate=0.0)[source]#

Bases: object

Container for a finite-difference derivative estimate.

Parameters:
error_estimate: float = 0.0#

Estimated error, when available (e.g. from Richardson extrapolation’s own convergence).

Type:

float

method: str = ''#

e.g. "forward", "backward", "central", "richardson".

Type:

str

step: float = 0.0#

Step size h used (the finest, for Richardson extrapolation).

Type:

float

value: float#

Estimated derivative.

Type:

float

class mathematicskit.calculus.Dual(real, dual=0.0)[source]#

Bases: object

A dual number \(\text{real} + \text{dual}\,\varepsilon\).

Parameters:

  • real (float) – The value part.

  • dual (float) – The derivative (tangent) part.

Examples

>>> x = Dual(3.0, 1.0)  # seed dx/dx = 1
>>> y = x * x + 2.0 * x  # f(x) = x^2 + 2x, f'(x) = 2x + 2
>>> y.real, y.dual
(15.0, 8.0)
cos()[source]#
dual#
exp()[source]#
log()[source]#
real#
sin()[source]#
sqrt()[source]#
class mathematicskit.calculus.ExhaustionResult(sides, lower, upper)[source]#

Bases: object

Container for Archimedes’ polygon bounds on \(\pi\).

Parameters:
lower: list#

Inscribed-polygon lower bounds on \(\pi\).

Type:

list of float

sides: list#

Number of polygon sides at each doubling (6, 12, 24, …).

Type:

list of int

upper: list#

Circumscribed-polygon upper bounds on \(\pi\).

Type:

list of float

class mathematicskit.calculus.GaussianQuadrature(n=5)[source]#

Bases: Quadrature

Gauss-Legendre quadrature: exact for polynomials up to degree 2n - 1.

Thin wrapper around scipy.integrate.fixed_quad(), which maps \([-1, 1]\) Gauss-Legendre nodes/weights onto [a, b]. For smooth integrands this needs far fewer evaluations than trapezoidal/Simpson for the same accuracy. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.7.

Parameters:

n (int) – Number of Gauss-Legendre nodes.

Examples

>>> result = GaussianQuadrature(n=5).integrate(lambda x: x**9, -1.0, 1.0)
>>> abs(round(result.value, 10))
0.0
integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.Quadrature[source]#

Bases: ABC

Common base for numerical-integration rules over [a, b].

abstractmethod integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.QuadratureResult(value, error_estimate=0.0, n_evaluations=0, method='', extra=<factory>)[source]#

Bases: object

Container for the output of a Quadrature rule.

Parameters:
error_estimate: float = 0.0#

Estimated absolute error (0 when the rule has no built-in estimate).

Type:

float

extra: dict#

Free-form diagnostics slot (e.g. adaptive quadrature’s subinterval history).

Type:

dict

method: str = ''#

e.g. "trapezoidal", "simpson", "gaussian", "adaptive".

Type:

str

n_evaluations: int = 0#

Number of function evaluations used.

Type:

int

value: float#

Estimated integral.

Type:

float

class mathematicskit.calculus.RiemannSum(n=100, rule='midpoint')[source]#

Bases: Quadrature

A Riemann sum on n equal subintervals, sampling each at its left end, right end, or midpoint.

Bernhard Riemann’s 1854 definition of the integral is the limit of such sums as the subintervals shrink. The left and right rules are \(O(h)\) accurate; the midpoint rule is \(O(h^2)\). Hand-rolled, since the sum is the definition being illustrated.

Parameters:

  • n (int) – Number of subintervals.

  • rule (str)

Examples

>>> RiemannSum(n=4, rule="left").integrate(lambda x: x, 0.0, 1.0).value
0.375
>>> RiemannSum(n=4, rule="midpoint").integrate(lambda x: x, 0.0, 1.0).value
0.5
integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.RombergQuadrature(levels=6)[source]#

Bases: Quadrature

Romberg integration: Richardson extrapolation applied to repeatedly halved trapezoidal rules.

Builds the triangular table

\[R_{i,0} = T_{2^i}, \qquad R_{i,j} = R_{i,j-1} + \frac{R_{i,j-1} - R_{i-1,j-1}}{4^j - 1},\]

where \(T_{2^i}\) is the trapezoidal rule on \(2^i\) subintervals. Each column cancels the next even power of \(h\) in the trapezoidal error expansion. Hand-rolled because scipy.integrate.romberg was removed in SciPy 1.15. See W. Romberg, “Vereinfachte numerische Integration,” Det Kongelige Norske Videnskabers Selskabs Forhandlinger 28(7) (1955), 30-36.

Parameters:

levels (int) – Number of halvings; the finest rule uses \(2^{\text{levels}}\) subintervals.

Examples

>>> result = RombergQuadrature(levels=5).integrate(np.exp, 0.0, 1.0)
>>> abs(result.value - (np.e - 1)) < 1e-12
True
integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.SimpsonsRule(n=100)[source]#

Bases: Quadrature

Composite Simpson’s rule, \(O(h^4)\) global error.

Fits a parabola through each pair of subintervals, via scipy.integrate.simpson() on n + 1 equally spaced samples (requires n even). See Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.4.

Parameters:

n (int) – Number of subintervals (must be even).

Examples

>>> result = SimpsonsRule(n=10).integrate(lambda x: x**3, 0.0, 1.0)
>>> round(result.value, 10)
0.25
integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.TanhSinhQuadrature(h=0.1, t_max=4.0)[source]#

Bases: Quadrature

Tanh-sinh (double-exponential) quadrature, robust to endpoint singularities.

Substitutes \(x = \tfrac{a+b}{2} + \tfrac{b-a}{2}\tanh(\tfrac{\pi}{2}\sinh t)\) and applies the trapezoidal rule in \(t\). The transformed integrand decays double-exponentially, so the trapezoidal rule converges very fast even when \(f\) blows up at an endpoint. The endpoints themselves are never evaluated. Nodes crowd toward the endpoints, so accuracy is limited by how precisely f can be evaluated there: an integrand such as \(1/\sqrt{1-x^2}\) loses digits once \(x\) is rounded near \(\pm1\). Hand-rolled for SciPy versions before 1.15, which lack scipy.integrate.tanhsinh. See H. Takahasi and M. Mori, “Double Exponential Formulas for Numerical Integration,” Publications of the Research Institute for Mathematical Sciences 9(3) (1974), 721-741.

Parameters:

  • h (float) – Step size in the transformed variable \(t\).

  • t_max (float) – The rule sums over \(t \in [-t_{\max}, t_{\max}]\).

Examples

>>> result = TanhSinhQuadrature(h=0.1).integrate(lambda x: 1 / np.sqrt(x), 0.0, 1.0)
>>> abs(result.value - 2.0) < 1e-10
True
integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.TrapezoidalRule(n=100)[source]#

Bases: Quadrature

Composite trapezoidal rule, \(O(h^2)\) global error.

\(\int_a^b f\,dx \approx h\left[\tfrac12 f(x_0) + f(x_1) + \dots + f(x_{n-1}) + \tfrac12 f(x_n)\right]\), \(h = (b-a)/n\), via scipy.integrate.trapezoid() on n + 1 equally spaced samples. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.3.

Parameters:

n (int) – Number of subintervals.

Examples

>>> result = TrapezoidalRule(n=1000).integrate(lambda x: x**2, 0.0, 1.0)
>>> round(result.value, 6)
0.333334
integrate(f, a, b)[source]#

Estimate \(\int_a^b f(x)\,dx\).

Parameters:
Return type:

QuadratureResult

Returns:

QuadratureResult

class mathematicskit.calculus.Variable(value, _parents=None)[source]#

Bases: object

A scalar node in a reverse-mode autodiff computation graph.

Parameters:

  • value (float) – The node’s forward value.

  • _parents (list[tuple[Variable, float]]) – Internal: (parent, local_partial_derivative) pairs recorded by the operation that created this node.

Examples

>>> x = Variable(3.0)
>>> y = x * x + 2.0 * x  # f(x) = x^2 + 2x
>>> y.backward()
>>> round(x.grad, 10)  # f'(3) = 2*3 + 2 = 8
8.0
backward()[source]#

Propagate gradients back through the computation graph.

Sets self.grad = 1 (the seed) and accumulates \(\partial \text{self}/\partial v\) into v.grad for every node v reachable as a parent, via reverse topological traversal.

cos()[source]#
exp()[source]#
grad#
log()[source]#
sin()[source]#
value#
mathematicskit.calculus.archimedes_pi_bounds(n_doublings=4)[source]#

Lower and upper bounds on \(\pi\) from regular polygons with \(6 \cdot 2^k\) sides.

Four doublings reach Archimedes’ 96-gon, which gives \(3\tfrac{10}{71} < \pi < 3\tfrac{1}{7}\).

Parameters:

n_doublings (int) – Number of side doublings after the initial hexagon.

Return type:

ExhaustionResult

Returns:

ExhaustionResult

Examples

>>> result = archimedes_pi_bounds(4)
>>> result.sides[-1]
96
>>> 3 + 10 / 71 < result.lower[-1] < math.pi < result.upper[-1] < 3 + 1 / 7
True
mathematicskit.calculus.backward_difference(f, x, h=1e-05)[source]#

Backward-difference derivative estimate, \(O(h)\).

\(f'(x) \approx \dfrac{f(x) - f(x-h)}{h}\).

Parameters:
Return type:

float

Returns:

float

Examples

>>> round(backward_difference(lambda t: t**2, 3.0, h=1e-4), 3)
6.0
mathematicskit.calculus.central_difference(f, x, h=1e-05)[source]#

Central-difference derivative estimate, \(O(h^2)\).

\(f'(x) \approx \dfrac{f(x+h) - f(x-h)}{2h}\) – one order of accuracy better than forward/backward differences at the same cost (2 evaluations), since the \(O(h)\) error terms in the underlying Taylor expansions cancel. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.1, eq. (4.4).

Parameters:
Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> round(float(central_difference(np.sin, 0.0, h=1e-4)), 8)
1.0
mathematicskit.calculus.clenshaw_curtis_nodes_and_weights(n)[source]#

Clenshaw-Curtis nodes \(x_k = \cos(k\pi/n)\) and weights on \([-1, 1]\).

The weights integrate exactly the Chebyshev interpolant through the \(n+1\) nodes. Computed with the explicit cosine-sum formula of L. N. Trefethen, Spectral Methods in MATLAB (Philadelphia: SIAM, 2000), program clencurt. Hand-rolled: SciPy has no Clenshaw-Curtis rule.

Parameters:

n (int) – Number of intervals; there are n + 1 nodes.

Returns:

nodes, weights (ndarray, shape (n + 1,))

Examples

>>> x, w = clenshaw_curtis_nodes_and_weights(4)
>>> round(float(w.sum()), 12)  # integrates 1 exactly: length of [-1, 1]
2.0
mathematicskit.calculus.complex_step_derivative(f, x, h=1e-20)[source]#

Complex-step derivative estimate, \(O(h^2)\) with no subtractive cancellation.

\(f'(x) \approx \operatorname{Im} f(x + ih) / h\). Because no two nearly equal numbers are subtracted, h can be taken absurdly small and the result is accurate to machine precision. Requires f to be real-analytic and implemented with complex-capable operations (e.g. numpy ufuncs). See J. N. Lyness and C. B. Moler, “Numerical Differentiation of Analytic Functions,” SIAM Journal on Numerical Analysis 4(2) (1967), 202-210, and W. Squire and G. Trapp, “Using Complex Variables to Estimate Derivatives of Real Functions,” SIAM Review 40(1) (1998), 110-112.

Parameters:

  • f (Callable) – Accepts a complex argument.

  • x (float)

  • h (float) – Imaginary step size.

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> abs(complex_step_derivative(np.exp, 1.0) - np.e) < 1e-15
True
mathematicskit.calculus.derivative(f, x)[source]#

Exact derivative of f at x via forward-mode dual numbers.

Parameters:

  • f (Callable[[Dual], Dual]) – A function built from Dual arithmetic and methods (+, -, *, /, **, .sin(), .cos(), .exp(), .log(), .sqrt()), taking and returning Dual.

  • x (float) – Point at which to differentiate.

Return type:

float

Returns:

float – \(f'(x)\), exact up to floating-point roundoff (no truncation error, unlike finite differences).

Examples

>>> f = lambda x: (x * x).sin()
>>> round(derivative(f, 1.0), 10) == round(2.0 * 1.0 * math.cos(1.0 * 1.0), 10)
True
mathematicskit.calculus.estimate_radius_of_convergence(coefficients, growth_factor=1.5)[source]#

Estimate a power series’ radius of convergence via the ratio test.

\(R = \lim_{n\to\infty}|a_n/a_{n+1}|\), estimated from the spacing-adjusted ratio \(|a_i/a_j|^{1/(j-i)}\) over consecutive nonzero coefficients (the \(1/(j-i)\) exponent is what makes series like "sin"/"cos", whose even or odd coefficients vanish identically, come out on the same footing as dense ones).

A truncated series can only ever sample that limit, so a series with infinite radius – "exp", "sin", "cos" – shows up not as a large number but as a ratio that keeps climbing with every term added. This function detects exactly that: if the estimate from the series’ last two nonzero coefficients exceeds the estimate from its first two by more than growth_factor, the ratio has not settled and inf is returned, rather than a number that is really just an artifact of where the series happened to be cut off (a degree-20 "exp" would otherwise report R = 20, and a degree-30 one R = 30). See any standard calculus text’s treatment of the ratio test for power series (e.g. Stewart, Calculus, 8th ed., Ch. 11.8).

Parameters:

  • coefficients (ndarray) – a[n] is the coefficient of \(x^n\) (as from maclaurin_coefficients()).

  • growth_factor (float) – How much larger the last ratio estimate must be than the first before the radius is judged unbounded. > 1.

Return type:

float

Returns:

float – The estimated radius, or inf when fewer than two nonzero coefficients are available to compare, or when the ratio is still growing (an unbounded radius of convergence).

Examples

>>> coeffs = maclaurin_coefficients("geometric", 30)
>>> round(estimate_radius_of_convergence(coeffs), 6)
1.0
>>> round(estimate_radius_of_convergence(maclaurin_coefficients("log1p", 30)), 2)
1.03
>>> # exp and sin converge everywhere: their ratio never settles.
>>> estimate_radius_of_convergence(maclaurin_coefficients("exp", 20))
inf
>>> estimate_radius_of_convergence(maclaurin_coefficients("sin", 20))
inf
mathematicskit.calculus.euler_maclaurin_trapezoid(f, a, b, n, odd_derivatives)[source]#

Trapezoidal rule plus Euler-Maclaurin endpoint corrections.

\[\int_a^b f\,dx \approx T_n - \sum_{k=1}^{m} \frac{B_{2k} h^{2k}}{(2k)!} \left(f^{(2k-1)}(b) - f^{(2k-1)}(a)\right),\]

with Bernoulli numbers \(B_{2k}\) from scipy.special.bernoulli(). Each correction removes one more even power of \(h\) from the trapezoidal error. Leonhard Euler (1735) and Colin Maclaurin (1742) found the formula independently. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.5.

Parameters:

  • f (Callable[[float], float])

  • a (float)

  • b (float)

  • n (int) – Number of trapezoidal subintervals.

  • odd_derivatives (list) – The odd-order derivatives of f, first, third, fifth and so on; one correction term is applied per derivative given.

Return type:

QuadratureResult

Returns:

QuadratureResult

Examples

>>> exact = np.e - 1
>>> plain = euler_maclaurin_trapezoid(np.exp, 0.0, 1.0, 8, [])
>>> corrected = euler_maclaurin_trapezoid(np.exp, 0.0, 1.0, 8, [np.exp, np.exp])
>>> abs(corrected.value - exact) < 1e-3 * abs(plain.value - exact)
True
mathematicskit.calculus.evaluate_series(coefficients, x)[source]#

Evaluate a truncated power series \(\sum_n a_n x^n\) (Horner’s method).

Parameters:
Returns:

float or ndarray

Examples

>>> import numpy as np
>>> coeffs = maclaurin_coefficients("exp", 15)
>>> round(float(evaluate_series(coeffs, 1.0)), 10) == round(np.e, 10)
True
mathematicskit.calculus.forward_difference(f, x, h=1e-05)[source]#

Forward-difference derivative estimate, \(O(h)\).

\(f'(x) \approx \dfrac{f(x+h) - f(x)}{h}\). See Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.1.

Parameters:
Return type:

float

Returns:

float

Examples

>>> round(forward_difference(lambda t: t**2, 3.0, h=1e-4), 3)
6.0
mathematicskit.calculus.gradient(f, x)[source]#

Gradient of a scalar function of several variables via reverse-mode AD.

Parameters:
Return type:

list[float]

Returns:

list of float – \(\nabla f(x)\).

Examples

>>> f = lambda x, y: x * x * y + y  # df/dx = 2xy, df/dy = x^2 + 1
>>> grad = gradient(f, [3.0, 2.0])
>>> [round(g, 10) for g in grad]
[12.0, 10.0]
mathematicskit.calculus.legendre_nodes_and_weights(n)[source]#

Gauss-Legendre nodes and weights on \([-1, 1]\), via numpy.polynomial.legendre.leggauss().

Nodes are the n roots of the degree-n Legendre polynomial \(P_n\); weights are \(w_i = \dfrac{2}{(1-x_i^2)[P_n'(x_i)]^2}\). leggauss computes both from the companion-matrix eigenvalues of \(P_n\), which is more robust than a from-scratch Newton iteration on the asymptotic initial guess. See Press et al., Numerical Recipes, 3rd ed., Sec. 4.6, and Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.7.

Parameters:

n (int) – Number of nodes (exact for polynomials up to degree 2n - 1).

Returns:

nodes, weights (ndarray, shape (n,))

Examples

>>> nodes, weights = legendre_nodes_and_weights(3)
>>> import numpy as np
>>> round(float(np.sum(weights)), 10)
2.0
mathematicskit.calculus.maclaurin_coefficients(name, order)[source]#

Maclaurin-series coefficients \(a_0, \dots, a_{\text{order}}\) for a standard function.

Supported name values and their series:

  • "exp": \(e^x = \sum_{n\geq0} x^n/n!\), \(R=\infty\).

  • "sin": \(\sin x = \sum_{n\geq0} (-1)^n x^{2n+1}/(2n+1)!\), \(R=\infty\).

  • "cos": \(\cos x = \sum_{n\geq0} (-1)^n x^{2n}/(2n)!\), \(R=\infty\).

  • "log1p": \(\ln(1+x) = \sum_{n\geq1} (-1)^{n+1} x^n/n\), \(R=1\).

  • "geometric": \(1/(1-x) = \sum_{n\geq0} x^n\), \(R=1\).

  • "arctan": \(\arctan x = \sum_{n\geq0} (-1)^n x^{2n+1}/(2n+1)\), \(R=1\).

Parameters:

  • name (str) – One of the supported function names above.

  • order (int) – Highest power to include (order + 1 coefficients returned).

Return type:

ndarray

Returns:

ndarray, shape (order + 1,) – a[n] is the coefficient of \(x^n\).

Examples

>>> import numpy as np
>>> coeffs = maclaurin_coefficients("exp", 4)
>>> np.allclose(coeffs, [1.0, 1.0, 0.5, 1.0 / 6.0, 1.0 / 24.0])
True
mathematicskit.calculus.partial_sums(coefficients, x)[source]#

Partial sums \(S_k(x) = \sum_{n=0}^{k} a_n x^n\) for every k.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (order + 1,) – result[k] is the degree-k partial sum.

Examples

>>> import numpy as np
>>> sums = partial_sums(np.array([1.0, 1.0, 0.5, 1.0 / 6.0]), 1.0)
>>> np.round(sums, 6)
array([1.      , 2.      , 2.5     , 2.666667])
mathematicskit.calculus.richardson_extrapolation(f, x, h=0.01, levels=4)[source]#

Richardson-extrapolate central differences to arbitrarily high order.

Central differences have an error series in even powers of h: \(D(h) = f'(x) + c_1 h^2 + c_2 h^4 + \dots\). Combining \(D(h)\) and \(D(h/2)\) as \(\dfrac{4 D(h/2) - D(h)}{3}\) cancels the \(h^2\) term, leaving \(O(h^4)\); repeating with successively halved step sizes in a triangular (Neville-style) table cancels one more error order each level. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 4.2, and the closely related idea in Ch. 4.5 (Romberg integration).

Parameters:

  • f (Callable[[float], float])

  • x (float)

  • h (float) – Coarsest step size; halved at each level.

  • levels (int) – Number of extrapolation levels (levels=1 is plain central difference at step h).

Return type:

DerivativeResult

Returns:

DerivativeResult – error_estimate is the change between the last two extrapolated values, a practical (not rigorous) error indicator.

Examples

>>> import numpy as np
>>> result = richardson_extrapolation(np.sin, 0.0, h=0.2, levels=5)
>>> round(float(result.value), 12)
1.0
mathematicskit.calculus.taylor_remainder_bound(max_derivative_bound, order, x, x0=0.0)[source]#

Lagrange remainder bound for a degree-order Taylor polynomial.

\(|R_n(x)| \leq \dfrac{M}{(n+1)!}|x-x_0|^{n+1}\), where M bounds \(|f^{(n+1)}|\) on the interval between x0 and x. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 1.1, Theorem 1.14 (Taylor’s Theorem).

Parameters:

  • max_derivative_bound (float) – A bound M on \(|f^{(\text{order}+1)}|\) on the relevant interval.

  • order (int) – Degree of the Taylor polynomial.

  • x (float) – Evaluation point and expansion center.

  • x0 (float) – Evaluation point and expansion center.

Return type:

float

Returns:

float

Examples

>>> # sin has |f^(n+1)| <= 1 everywhere, so this bounds the 5th-order
>>> # Maclaurin remainder at x=0.5.
>>> round(taylor_remainder_bound(1.0, 5, 0.5), 8)
2.17e-05
mathematicskit.calculus.truncation_error(f, coefficients, x)[source]#

Empirical truncation error |f(x) - S_k(x)| for every partial sum.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (order + 1,)

Examples

>>> import numpy as np
>>> from mathematicskit.calculus.systems.taylor_series import maclaurin_coefficients
>>> coeffs = maclaurin_coefficients("exp", 15)
>>> errors = truncation_error(np.exp, coeffs, 1.0)
>>> bool(errors[-1] < errors[0])
True