mathematicskit.numerical_analysis#

mathematicskit.numerical_analysis hub · History · Examples · API reference

Root finding (bisection, Newton-Raphson, secant, fixed-point iteration; hand-rolled for their per-iterate convergence history, cross-checked against scipy.optimize) with convergence-order verification; Lagrange and Newton divided-difference polynomial interpolation (hand-rolled, cross-checked against scipy.interpolate); cubic spline interpolation via scipy.interpolate.CubicSpline; Chebyshev interpolation nodes (numpy.polynomial.chebyshev.chebpts2) and the Runge phenomenon; least-squares polynomial regression via numpy.linalg.lstsq; and shared error/stability-analysis utilities (numpy.linalg.cond-based condition number, the hand-rolled Lebesgue constant and empirical convergence order).

mathematicskit.numerical_analysis: root finding and polynomial approximation.

Root finding (bisection, Newton-Raphson, secant, fixed-point iteration, hand-rolled to expose per-iterate convergence history, cross-checked against scipy.optimize in tests) with convergence-order verification; Lagrange and Newton divided-difference polynomial interpolation (hand-rolled, cross-checked against scipy.interpolate in tests); cubic spline interpolation (natural and clamped) built on scipy.interpolate.CubicSpline; Chebyshev interpolation nodes (numpy.polynomial.chebyshev.chebpts2) and the Runge phenomenon; least-squares polynomial regression via numpy.linalg.lstsq; Halley’s and Steffensen’s methods and Aitken’s delta-squared acceleration; Hermite interpolation (scipy.interpolate.KroghInterpolator); Bernstein polynomials, Padé approximants, and Remez minimax approximation; Horner’s scheme and Wilkinson’s polynomial with root condition numbers; Kahan’s compensated summation; and shared error/stability-analysis utilities (empirical convergence order, the Lebesgue constant as an interpolation problem’s condition number, and a general condition-number wrapper around numpy.linalg.cond).

class mathematicskit.numerical_analysis.Bisection(f, a, b, tol=1e-08, max_iter=1000)[source]#

Bases: IterativeRootFinder

Bisection method on a bracket [a, b] with f(a) * f(b) < 0.

Halves the bracket each iteration, keeping the half where the sign change persists; converges linearly (order 1) with error bound \(|x_n - \text{root}| \leq (b - a) / 2^{n+1}\). See Burden & Faires, Numerical Analysis, 10th ed., Ch. 2.1.

Parameters:

f (Callable[[float], float]) – Continuous function f(x) -> float.
a (float) – Bracket endpoints with f(a) and f(b) of opposite sign.
b (float) – Bracket endpoints with f(a) and f(b) of opposite sign.
tol (float) – Stop once the bracket width is below tol.
max_iter (int) – Maximum number of bisections.

Examples

>>> result = Bisection(lambda x: x**2 - 2.0, 0.0, 2.0, tol=1e-10).solve()
>>> round(result.root, 6)
1.414214
>>> result.converged
True
>>> # ``history`` holds exactly one midpoint per iteration performed.
>>> len(result.history) == result.iterations
True

solve()[source]#

Run the iteration to convergence (or max_iter).

Return type:: RootResult
Returns:: RootResult

class mathematicskit.numerical_analysis.ChebyshevInterpolant(f, n, a=-1.0, b=1.0)[source]#

Bases: Interpolant

Polynomial interpolant of a function sampled at Chebyshev nodes.

Evaluated with the barycentric interpolation formula specialized to Chebyshev-of-the-second-kind nodes (weights \(w_k = (-1)^k \delta_k\), \(\delta_k = 1/2\) at the endpoints and 1 otherwise), which is both faster (\(O(n)\) per point after setup) and far better conditioned than the raw Lagrange form. See Trefethen, Approximation Theory and Approximation Practice, 2013, Ch. 5, eq. (5.3)-(5.4).

Parameters:

f (callable) – Function to interpolate, f(x) -> float.
n (int) – Number of Chebyshev nodes.
a (float) – Interval endpoints.
b (float) – Interval endpoints.

Examples

>>> import numpy as np
>>> p = ChebyshevInterpolant(np.sin, n=20, a=-3.0, b=3.0)
>>> bool(abs(p.evaluate(1.0) - np.sin(1.0)) < 1e-10)
True

evaluate(x_new)[source]#

Evaluate the interpolant at new point(s).

Parameters:: x_new (float or array-like of float)
Returns:: float or ndarray

class mathematicskit.numerical_analysis.CubicSpline(x, y, boundary='natural', fpa=None, fpb=None)[source]#

Bases: Interpolant

Piecewise cubic spline through (x_i, y_i), natural or clamped.

On each subinterval \([x_i, x_{i+1}]\),

\[S_i(x) = a_i + b_i (x - x_i) + c_i (x - x_i)^2 + d_i (x - x_i)^3\]

with \(S, S', S''\) continuous across nodes. A natural spline fixes \(S''(x_0) = S''(x_n) = 0\); a clamped spline instead matches prescribed end-slopes \(S'(x_0) = f'_a\), \(S'(x_n) = f'_b\). Delegates to scipy.interpolate.CubicSpline with bc_type="natural" or bc_type=((1, fpa), (1, fpb)) respectively. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 3.5, Algorithms 3.4-3.5.

Parameters:

x (array-like, shape (n + 1,)) – Interpolation nodes (strictly increasing) and values.
y (array-like, shape (n + 1,)) – Interpolation nodes (strictly increasing) and values.
boundary (str) – Boundary-condition type.
fpa (float | None) – End-slopes \(S'(x_0)\), \(S'(x_n)\); required when boundary="clamped".
fpb (float | None) – End-slopes \(S'(x_0)\), \(S'(x_n)\); required when boundary="clamped".

Examples

>>> import numpy as np
>>> x = np.array([0.0, 1.0, 2.0, 3.0])
>>> y = np.array([0.0, 1.0, 0.0, 1.0])
>>> spline = CubicSpline(x, y, boundary="natural")
>>> # Interpolates exactly at the nodes themselves.
>>> np.allclose([spline.evaluate(xi) for xi in x], y)
True
>>> # A straight line is reproduced exactly by a clamped spline whose
>>> # end-slopes match the line's slope.
>>> line = CubicSpline(x, 2.0 * x + 1.0, boundary="clamped", fpa=2.0, fpb=2.0)
>>> round(float(line.evaluate(1.7)), 10)
4.4

evaluate(x_new)[source]#

Evaluate the interpolant at new point(s).

Parameters:: x_new (float or array-like of float)
Returns:: float or ndarray

second_derivative_at_nodes()[source]#

Return \(S''(x_i)\) at every node.

Return type:: ndarray
Returns:: ndarray, shape (n + 1,)

class mathematicskit.numerical_analysis.FixedPointIteration(g, x0, tol=1e-08, max_iter=1000)[source]#

Bases: IterativeRootFinder

Fixed-point iteration: \(x_{n+1} = g(x_n)\), converging to a fixed point x* = g(x*) (equivalently a root of f(x) = g(x) - x).

Converges linearly if |g'(x*)| < 1 in a neighborhood of the fixed point (the contraction condition); see Burden & Faires, Numerical Analysis, 10th ed., Ch. 2.2, Theorem 2.4 (Fixed-Point Theorem) and Corollary 2.5.

Parameters:

g (Callable[[float], float]) – Iteration function g(x) -> float.
x0 (float) – Initial guess.
tol (float) – Stop once |x_{n+1} - x_n| < tol.
max_iter (int) – Maximum number of iterations.

Examples

>>> # x* solves x = cos(x); g(x) = cos(x) is a contraction near x*.
>>> import numpy as np
>>> result = FixedPointIteration(np.cos, x0=0.5, tol=1e-10).solve()
>>> bool(abs(result.root - np.cos(result.root)) < 1e-9)
True

solve()[source]#

Run the iteration to convergence (or max_iter).

Return type:: RootResult
Returns:: RootResult

class mathematicskit.numerical_analysis.Halley(f, fprime, fprime2, x0, tol=1e-08, max_iter=1000)[source]#

Bases: IterativeRootFinder

Halley’s method: a cubically convergent refinement of Newton’s method.

\[x_{n+1} = x_n - \frac{2 f(x_n) f'(x_n)}{2 f'(x_n)^2 - f(x_n) f''(x_n)}\]

Equivalent to applying Newton’s method to \(f/\sqrt{|f'|}\), or to stepping to the root of the osculating hyperbola at \(x_n\). Converges with order 3 near a simple root. See E. Halley, “Methodus nova accurata & facilis inveniendi radices aequationum quarumcumque generaliter, sine praevia reductione,” Philosophical Transactions 18 (1694), 136-148; and T. R. Scavo and J. B. Thoo, “On the Geometry of Halley’s Method,” American Mathematical Monthly 102 (1995), 417-426. The same iteration is what scipy.optimize.newton() runs when given fprime2; it is hand-rolled here to expose the iterates.

Parameters:

f (Callable[[float], float]) – Function and its first and second derivatives, all float -> float.
fprime (Callable[[float], float]) – Function and its first and second derivatives, all float -> float.
fprime2 (Callable[[float], float]) – Function and its first and second derivatives, all float -> float.
x0 (float) – Initial guess.
tol (float) – Stop once |x_{n+1} - x_n| < tol.
max_iter (int) – Maximum number of iterations.

Examples

>>> result = Halley(lambda x: x**2 - 2.0, lambda x: 2.0 * x, lambda x: 2.0, x0=1.0).solve()
>>> round(result.root, 12)
1.414213562373
>>> result.iterations <= 4
True

solve()[source]#

Run the iteration to convergence (or max_iter).

Return type:: RootResult
Returns:: RootResult

class mathematicskit.numerical_analysis.HermiteInterpolant(x, y, dydx)[source]#

Bases: Interpolant

Hermite interpolating polynomial matching values and slopes.

The unique polynomial \(H\) of degree \(\le 2n + 1\) with \(H(x_i) = y_i\) and \(H'(x_i) = y'_i\) at \(n + 1\) distinct nodes. In Newton form it is the divided-difference interpolant on the doubled node list \(x_0, x_0, x_1, x_1, \ldots\), with each repeated first difference \(f[x_i, x_i]\) replaced by \(y'_i\). The error is

\[f(x) - H(x) = \frac{f^{(2n+2)}(\xi)}{(2n+2)!} \prod_{i=0}^{n} (x - x_i)^2 .\]

See C. Hermite, “Sur la formule d’interpolation de Lagrange,” Journal für die reine und angewandte Mathematik 84 (1878), 70-79; Burden & Faires, Numerical Analysis, 10th ed., Ch. 3.4. Built on scipy.interpolate.KroghInterpolator, which accepts repeated nodes as derivative conditions (F. T. Krogh, Mathematics of Computation 24 (1970), 185-190).

Parameters:

x (array-like, shape (n + 1,)) – Distinct interpolation nodes and values.
y (array-like, shape (n + 1,)) – Distinct interpolation nodes and values.
dydx (array-like, shape (n + 1,)) – First derivatives at the nodes.

Examples

>>> # Two nodes with slopes determine a cubic: recover x^3 exactly.
>>> H = HermiteInterpolant([0.0, 1.0], [0.0, 1.0], dydx=[0.0, 3.0])
>>> round(H(0.5), 12)
0.125
>>> round(float(H.derivative(1.0)), 12)
3.0

derivative(x_new)[source]#: Evaluate \(H'\) at x_new.

evaluate(x_new)[source]#

Evaluate the interpolant at new point(s).

Parameters:: x_new (float or array-like of float)
Returns:: float or ndarray

class mathematicskit.numerical_analysis.HornerResult(value, derivative, quotient)[source]#

Bases: object

Container for one step of Horner’s scheme: \(p(x_0)\) and the deflated quotient \(q(x)\) with \(p(x) = (x - x_0)\,q(x) + p(x_0)\).

Parameters:

value (float)
derivative (float)
quotient (ndarray)

derivative: float#

\(p'(x_0) = q(x_0)\).

Type:: float

quotient: ndarray#

Coefficients of \(q\), highest power first (numpy.polyval convention).

Type:: ndarray, shape (degree,)

value: float#

\(p(x_0)\), the remainder of dividing by \(x - x_0\).

Type:: float

class mathematicskit.numerical_analysis.Interpolant(x, y)[source]#

Bases: ABC

Common base for algorithms that build an approximation to tabulated data (x_i, y_i) and can be evaluated at new points.

Parameters:

x (array-like, shape (n,)) – Interpolation nodes and values, respectively.
y (array-like, shape (n,)) – Interpolation nodes and values, respectively.

abstractmethod evaluate(x_new)[source]#

Evaluate the interpolant at new point(s).

Parameters:: x_new (float or array-like of float)
Returns:: float or ndarray

class mathematicskit.numerical_analysis.IterativeRootFinder(tol=1e-08, max_iter=1000)[source]#

Bases: ABC

Common base for scalar root-finding algorithms.

Subclasses store whatever state they need in __init__ (a function, a bracket or initial guess(es)) and implement solve().

Parameters:

tol (float) – Convergence tolerance (interpretation is method-specific: an absolute change in the iterate, or a bracket width).
max_iter (int) – Maximum number of iterations before giving up.

abstractmethod solve()[source]#

Run the iteration to convergence (or max_iter).

Return type:: RootResult
Returns:: RootResult

class mathematicskit.numerical_analysis.LagrangeInterpolant(x, y)[source]#

Bases: Interpolant

Lagrange-form interpolating polynomial through (x_i, y_i).

\[p(x) = \sum_{i=0}^{n} y_i L_i(x), \qquad L_i(x) = \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}\]

Evaluated directly from this basis-function sum (\(O(n^2)\) per evaluation point); see Burden & Faires, Numerical Analysis, 10th ed., Ch. 3.1, Theorem 3.2.

Parameters:

x (array-like, shape (n + 1,)) – Interpolation nodes (distinct) and values.
y (array-like, shape (n + 1,)) – Interpolation nodes (distinct) and values.

Examples

>>> import numpy as np
>>> # Exact for polynomials up to the interpolation degree.
>>> x = np.array([0.0, 1.0, 2.0, 3.0])
>>> y = x**3 - 2.0 * x + 1.0
>>> p = LagrangeInterpolant(x, y)
>>> round(float(p.evaluate(1.5)), 10)
1.375
>>> round(1.5**3 - 2.0 * 1.5 + 1.0, 10)
1.375

evaluate(x_new)[source]#

Evaluate the interpolant at new point(s).

Parameters:: x_new (float or array-like of float)
Returns:: float or ndarray

class mathematicskit.numerical_analysis.MinimaxResult(coefficients, max_error, reference, iterations, converged)[source]#

Bases: object

Container for the output of the Remez minimax approximation.

Parameters:

coefficients (ndarray)
max_error (float)
reference (ndarray)
iterations (int)
converged (bool)

coefficients: ndarray#

Best-approximation polynomial in the power basis, highest power first (numpy.polyval convention).

Type:: ndarray, shape (degree + 1,)

converged: bool#

Whether the error equioscillated to within the tolerance.

Type:: bool

evaluate(x)[source]#: Evaluate the minimax polynomial at x.

iterations: int#

Number of exchange steps performed.

Type:: int

max_error: float#

\(\max |f - p|\) over the interval, measured on a fine grid.

Type:: float

reference: ndarray#

Final reference (alternation) points.

Type:: ndarray, shape (degree + 2,)

class mathematicskit.numerical_analysis.NewtonDividedDifference(x, y)[source]#

Bases: Interpolant

Newton divided-difference form of the interpolating polynomial.

\[p(x) = f[x_0] + f[x_0,x_1](x-x_0) + f[x_0,x_1,x_2](x-x_0)(x-x_1) + \dots\]

built from the divided-difference table \(f[x_i,\dots,x_{i+k}] = \dfrac{f[x_{i+1},\dots,x_{i+k}] - f[x_i,\dots,x_{i+k-1}]}{x_{i+k}-x_i}\) and evaluated with nested (Horner-like) multiplication, \(O(n)\) per point after an \(O(n^2)\) table build – cheaper to extend with a new node than rebuilding a Lagrange form from scratch. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 3.3.

Parameters:

x (array-like, shape (n + 1,)) – Interpolation nodes (distinct) and values.
y (array-like, shape (n + 1,)) – Interpolation nodes (distinct) and values.

Examples

>>> import numpy as np
>>> x = np.array([0.0, 1.0, 2.0, 3.0])
>>> y = x**3 - 2.0 * x + 1.0
>>> p = NewtonDividedDifference(x, y)
>>> round(float(p.evaluate(1.5)), 10)
1.375

coefficients#

divided-difference coefficients f[x0], f[x0,x1], …

Type:: ndarray, shape (n,)

property divided_difference_table: ndarray#

full divided-difference table (upper triangle populated).

Type:: ndarray, shape (n, n)

evaluate(x_new)[source]#

Evaluate the interpolant at new point(s).

Parameters:: x_new (float or array-like of float)
Returns:: float or ndarray

class mathematicskit.numerical_analysis.NewtonRaphson(f, fprime, x0, tol=1e-08, max_iter=1000)[source]#

Bases: IterativeRootFinder

Newton-Raphson method: \(x_{n+1} = x_n - f(x_n) / f'(x_n)\).

Converges quadratically (order 2) near a simple root, provided \(f'\) doesn’t vanish there. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 2.3, Theorem 2.9.

Parameters:

f (Callable[[float], float]) – Function and its derivative, both float -> float.
fprime (Callable[[float], float]) – Function and its derivative, both float -> float.
x0 (float) – Initial guess.
tol (float) – Stop once |x_{n+1} - x_n| < tol.
max_iter (int) – Maximum number of iterations.

Examples

>>> result = NewtonRaphson(lambda x: x**2 - 2.0, lambda x: 2.0 * x, x0=1.0).solve()
>>> round(result.root, 10)
1.4142135624

solve()[source]#

Run the iteration to convergence (or max_iter).

Return type:: RootResult
Returns:: RootResult

class mathematicskit.numerical_analysis.PadeApproximant(taylor_coefficients, m, n)[source]#

Bases: object

Padé approximant \([m/n]\) of a power series: the rational function \(p(x)/q(x)\), \(\deg p \le m\), \(\deg q \le n\), \(q(0) = 1\), whose Taylor series agrees with \(f\) through \(x^{m+n}\).

With Taylor coefficients \(c_i\) (and \(c_i = 0\) for \(i < 0\)), the denominator solves the \(n \times n\) linear system \(\sum_{j=0}^{n} q_j c_{m+k-j} = 0\), \(k = 1, \ldots, n\), and then \(p_i = \sum_{j=0}^{\min(i, n)} q_j c_{i-j}\). See H. Padé, “Sur la représentation approchée d’une fonction par des fractions rationnelles,” Annales scientifiques de l’École Normale Supérieure (3) 9 (1892), supplement, 3-93; G. A. Baker and P. Graves-Morris, Padé Approximants, 2nd ed. (Cambridge University Press, 1996), Ch. 1.

Parameters:

taylor_coefficients (array-like) – \(c_0, c_1, \ldots\), lowest order first; at least m + n + 1 of them.
m (int) – Numerator and denominator degrees.
n (int) – Numerator and denominator degrees.

numerator, denominator

Coefficients of \(p\) and \(q\), lowest order first.

Type:: ndarray

Examples

>>> import math
>>> c = [1.0 / math.factorial(k) for k in range(3)]
>>> r = PadeApproximant(c, 1, 1)  # exp(x) ~ (1 + x/2) / (1 - x/2)
>>> r.numerator.tolist(), r.denominator.tolist()
([1.0, 0.5], [1.0, -0.5])
>>> round(r(1.0), 10)
3.0

evaluate(x)[source]#: Evaluate \(p(x)/q(x)\).

class mathematicskit.numerical_analysis.PolynomialRegression(x, y, degree)[source]#

Bases: object

Least-squares polynomial fit of degree degree to (x_i, y_i).

Parameters:

x (array-like, shape (n,)) – Data points.
y (array-like, shape (n,)) – Data points.
degree (int) – Polynomial degree (degree < n is required so the fit isn’t underdetermined into exact interpolation).

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> x = np.linspace(-1.0, 1.0, 50)
>>> y = 3.0 * x**2 - 2.0 * x + 1.0  # noiseless: fit should be essentially exact
>>> model = PolynomialRegression(x, y, degree=2)
>>> result = model.fit()
>>> np.allclose(result.coefficients, [3.0, -2.0, 1.0], atol=1e-8)
True
>>> round(result.r_squared, 6)
1.0

fit()[source]#

Fit the model and return the result.

Return type:: RegressionResult
Returns:: RegressionResult

predict(x_new)[source]#

Evaluate the fitted polynomial at new points.

Fits the model on first use and caches the coefficients, so calling fit() beforehand is optional; the cache is never invalidated, since a PolynomialRegression instance’s data and degree are fixed at construction.

Parameters:: x_new (array-like of float)
Return type:: ndarray
Returns:: ndarray

class mathematicskit.numerical_analysis.RegressionResult(coefficients, fitted_values, residuals, r_squared, adjusted_r_squared, condition_number=None)[source]#

Bases: object

Container for the output of a least-squares regression fit.

Parameters:

coefficients (ndarray)
fitted_values (ndarray)
residuals (ndarray)
r_squared (float)
adjusted_r_squared (float)
condition_number (float | None)

adjusted_r_squared: float#

\(R^2\) adjusted for the number of predictors.

Type:: float

coefficients: ndarray#

Fitted coefficients, highest power first (numpy.polyval convention).

Type:: ndarray, shape (degree + 1,)

condition_number: float | None = None#

Condition number of the design (Vandermonde) matrix used for the fit, a diagnostic for numerical stability at high polynomial degree (see mathematicskit.numerical_analysis.utils.error_analysis).

Type:: float, optional

fitted_values: ndarray#

Model predictions at the training points.

Type:: ndarray, shape (n,)

r_squared: float#

Coefficient of determination \(R^2\).

Type:: float

residuals: ndarray#

y - fitted_values.

Type:: ndarray, shape (n,)

class mathematicskit.numerical_analysis.RootResult(root, converged, iterations, history, method='', extra=<factory>)[source]#

Bases: object

Container for the output of an IterativeRootFinder.

Parameters:

root (float)
converged (bool)
iterations (int)
history (ndarray)
method (str)
extra (dict)

converged: bool#

Whether the stopping tolerance was met before max_iter.

Type:: bool

extra: dict#

Free-form slot for method-specific diagnostics (e.g. the final bracket for bisection).

Type:: dict

history: ndarray#

The sequence of iterates, starting from the initial guess(es), used for convergence-order analysis (see mathematicskit.numerical_analysis.utils.error_analysis).

Type:: ndarray, shape (iterations + 1,)

iterations: int#

Number of iterations actually performed.

Type:: int

method: str = ''#

Name of the method used (e.g. "bisection").

Type:: str

root: float#

The final estimate of the root.

Type:: float

class mathematicskit.numerical_analysis.Secant(f, x0, x1, tol=1e-08, max_iter=1000)[source]#

Bases: IterativeRootFinder

Secant method: Newton’s method with the derivative replaced by a finite-difference slope through the two most recent iterates.

\(x_{n+1} = x_n - f(x_n) \dfrac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}\). Converges superlinearly with order \((1+\sqrt5)/2 \approx 1.618\) (the golden ratio) near a simple root – see Burden & Faires, Numerical Analysis, 10th ed., Ch. 2.3, Theorem 2.10 – without requiring an analytic derivative.

Parameters:

f (Callable[[float], float]) – Function f(x) -> float.
x0 (float) – Two initial guesses.
x1 (float) – Two initial guesses.
tol (float) – Stop once |x_{n+1} - x_n| < tol.
max_iter (int) – Maximum number of iterations.

Examples

>>> result = Secant(lambda x: x**2 - 2.0, x0=1.0, x1=2.0).solve()
>>> round(result.root, 8)
1.41421356

solve()[source]#

Run the iteration to convergence (or max_iter).

Return type:: RootResult
Returns:: RootResult

class mathematicskit.numerical_analysis.Steffensen(g, x0, tol=1e-08, max_iter=1000)[source]#

Bases: IterativeRootFinder

Steffensen’s method: fixed-point iteration accelerated by Aitken’s \(\Delta^2\) process at every step.

From \(x_n\), compute \(g(x_n)\) and \(g(g(x_n))\), then jump to their Aitken extrapolation

\[x_{n+1} = x_n - \frac{\bigl(g(x_n) - x_n\bigr)^2}{g(g(x_n)) - 2 g(x_n) + x_n}.\]

Converges quadratically to a fixed point with \(g'(x^*) \neq 1\), even where plain iteration converges only linearly (or diverges), and without any derivative. See J. F. Steffensen, “Remarks on iteration,” Skandinavisk Aktuarietidskrift 16 (1933), 64-72; Burden & Faires, Numerical Analysis, 10th ed., Ch. 2.5, Algorithm 2.6. The same iteration is scipy.optimize.fixed_point() with method="del2"; it is hand-rolled here to expose the iterates.

Parameters:

g (Callable[[float], float]) – Iteration function g(x) -> float whose fixed point is sought.
x0 (float) – Initial guess.
tol (float) – Stop once |x_{n+1} - x_n| < tol.
max_iter (int) – Maximum number of iterations.

Examples

>>> import numpy as np
>>> result = Steffensen(np.cos, x0=0.5, tol=1e-12).solve()
>>> round(result.root, 10)
0.7390851332
>>> result.iterations < 6
True

solve()[source]#

Run the iteration to convergence (or max_iter).

Return type:: RootResult
Returns:: RootResult

mathematicskit.numerical_analysis.aitken_delta_squared(sequence)[source]#

Apply Aitken’s \(\Delta^2\) process to a sequence.

\[\hat a_n = a_n - \frac{(\Delta a_n)^2}{\Delta^2 a_n} = a_n - \frac{(a_{n+1} - a_n)^2}{a_{n+2} - 2 a_{n+1} + a_n}\]

The transform is exact for sequences of the form \(a_n = L + c\,r^n\) (\(r \neq 1\)), and for any linearly convergent sequence it converges to the limit faster than the original. Where \(\Delta^2 a_n = 0\), \(a_{n+2}\) is returned.

Parameters:: sequence (array-like, shape (n,)) – Terms \(a_0, \ldots, a_{n-1}\), with n >= 3.
Return type:: ndarray
Returns:: ndarray, shape (n - 2,) – The accelerated sequence \(\hat a_0, \ldots, \hat a_{n-3}\).

Examples

>>> import numpy as np
>>> # a_n = 2 + 3 * 0.5**n: one Aitken step recovers the limit exactly.
>>> a = 2.0 + 3.0 * 0.5 ** np.arange(6)
>>> np.allclose(aitken_delta_squared(a), 2.0)
True

mathematicskit.numerical_analysis.bernstein_polynomial(f, n, x)[source]#

Evaluate the degree-n Bernstein polynomial of f on \([0, 1]\).

\[(B_n f)(x) = \sum_{k=0}^{n} f\!\left(\tfrac{k}{n}\right) \binom{n}{k} x^k (1 - x)^{n-k}\]

The weights are the binomial probabilities \(P(K = k)\) for \(K \sim \mathrm{Bin}(n, x)\), so \(B_n f(x) = \mathbb{E}[f(K/n)]\), and the law of large numbers forces \(B_n f \to f\) uniformly for every continuous f: a constructive proof of the Weierstrass approximation theorem (S. Bernstein, Comm. Soc. Math. Kharkov (2) 13 (1912), 1-2). Convergence is slow: \(B_n(x^2) = x^2 + x(1 - x)/n\).

Parameters:

f (Callable) – Vectorized function on \([0, 1]\).
n (int) – Degree, n >= 1.
x (float or array-like) – Evaluation point(s) in \([0, 1]\).

Returns:

float or ndarray

Examples

>>> round(bernstein_polynomial(lambda t: t**2, 4, 0.5), 10)  # 0.25 + 0.25/4
0.3125

mathematicskit.numerical_analysis.chebyshev_nodes(n, a=-1.0, b=1.0)[source]#

Return the n Chebyshev points of the second kind on [a, b].

\[x_k = \frac{a+b}{2} + \frac{b-a}{2} \cos\!\left(\frac{k\pi}{n-1}\right), \qquad k=0,\dots,n-1\]

the extrema of the Chebyshev polynomial \(T_{n-1}\), returned in increasing order. numpy.polynomial.chebyshev.chebpts2() computes these on [-1, 1]; this function affinely maps them onto [a, b]. See Trefethen, Approximation Theory and Approximation Practice, 2013, Ch. 12.

Parameters:

n (int) – Number of nodes.
a (float) – Interval endpoints.
b (float) – Interval endpoints.

Return type:

ndarray

Returns:

ndarray, shape (n,)

Examples

>>> nodes = chebyshev_nodes(5)
>>> float(nodes[0]), float(nodes[-1])
(-1.0, 1.0)
>>> import numpy as np
>>> bool(np.all(np.diff(nodes) > 0))
True

mathematicskit.numerical_analysis.condition_number(matrix)[source]#

2-norm condition number \(\kappa_2(A) = \sigma_{\max}/\sigma_{\min}\), via numpy.linalg.cond().

Used to flag numerically unstable fits/interpolation problems, e.g. PolynomialRegression()’s Vandermonde design matrix, whose condition number grows quickly with polynomial degree. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 7.5, and Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 12.

Parameters:: matrix (ndarray)
Return type:: float
Returns:: float

Examples

>>> import numpy as np
>>> A = np.diag([1.0, 2.0, 100.0])
>>> round(condition_number(A), 4)
100.0

mathematicskit.numerical_analysis.estimate_convergence_order(history, root)[source]#

Estimate an iterative sequence’s empirical order of convergence.

Given the last four iterates \(x_{n-2}, x_{n-1}, x_n, x_{n+1}\) and errors \(e_k = |x_k - \text{root}|\), the order \(p\) satisfying \(e_{k+1} \approx C e_k^p\) is estimated from

\[p \approx \frac{\ln(e_{n+1} / e_n)}{\ln(e_n / e_{n-1})}\]

using the two most recent error ratios. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 2.4 (“Error Analysis for Iterative Methods”), Definition 2.8 and the surrounding discussion.

Parameters:

history (ndarray) – Sequence of iterates, e.g. RootResult.history; needs at least 4 entries.
root (float) – The (numerically) exact root the sequence converges to, e.g. RootResult.root or a known closed-form value.

Return type:

float

Returns:

float – The estimated order \(p\). The textbook values it should reproduce are 1.0 for linearly convergent methods (bisection, fixed-point iteration under a contraction), ~1.618 – the golden ratio – for the secant method’s superlinear convergence, and 2.0 for Newton-Raphson’s quadratic convergence at a simple root.

Examples

>>> import numpy as np
>>> # A synthetic quadratically-convergent sequence: e_{k+1} = e_k^2.
>>> errors = [0.1, 0.01, 0.0001, 1e-8]
>>> history = np.array([1.0 - e for e in errors])
>>> round(estimate_convergence_order(history, 1.0), 2)
2.0

mathematicskit.numerical_analysis.horner(coefficients, x0)[source]#

Evaluate a polynomial and its derivative at x0 by Horner’s scheme.

Writing \(p(x) = a_n x^n + \dots + a_1 x + a_0\) in nested form \((\cdots((a_n x + a_{n-1}) x + a_{n-2}) \cdots) x + a_0\) needs only \(n\) multiplications and \(n\) additions, which is optimal (Ostrowski 1954, Pan 1966). The intermediate values \(b_k = b_{k+1} x_0 + a_k\) are the coefficients of the quotient \(q\) in \(p(x) = (x - x_0) q(x) + p(x_0)\) (synthetic division), and a second pass gives \(p'(x_0) = q(x_0)\).

Parameters:

coefficients (array-like, shape (n + 1,)) – Coefficients, highest power first (numpy.polyval convention).
x0 (float) – Evaluation point.

Return type:

HornerResult

Returns:

HornerResult

Examples

>>> # p(x) = 2x^3 - 6x^2 + 2x - 1 at x = 3
>>> result = horner([2.0, -6.0, 2.0, -1.0], 3.0)
>>> result.value, result.derivative
(5.0, 20.0)
>>> result.quotient.tolist()  # 2x^2 + 0x + 2
[2.0, 0.0, 2.0]

mathematicskit.numerical_analysis.kahan_sum(values)[source]#

Sum floating-point numbers with Kahan’s compensated summation.

A running correction c captures the low-order bits lost when each term is added to the (much larger) partial sum, and feeds them back into the next term:

y = x_i - c
t = s + y
c = (t - s) - y
s = t

Parameters:: values (array-like of float) – The terms to add.
Return type:: float
Returns:: float

Examples

>>> values = [1.0] + [1e-16] * 10_000
>>> naive = 0.0
>>> for v in values:  # left to right: every tiny term is rounded away
...     naive += v
>>> naive
1.0
>>> round(kahan_sum(values), 15)
1.000000000001

mathematicskit.numerical_analysis.lebesgue_constant(nodes, n_eval=2000)[source]#

Estimate the Lebesgue constant of a set of interpolation nodes.

The Lebesgue constant \(\Lambda_n = \max_x \sum_i |L_i(x)|\) (where \(L_i\) are the Lagrange basis polynomials for nodes) bounds how much interpolation error can be amplified relative to the best possible polynomial approximation, \(\|f - p_n\|_\infty \leq (1 + \Lambda_n) \|f - p_n^*\|_\infty\), and so serves as the interpolation problem’s condition number: large \(\Lambda_n\) (as for equally spaced nodes, where it grows exponentially in n) signals the Runge-phenomenon-prone regime, while Chebyshev nodes keep it growing only like \(O(\log n)\). See Trefethen, Approximation Theory and Approximation Practice, 2013, Ch. 15.

Parameters:

nodes (ndarray) – Interpolation nodes.
n_eval (int) – Number of evaluation points spanning [min(nodes), max(nodes)] used to approximate the max.

Return type:

float

Returns:

float – The (approximate) Lebesgue constant.

Examples

>>> import numpy as np
>>> equally_spaced = np.linspace(-1, 1, 15)
>>> chebyshev = np.cos(np.pi * (2 * np.arange(15) + 1) / (2 * 15))
>>> lebesgue_constant(equally_spaced) > lebesgue_constant(chebyshev)
True

mathematicskit.numerical_analysis.remez_minimax(f, degree, a=-1.0, b=1.0, tol=1e-10, max_iter=1000, n_grid=20001)[source]#

Best uniform (minimax) polynomial approximation by the Remez exchange algorithm.

By Chebyshev’s equioscillation theorem, \(p^*\) of degree \(\le n\) minimizes \(\max_{[a,b]} |f - p|\) if and only if the error attains its maximum magnitude with alternating sign at \(n + 2\) points. Starting from the Chebyshev extrema, each step solves the linear system

\[p(x_i) + (-1)^i E = f(x_i), \qquad i = 0, \ldots, n + 1,\]

for \(p\) and the levelled error \(E\), then moves the reference to the alternating extrema of the new error curve (located on a fine grid, then polished with scipy.optimize.minimize_scalar()), until they are equal to within tol. See E. Ya. Remez, “Sur la détermination des polynômes d’approximation de degré donnée,” Comm. Soc. Math. Kharkov (4) 10 (1934), 41-63; L. N. Trefethen, Approximation Theory and Approximation Practice (SIAM, 2013), Ch. 10.

Parameters:

f (Callable) – Vectorized continuous function on \([a, b]\).
degree (int) – Polynomial degree \(n\).
a (float) – Interval endpoints.
b (float) – Interval endpoints.
tol (float) – Stop when \((\max|e| - \min_i |e(x_i)|) / \max|e| < ` `tol\).
max_iter (int) – Maximum number of exchange steps.
n_grid (int) – Size of the grid used to locate error extrema.

Return type:

MinimaxResult

Returns:

MinimaxResult

Examples

>>> # Best cubic to x^4 on [-1, 1] is x^2 - 1/8, error 1/8 (= 2^{-3}).
>>> res = remez_minimax(lambda x: x**4, 3)
>>> bool(np.allclose(res.coefficients, [0.0, 1.0, 0.0, -0.125]))
True
>>> round(res.max_error, 8)
0.125

mathematicskit.numerical_analysis.root_condition_numbers(coefficients, root)[source]#

Relative condition number of a simple root with respect to each coefficient.

Perturbing \(a_k\) to \(a_k(1 + \varepsilon)\) moves a simple root \(r\) by \(\delta r \approx -\varepsilon\,a_k r^k / p'(r)\), so the relative condition number is

\[\kappa_k = \frac{|a_k|\,|r|^k}{|r|\,|p'(r)|}.\]

See Wilkinson, Rounding Errors in Algebraic Processes (1963), Ch. 2; Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed. (SIAM, 2002), Ch. 26.

Parameters:

coefficients (array-like, shape (n + 1,)) – Coefficients, highest power first.
root (float) – A simple, nonzero root of the polynomial.

Return type:

ndarray

Returns:

ndarray, shape (n + 1,) – \(\kappa_k\), in the same order as coefficients.

Examples

>>> # (x - 1)(x - 2) = x^2 - 3x + 2, root 2, p'(2) = 1
>>> root_condition_numbers([1.0, -3.0, 2.0], 2.0).tolist()
[2.0, 3.0, 1.0]

mathematicskit.numerical_analysis.runge_function(x)[source]#

Runge’s classic example, \(f(x) = 1/(1+25x^2)\) on \([-1,1]\).

Parameters:: x (float or array-like of float)
Returns:: float or ndarray

mathematicskit.numerical_analysis.runge_phenomenon_errors(degrees, n_eval=1000)[source]#

Compare equally spaced vs. Chebyshev-node interpolation error for runge_function(), across a range of polynomial degrees.

Demonstrates the Runge phenomenon: max error for equally spaced nodes grows (eventually diverging) with degree, while Chebyshev-node error shrinks.

Parameters:

degrees (array-like of int) – Polynomial degrees (number of nodes minus 1) to test.
n_eval (int) – Number of evaluation points for the max-error estimate.

Returns:

equal_errors, chebyshev_errors (ndarray, shape (len(degrees),)) – Max absolute interpolation error over [-1, 1] for each degree.

Examples

>>> equal_err, cheb_err = runge_phenomenon_errors([5, 10, 15, 20])
>>> bool(equal_err[-1] > equal_err[0])
True
>>> bool(cheb_err[-1] < cheb_err[0])
True

mathematicskit.numerical_analysis.wilkinson_polynomial(n=20)[source]#

Coefficients of Wilkinson’s polynomial \(w(x) = \prod_{k=1}^{n} (x - k)\).

Built with numpy.poly(). For n = 20 several coefficients exceed \(2^{53}\) and are already rounded in double precision, and the roots are so sensitive that changing the \(x^{19}\) coefficient \(-210\) by \(2^{-23}\) sends ten of them into the complex plane (Wilkinson 1959, 1963).

Parameters:: n (int) – Degree (number of roots \(1, \ldots, n\)).
Return type:: ndarray
Returns:: ndarray, shape (n + 1,) – Coefficients, highest power first.

Examples

>>> wilkinson_polynomial(3).tolist()  # (x-1)(x-2)(x-3)
[1.0, -6.0, 11.0, -6.0]
>>> float(wilkinson_polynomial(20)[1])
-210.0