ore_algebra.analytic.bounds

Error bounds

Reference:

Marc Mezzarobba, Truncation Bounds for Differentially Finite Series. 2018. Annales Henri Lebesgue, to appear. <https://hal.archives-ouvertes.fr/hal-01817568>

Functions

abs_min_nonzero_root(pol[, tol, min_log, prec])	Compute an enclosure of the absolute value of the nonzero complex root of pol closest to the origin.
bound_polynomials(pols[, Poly])	Compute a common majorant polynomial for the polynomials in pol.
graeffe(pol)	Compute the Graeffe iterate of a polynomial.
growth_parameters(dop)	Find κ, α such that the solutions of dop grow at most like sum(α^n*x^n/n!^κ) ≈ exp(κ*(α·x)^(1/κ)).

Classes

DiffOpBound(dop[, leftmost, special_shifts, …])	A “bound” on the “inverse” of a differential operator at a regular point.
HyperexpMajorant(integrand, num, den[, shift])	A formal power series of the form rat1(z)·exp(int(rat2(ζ), ζ=0..z)), with nonnegative coefficients.
MajorantSeries(variable_name[, cvrad])	A formal power series with nonnegative coefficients
MultiDiffOpBound(majs)	Ad hoc wrapper for passing several DiffOpBounds to StoppingCriterion.
RatSeqBound(nums, den[, exceptional_indices])	Bounds on the tails of a.e.
RationalMajorant(fracs)	A rational power series with nonnegative coefficients, represented as an unevaluated sum of rational fractions with factored denominators.

Exceptions

BoundPrecisionError

class ore_algebra.analytic.bounds.DiffOpBound(dop, leftmost=0, special_shifts=None, max_effort=2, pol_part_len=None, bound_inverse='simple')

A “bound” on the “inverse” of a differential operator at a regular point.

A DiffOpBound may be thought of as a sequence of formal power series

v[n](z) = 1/den(z) · exp ∫ (pol[n](z) + cst·z^ℓ·num[n](z)/den(z))

where

• cst is a real number,
• den(z) is a polynomial with constant coefficients,
• pol[n](z) and num[n](z) are polynomials with coefficients depending on n (given by RatSeqBound objects), and ℓ >= deg(pol[n]).

These series can be used to bound the tails of logarithmic power series solutions y(z) of dop(y) = 0 belonging to a certain subspace (see the documentation of __init__() for details). More precisely, write

y(z) - ỹ(z) = z^λ·(u[0](z)/0! + u[1](z)·log(z)/1! + ···)

where y(z) is a solution of self.dop (in the correct subspace, with λ = self.leftmost) and ỹ(z) is its truncation to O(z^n1). Then, for suitable n0 ∈ ℕ and p(z) ∈ ℝ_+[z], the series ŷ(z) = v[n0](z)·p(z) is a common majorant of u[0], u[1], …

In the typical situation where n0 ≤ n1 and y(z) does not depend on initial conditions “past” n1, a polynomial p(z) of valuation at least n1 with this property can be computed using the methods normalized_residual() and rhs(). Variants with different p hold in more general settings. See the documentation of normalized_residual() and rhs() for more information.

Note that multiplying dop by a rational function changes p(z).

DiffOpBounds are refinable: calling the method refine() will try to replace the parametrized series v[n](z) by one giving tighter bounds. The main effect of refinement is to increase the degree of the polynomial part. This can be done several times, but repeated calls to refine() quickly become expensive.

EXAMPLES:

sage: from ore_algebra import *
sage: from ore_algebra.analytic.bounds import *
sage: Dops, x, Dx = DifferentialOperators()

A majorant sequence:

sage: maj = DiffOpBound((x^2 + 1)*Dx^2 + 2*x*Dx, pol_part_len=0)
sage: print(maj.__repr__(asympt=False))
1.000.../((-x + [0.994...])^2)*exp(int(POL+1.000...*NUM/(-x + [0.994...])^2))
where
POL=0,
NUM=bound(0, ord≤1)*z^0 +
bound((-2.000...*n + 2.000...)/(n^2 - n), ord≤1)*z^1

A majorant series extracted from that sequence:

sage: maj(3)
(1.00... * (-x + [0.994...])^-2)*exp(int([3.000...]...)^2)))

An example with a nontrivial polynomial part:

sage: dop = (x+1)*(x^2+1)*Dx^3-(x-1)*(x^2-3)*Dx^2-2*(x^2+2*x-1)*Dx
sage: DiffOpBound(dop, pol_part_len=3)
1.000.../((-x + [0.99...])^3)*exp(int(POL+1.000...*NUM/(-x + [0.99...])^3)) where
POL=~6.000...*z^0 + ~3.000...*z^1 + ~5.000...*z^2,
NUM=~7.000...*z^3 + ~2.000...*z^4 + ~5.000...*z^5

Refining:

sage: from ore_algebra.examples import fcc
sage: maj = DiffOpBound(fcc.dop5, special_shifts=[(0, 1)])
sage: maj.maj_den
(-z + [0.2047...])^13
sage: maj.maj_den.value()(1/10)
[1.825...]
sage: maj.refine()
sage: maj.maj_den.value()(1/10)
[436565.0...]
sage: maj.majseq_pol_part(10)
[[41.256...], [188.43...]]
sage: maj.refine()
sage: maj.majseq_pol_part(10)
[[41.256...], [188.43...], [920.6...], [4518.9...]]

TESTS:

sage: print(DiffOpBound(Dx - 1, pol_part_len=0).__repr__(asympt=False))
1.000.../(1.000...)*exp(int(POL+1.000...*NUM/1.000...))
where
POL=0,
NUM=bound(-1.000.../n, ord≤1)*z^0

sage: QQi.<i> = QuadraticField(-1)
sage: for dop in [
....:     # orders <= 1 are not supported
....:     Dx, Dx - 1, 1/1000*Dx - 1, i*Dx, Dx + i, Dx^2,
....:     (x^2 + 1)*Dx^2 + 2*x*Dx,
....:     Dx^2 - x*Dx
....: ]:
....:     DiffOpBound(dop)._test()

sage: for l in range(10):
....:     DiffOpBound(Dx - 5*x^4, pol_part_len=l)._test()
....:     DiffOpBound((1-x^5)*Dx - 5*x^4, pol_part_len=l)._test()

sage: from ore_algebra.analytic.bounds import _test_diffop_bound
sage: _test_diffop_bound() # long time

bwrec

File: /home/staff/mkauers/software/own/symbCompPackages/sage/ore_algebra/src/ore_algebra/analytic/bounds.py (starting at line 1771)

normalized_residual(n, last, bwrec_nplus=None, Ring=None)

Compute the “normalized residual” associated to a truncated solution of dop(y) = 0.

Consider a solution

y(z) = z^λ·sum[i,k](y[i,k]·z^i·log(z)^k/k!)

of self.dop(y) = 0, and its truncated series expansion

ỹ(z) = z^λ·sum[i<n,k](y[i,k]·z^i·log(z)^k/k!).

Denote s = deg[z](dop(z,θ)). The equation

monic(dop(z=0,θ))(f(z)) = dop(ỹ)

has at least one solution of the form

f(z) = z^(λ+n)·sum[k](f[k](z)·log(z)^k/k!)

for a finite list [f[0], f[1], …] of polynomials of degree ≤ s-1 (exactly one when none of λ+n, λ+n+1, …, λ+n+s-1 is a root of the indicial polynomial dop(z=0,n)).

This method takes as input the truncation order n and the coefficients

last = [[y[n-1,0], y[n-1,1], ...],
        [y[n-2,0], y[n-2,1], ...],
        ...,
        [y[n-s,0], y[n-s,1], ...]],

and returns a list [f[0], f[1], …] as above.

In order to avoid redundant computations, is possible to pass as additional input the series expansions around λ+n+j (0≤j≤s) of the coefficients of the recurrence operator dop(S⁻¹,ν) = sum[0≤i≤s](b[i](ν)·S⁻¹) associated to dop.

The optional Ring parameter makes it possible to choose the coefficient domain. It is there for debugging purposes.

Warning

The bound holds for the normalized residual computed using the operator self.dop, not the one given as input to __init__. These operators differ by a power-of-x factor, which may change the normalized residual.

EXAMPLES:

sage: from ore_algebra import *
sage: from ore_algebra.analytic.bounds import *
sage: Dops, t, Dt = DifferentialOperators(QQ, 't')

Compute the normalized residual associated to a truncation of the exponential series:

sage: trunc = t._exp_series(5); trunc
1/24*t^4 + 1/6*t^3 + 1/2*t^2 + t + 1
sage: maj = DiffOpBound(Dt - 1)
sage: nres = maj.normalized_residual(5, [[trunc[4]]]); nres
[[-0.00833333333333333 +/- 5.77e-18]]

Check that it has the expected properties:

sage: dopT = (Dt - 1).to_T('Tt'); dopT
Tt - t
sage: dopT.map_coefficients(lambda pol: pol[0])(nres[0]*t^5)
([-0.0416666666666667 +/- 6.40e-17])*t^5
sage: (Dt - 1).to_T('Tt')(trunc).change_ring(CBF)
([-0.0416666666666667 +/- 4.26e-17])*t^5

Note that using Dt - 1 instead of θt - t makes a difference in the result, since it amounts to a division by t:

sage: (Dt - 1)(trunc).change_ring(CBF)
([-0.0416666666666667 +/- 4.26e-17])*t^4

TESTS:

sage: maj = DiffOpBound(Dt^2 + 1)
sage: trunc = t._sin_series(5) + t._cos_series(5)
sage: maj._check_normalized_residual(5, [trunc], ZZ.zero(), QQ)
0

sage: Pol.<n> = CBF[]
sage: Jets.<eta> = CBF[]
sage: bwrec = [n*(n-1), Pol(0), Pol(1)]
sage: bwrec_nplus = [[Jets(pol(5+i)) for pol in bwrec]
....:                for i in [0,1]]
sage: last = [[trunc[4]], [trunc[3]]]
sage: res1 = maj.normalized_residual(5, last, bwrec_nplus)
sage: res2 = maj.normalized_residual(5, last)
sage: len(res1) == len(res2) == 1
True
sage: res1[0] - res2[0]
([+/- ...e-18])*t + [+/- ...e-18]

This operator annihilates t^(1/3)*[1/(1-t)+log(t)^2*exp(t)]+exp(t):

sage: dop = ((81*(-1+t))*t^4*(3*t^6-19*t^5+61*t^4-85*t^3+106*t^2
....: -22*t+28)*Dt^5-27*t^3*(36*t^8-315*t^7+1346*t^6-3250*t^5
....: +4990*t^4-5545*t^3+2788*t^2-1690*t+560)*Dt^4+27*t^2*(54*t^9
....: -555*t^8+2678*t^7-7656*t^6+13370*t^5-17723*t^4+13070*t^3
....: -6254*t^2+4740*t-644)*Dt^3-3*t*(324*t^10-3915*t^9+20871*t^8
....: -67614*t^7+130952*t^6-190111*t^5+180307*t^4-71632*t^3
....: +73414*t^2-26368*t-868)*Dt^2+(243*t^11-3645*t^10+21276*t^9
....: -77346*t^8+163611*t^7-249067*t^6+297146*t^5-83366*t^4
....: +109352*t^3-97772*t^2-4648*t+896)*Dt+162*t^10-1107*t^9
....: +5292*t^8-12486*t^7+17908*t^6-37889*t^5-6034*t^4-1970*t^3
....: +36056*t^2+2044*t-896)

We check that the residuals corresponding to various truncated solutions (both without and with logs, with lefmost=1/3 and leftmost=0) are correctly computed:

sage: n = 20
sage: zero = t.parent().zero()

sage: maj = DiffOpBound(dop, leftmost=0, special_shifts=[(0, 1)])
sage: trunc = [t._exp_series(n), zero, zero]
sage: maj._check_normalized_residual(n, trunc, 0, QQ)
0

sage: maj = DiffOpBound(dop, leftmost=1/3, special_shifts=[(0, 1)])
sage: trunc = [(1-t).inverse_series_trunc(n), zero, zero]
sage: maj._check_normalized_residual(n, trunc, 1/3, QQ)
0
sage: trunc = [(1-t).inverse_series_trunc(n), zero, 2*t._exp_series(n)]
sage: maj._check_normalized_residual(n, trunc, 1/3, QQ)
0

plot(ini=None, pt=None, eps=[1.000000000000000e-50 +/- 7.63e-68], tails=None, color='blue', title=True, intervals=True, **opts)

EXAMPLES:

sage: from ore_algebra import *
sage: from ore_algebra.analytic.bounds import DiffOpBound
sage: from ore_algebra.analytic.local_solutions import LogSeriesInitialValues
sage: Dops, x, Dx = DifferentialOperators()

sage: DiffOpBound(Dx - 1).plot([CBF(1)], CBF(i)/2, RBF(1e-20))
Graphics object consisting of 4 graphics primitives

sage: DiffOpBound(x*Dx^3 + 2*Dx^2 + x*Dx).plot(eps=1e-8)
Graphics object consisting of 4 graphics primitives

sage: dop = x*Dx^2 + Dx + x
sage: DiffOpBound(dop, 0, [(0,2)]).plot(eps=1e-8,
....:       ini=LogSeriesInitialValues(0, {0: (1, 0)}, dop))
Graphics object consisting of 4 graphics primitives

sage: dop = ((x^2 + 10*x + 50)*Dx^10 + (5/9*x^2 + 50/9*x + 155/9)*Dx^9
....: + (-10/3*x^2 - 100/3*x - 190/3)*Dx^8 + (30*x^2 + 300*x + 815)*Dx^7
....: + (145*x^2 + 1445*x + 3605)*Dx^6 + (5/2*x^2 + 25*x + 115/2)*Dx^5
....: + (20*x^2 + 395/2*x + 1975/4)*Dx^4 + (-5*x^2 - 50*x - 130)*Dx^3
....: + (5/4*x^2 + 25/2*x + 105/4)*Dx^2 + (-20*x^2 - 195*x - 480)*Dx
....: + 5*x - 10)
sage: DiffOpBound(dop, 0, [], pol_part_len=4, # not tested
....:         bound_inverse="solve").plot(eps=1e-10)
Graphics object consisting of 4 graphics primitives

TESTS:

sage: DiffOpBound(Dx - 1).plot()
Graphics object consisting of 4 graphics primitives

rhs(n1, normalized_residuals, maj=None)

Compute the right-hand side of a majorant equation valid for each of the given normalized residuals.

INPUT:

A list of normalized residuals q (as computed by normalized_residual() i.e., in particular, with an implicit z^n factor) corresponding to solutions y of self.dop truncated to a same order n1. Optionally, a HyperexpMajorant maj = self(n0) for some n0 ≤ n1.

OUTPUT:

A polynomial (q#)(z) such that, with (q^)(z) = z^n1·(q#)(z),

z·ŷ'(z) - ŷ(z) = (q^)(z)·v[n0](z)·den(z)                     (*)

is a majorant equation of self.dop(ỹ) = Q₀(θ)·q(z) (where Q₀ = monic indicial polynomial) for all q ∈ normalized_residuals. More precisely, if y(z) is a solution of dop(y) = 0 associated to one of the q’s, if ŷ(z) is a solution of (*), and if

|y[λ+n,k]| ≤ ŷ[n]   for   n ≥ n1,   0 ≤ k < mult(n, Q₀),     (**)

then |y[λ+n,k]| ≤ ŷ[n] for all n ≥ n1, k ≥ 0. If maj is omitted, the bound will hold for any choice of n0 ≤ n1 in (*), but may be coarser than that corresponding to a particular n0.

The typical application is with n0 = n1 larger than the n’s corresponding to roots λ+n of Q₀ where the y have nonzero initial values. In this case, one can take

ŷ(z) = v[n0](z)·∫(w⁻¹·(q^)(w)·dw, w=0..z)

and the conditions (**) trivially hold true. (In general, one would need to adjust the integration constant so that they do.)

Note: Some of the above actually makes sense for n1 < n0 as well, provided that (**) also hold for n1 ≤ n < n0 and k ≥ 0 and that q^ be suitably modified.

tail_majorant(n, normalized_residuals)

Bound the tails of order n of solutions of self.dop(y) == 0.

INPUT:

A list of normalized residuals q (as computed by normalized_residual(), i.e., in particular, with an implicit z^n factor) corresponding to solutions y of self.dop truncated to a same order n.

The truncation order n is required to be larger than all n’ such that self.leftmost + n’ is a root of the indicial polynomial of self.dop, and the solution of interest has nonzero initial values there.

OUTPUT:

A HyperexpMajorant representing a common majorant series for the tails y[n:](z) of the corresponding solutions.

class ore_algebra.analytic.bounds.HyperexpMajorant(integrand, num, den, shift=0)

A formal power series of the form rat1(z)·exp(int(rat2(ζ), ζ=0..z)), with nonnegative coefficients.

The fraction rat1 is represented in the form z^shift*num(z)/den(z).

TESTS:

sage: from ore_algebra.analytic.bounds import *
sage: Pol.<z> = RBF[]
sage: one = Pol.one().factor()
sage: den0 = Factorization([(1-z,1)])
sage: integrand = RationalMajorant([(4+4*z, one), (z^2, den0)])
sage: den = Factorization([(1/3-z, 1)])
sage: maj = HyperexpMajorant(integrand, Pol.one(), den); maj
(1.00... * (-z + [0.333...])^-1)*exp(int(4.0...
                                        + 4.0...*z + z^2/(-z + 1.0...)))
sage: maj.cvrad
[0.333...]
sage: maj.bound_series(0, 4)
([336.000...])*z^3 + ([93.000...])*z^2 + ([21.000...])*z + [3.000...]
sage: maj._test()
sage: maj*=z^20
sage: maj
(z^20*1.00... * (-z + [0.333...])^-1)*exp(int(4.000...
                                    + 4.000...*z + z^2/(-z + 1.000...)))
sage: maj._test()

bound_series(rad, ord)

TESTS:

sage: from ore_algebra import *
sage: from ore_algebra.analytic.bounds import DiffOpBound
sage: Dops, x, Dx = DifferentialOperators()
sage: maj = DiffOpBound(Dx-1)(10)
sage: maj.bound(RBF(1000))
[1.97007111401...e+434 +/- ...]

bound_tail_series(rad, ord, n)

Compute a termwise bound on the first ord derivatives of the remainder of order start_index of this majorant series.

exp_part_coeffs_lbounds()

Return a lower bound on the coefficients of the series expansion at the origin of the exponential part.

The output is a generator yielding (index, coeff lower bound) pairs.

It works by finding a “minorant series” for the asymptotically dominant term of the integrand. In the simplest case, the minorant series is generated by a two-term recurrence. (We could get tighter bounds at a higher cost by repeating the process for other terms and multiplying the resulting bounds, and/or using higher-order recurrences in more cases while still avoiding subtractions.)

Like with inv_exp_part_series0, the rational part of the majorant is ignored, and the returned bound makes use of the fact that the integration bounds are chosen so that the exponential part is 1 at 0.

class ore_algebra.analytic.bounds.MajorantSeries(variable_name, cvrad=0)

A formal power series with nonnegative coefficients

bound(rad, rows=1, cols=1, tail=None)

Bound the Frobenius norm of the matrix of the given dimensions whose columns are all equal to

[g(rad), g’(rad), g’‘(rad)/2, …, 1/(r-1)!·g^(r-1)(rad)]

and g is this majorant series. Typically, g(z) is a common majorant series of the elements of a basis of solutions of some differential equation, and the result is then a bound on the corresponding fundamental matrix Y(ζ) for all for all ζ with |ζ| ≤ rad.

bound_series(rad, ord)

Compute a termwise bound on the series expansion of self at rad to order O(x^ord).

More precisely, the upper bound of each interval coefficient is a bound on the corresponding coefficient of self (which itself is a bound on the absolute value of the corresponding coefficient of the series this object is intended to bound).

series(rad, ord)

Compute the series expansion of self at rad to order O(x^ord).

With rad = 0, this returns the majorant series itself. More generally, this can be used to obtain bounds on the derivatives of the series this majorant bounds on disks contained within its disk of convergence.

class ore_algebra.analytic.bounds.MultiDiffOpBound(majs)

Ad hoc wrapper for passing several DiffOpBounds to StoppingCriterion.

(Useful for handling several valuation groups at once.)

class ore_algebra.analytic.bounds.RatSeqBound(nums, den, exceptional_indices={-1: 1})

Bounds on the tails of a.e. rational sequences and their derivatives.

Consider a vector of rational sequences sharing the same denominator,

f(n) = nums(n)/den(n) = [num[i](n)/den(n)]_i.

We assume that den is monic and deg(nums) < deg(den). Let

          ⎧ f^(t)(n),                       n ∉ exceptional_indices,
F[t](n) = ⎨ (d/dX)^t(num(n+X)/(X^{-m}·den(n+X)))(X=0)),
          ⎩                                 exceptional_indices[n] = m

            (the first formula is the specialization to m = 0 of the
            second one),

ref(n, ord) = sum[t=0..ord-1](|n*F[t](n)/t!|),

τ(n) = sum[k=-∞..n](exceptional_indices[k]).

An instance of this class represents a vector of bounds b(n) such that

∀ k ≥ n,   |ref(k, τ(k))| ≤ b(n)  (componentwise).

The bounds are not guaranteed to be nonincreasing.

Such bounds appear as coefficients in the parametrized majorant series associated to differential operators, see DiffOpBound. The ability to bound a sum of derivatives rather than a single rational function is useful to support logarithmic solutions at regular singular points. Vectors of bounds are supported purely for performance reasons, to avoid redundant computations on the indices and denominators.

ALGORITHM:

This version essentially bounds the numerators (from above) and the denominator (from below) separately. This simple strategy works well in the typical case where the indicial equation has only small roots, and makes it easy to share part of the computation over a vector of bounds. In the presence of, e.g., large real roots, however, it is not much better than waiting to get past the largest root.

See the git history for a tighter but more expensive alternative.

EXAMPLES:

sage: Pols.<n> = QQ[]
sage: from ore_algebra.analytic.bounds import RatSeqBound

sage: bnd = RatSeqBound([Pols(1)], n*(n-1)); bnd
bound(1/(n^2 - n), ord≤1)
    = +infinity, +infinity, 1.0000, 0.50000, 0.33333, 0.25000, 0.20000,
    0.16667, ..., ~1.00000*n^-1
sage: [bnd(k)[0] for k in range(5)]
[[+/- inf], [+/- inf], [1.000...], [0.500...], [0.333...]]
sage: bnd._test()
sage: bnd.plot()
Graphics object...

sage: bnd = RatSeqBound([-n], n*(n-3), {-1: 1, 0:1, 3:1}); bnd
bound(-1/(n - 3), ord≤3)
    = 74.767, 74.767, 22.439, 12.000, 12.000, 4.3750, 2.8889, 2.2969,
    ..., ~1.00000
    [74.7...]    for  n <= 0,
    [12.0...]    for  n <= 3
sage: [(bnd.ref(k, bnd.ord(k))[0], bnd(k)[0]) for k in range(5)]
[(0,          [74.767...]),
 (0.750...,   [74.767...]),
 (4.000...,   [22.439...]),
 ([3.000...], [12.000...]),
 (12.000...,  [12.000...])]
sage: bnd._test()

sage: RatSeqBound([n], n, {})
Traceback (most recent call last):
...
ValueError: expected deg(num) < deg(den)

sage: bnd = RatSeqBound([n^5-100*n^4+2], n^3*(n-1/2)*(n-2)^2,
....:                   {0:3, 2:2})
sage: bnd._test(200)
sage: bnd.plot()
Graphics object...

sage: bndvec = RatSeqBound([n, n^2, n^3], (n+1)^4, {-1: 1})
sage: for bnd in bndvec:
....:     bnd._test()

TESTS:

sage: RatSeqBound([Pols(3)], n)(10)
[3.000...]
sage: QQi.<i> = QuadraticField(-1, 'i')
sage: RatSeqBound([Pols(1)], n+i)._test()
sage: RatSeqBound([-n], n*(n-3), {3:1})._test()
sage: RatSeqBound([-n], n*(n-3), {0:1})._test()
sage: RatSeqBound([-n], n*(n-3), {0:1,3:1})._test()
sage: RatSeqBound([CBF(i)*n], n*(n-QQbar(i)), {0:1})._test()
sage: RatSeqBound([QQi['n'](3*i+1)], n + (i-1)/3, {-1: 1})._test()

sage: from ore_algebra.analytic.bounds import _test_RatSeqBound
sage: _test_RatSeqBound() # long time
sage: _test_RatSeqBound(base=QQi, number=3, deg=3) # long time

extend(nums)

Add new sequences to this bound, without changing the rest of the data.

Use with care!

plot(rng=xrange(40), tight=None)

Plot this bound and its reference function.

The vector of nums/bounds must have length one.

EXAMPLES:

sage: from ore_algebra.analytic.bounds import RatSeqBound
sage: Pols.<n> = QQ[]
sage: i = QuadraticField(-1).gen()
sage: bnd = RatSeqBound(
....:     [CBF(i)*n+42], n*(n-3)*(n-i-20), {0:1,3:1})
sage: bnd.plot()
Graphics object consisting of ... graphics primitives
sage: bnd.plot(range(30))
Graphics object consisting of ... graphics primitives

ref(n, ord)

Reference value for a single n.

class ore_algebra.analytic.bounds.RationalMajorant(fracs)

A rational power series with nonnegative coefficients, represented as an unevaluated sum of rational fractions with factored denominators.

The terms must be ordered by decreasing radii of convergence as estimated from the denominator (some numerators may be zero, or, more generally, the numerator and denominator may have common factors).

TESTS:

sage: from ore_algebra.analytic.bounds import *
sage: Pol.<z> = RBF[]
sage: den = Factorization([(1-z, 2), (2-z, 1)])
sage: one = Pol.one().factor()
sage: maj = RationalMajorant([(1 + z, one), (z^2, den), (Pol(0), den)])
sage: maj
1.000... + 1.000...*z + z^2/((-z + 2.000...) * (-z + 1.000...)^2)
sage: maj(1/2)
[2.166...]
sage: maj*(z^10)
1.000...*z^10 + 1.000...*z^11 + z^12/((-z + 2.000...) * (-z + 1.000...)^2)
sage: maj.cvrad
1.000000000000000
sage: maj.series(0, 4)
1.250000000000000*z^3 + 0.5000000000000000*z^2 + z + 1.000000000000000
sage: maj._test()
sage: maj._test(1 + z + z^2/((1-z)^2*(2-z)), return_difference=True)
[0, 0, 0, ...]
sage: maj._test(1 + z + z^2/((1-z)*(2-z)), return_difference=True)
[0, 0, 0, 0.5000000000000000, 1.250000000000000, ...]

sage: RationalMajorant([(Pol(0), den), (Pol(0), den)]).cvrad
[+/- inf]

bound_integral(rad, ord)

Compute a termwise bound on the series expansion of int(self, 0..z) at z = rad, to order O(z^ord).

ore_algebra.analytic.bounds.abs_min_nonzero_root(pol, tol=0.0100000000000000, min_log=-infinity, prec=None)

Compute an enclosure of the absolute value of the nonzero complex root of pol closest to the origin.

INPUT:

• pol – Nonzero polynomial.
• tol – An indication of the required relative accuracy (interval width over exact value). It is currently not guaranteed that the relative accuracy will be smaller than tol.
• min_log – Return a bound larger than 2^min_log. The function may loop if there is a nonzero root of modulus bounded by that value.
• prec – working precision.

ALGORITHM:

Essentially the method of Davenport & Mignotte (1990).

EXAMPLES:

sage: from ore_algebra.analytic.bounds import abs_min_nonzero_root
sage: Pol.<z> = QQ[]
sage: pol = 1/10*z^3 + z^2 + 1/7
sage: sorted(z[0].abs() for z in pol.roots(CC))
[0.377695553183559, 0.377695553183559, 10.0142451007998]

sage: abs_min_nonzero_root(pol)
[0.38 +/- 3.31e-3]

sage: abs_min_nonzero_root(pol, tol=1e-10)
[0.3776955532 +/- 2.41e-11]

sage: abs_min_nonzero_root(pol, min_log=-1.4047042967)
[0.3776955532 +/- 2.41e-11]

sage: abs_min_nonzero_root(pol, min_log=-1.4047042966)
Traceback (most recent call last):
...
ValueError: there is a root smaller than 2^(-1.40470429660000)

sage: abs_min_nonzero_root(pol, tol=1e-50)
[0.3776955531835593496507263902642801708344727099333...]

sage: abs_min_nonzero_root(Pol.zero())
Traceback (most recent call last):
...
ValueError: expected a nonzero polynomial

TESTS:

sage: abs_min_nonzero_root(CBF['x'].one())
+Infinity
sage: abs_min_nonzero_root(CBF['x'].gen())
+Infinity
sage: abs_min_nonzero_root(CBF['x'].gen() - 1/3)
[0.33 +/- 3.34e-3]

An example where the ability to increase the precision is used:

sage: import logging; logging.basicConfig()
sage: logger = logging.getLogger('ore_algebra.analytic.bounds')
sage: logger.setLevel(logging.DEBUG)

sage: abs_min_nonzero_root(z^7 - 2*(1000*z^2-1), tol=1e-20)
DEBUG:ore_algebra.analytic.bounds:failed to bound the roots...
[0.03162277660143379332...]

And one where that used to be the case:

sage: from ore_algebra import DifferentialOperators
sage: from ore_algebra.analytic.bounds import DiffOpBound
sage: Dops, x, Dx = DifferentialOperators()
sage: dop = (x^2 + 10*x + 50)*Dx^2 + Dx + 1
sage: maj = DiffOpBound(dop, bound_inverse="simple")
INFO:ore_algebra.analytic.bounds:bounding local operator (simple, pol_part_len=None, max_effort=2)...

sage: logger.setLevel(logging.WARNING)

ore_algebra.analytic.bounds.bound_polynomials(pols, Poly=None)

Compute a common majorant polynomial for the polynomials in pol.

Note that this returns a majorant, not some kind of enclosure.

TESTS:

sage: from ore_algebra.analytic.bounds import bound_polynomials
sage: Pol.<z> = PolynomialRing(QuadraticField(-1, 'i'), sparse=True)
sage: bound_polynomials([(-1/3+z) << (10^10), (-2*z) << (10^10)])
2.000...*z^10000000001 + [0.333...]*z^10000000000
sage: bound_polynomials([Pol(0)])
0
sage: bound_polynomials([])
Traceback (most recent call last):
...
IndexError: list index out of range

ore_algebra.analytic.bounds.graeffe(pol)

Compute the Graeffe iterate of a polynomial.

EXAMPLES:

sage: from ore_algebra.analytic.bounds import graeffe
sage: Pol.<x> = QQ[]

sage: pol = 6*x^5 - 2*x^4 - 2*x^3 + 2*x^2 + 1/12*x^2^2
sage: sorted(graeffe(pol).roots(CC))
[(0.000000000000000, 2), (0.110618733062304 - 0.436710223946931*I, 1),
(0.110618733062304 + 0.436710223946931*I, 1), (0.547473953628478, 1)]
sage: sorted([(z^2, m) for z, m in pol.roots(CC)])
[(0.000000000000000, 2), (0.110618733062304 - 0.436710223946931*I, 1),
(0.110618733062304 + 0.436710223946931*I, 1), (0.547473953628478, 1)]

TESTS:

sage: graeffe(CIF['x'].zero())
0
sage: graeffe(RIF['x'](-1/3))
0.1111111111111111?

ore_algebra.analytic.bounds.growth_parameters(dop)

Find κ, α such that the solutions of dop grow at most like sum(α^n*x^n/n!^κ) ≈ exp(κ*(α·x)^(1/κ)).

EXAMPLES:

sage: from ore_algebra import *
sage: DiffOps, x, Dx = DifferentialOperators()
sage: from ore_algebra.analytic.bounds import growth_parameters
sage: growth_parameters(Dx^2 + 2*x*Dx) # erf(x)
(1/2, [1.4...])
sage: growth_parameters(Dx^2 + 8*x*Dx) # erf(2*x)
(1/2, [2.8...])
sage: growth_parameters(Dx^2 - x) # Airy
(2/3, [1.0...])
sage: growth_parameters(x*Dx^2 + (1-x)*Dx) # Ei(1, -x)
(1, [1.0...])
sage: growth_parameters((Dx-1).lclm(Dx-2))
(1, [2.0...])
sage: growth_parameters((Dx - x).lclm(Dx^2 - 1))
(1/2, [1.0...])
sage: growth_parameters(x^2*Dx^2 + x*Dx + 1)
(+Infinity, 0)