Reference

contquadgk(f, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm, segbuf=nothing)

Numerically integrate the function f(x) from a to b, and optionally over additional intervals b to c and so on. Keyword options include a relative error tolerance rtol (if atol==0, defaults to sqrt(eps) in the precision of the endpoints), an absolute error tolerance atol (defaults to 0), a maximum number of function evaluations maxevals (defaults to 10^7), and the order of the integration rule (defaults to 7).

Returns a pair (I,E) of the estimated integral I and an estimated upper bound on the absolute error E. If maxevals is not exceeded then E <= max(atol, rtol*norm(I)) will hold. (Note that it is useful to specify a positive atol in cases where norm(I) may be zero.)

The endpoints a et cetera can also be complex (in which case the integral is performed over straight-line segments in the complex plane). If the endpoints are BigFloat, then the integration will be performed in BigFloat precision as well.

More generally, the precision is set by the precision of the integration endpoints (promoted to floating-point types).

The integrand f(x) can return any numeric scalar, vector, or matrix type, or in fact any type supporting +, -, multiplication by real values, and a norm (i.e., any normed vector space). Alternatively, a different norm can be specified by passing a norm-like function as the norm keyword argument (which defaults to norm).

The algorithm is an adaptive Gauss-Kronrod integration technique: the integral in each interval is estimated using a Kronrod rule (2*order+1 points) and the error is estimated using an embedded Gauss rule (order points). The interval with the largest error is then subdivided into two intervals and the process is repeated until the desired error tolerance is achieved.

These quadrature rules work best for smooth functions within each interval, so if your function has a known discontinuity or other singularity, it is best to subdivide your interval to put the singularity at an endpoint. For example, if f has a discontinuity at x=0.7 and you want to integrate from 0 to 1, you should use quadgk(f, 0,0.7,1) to subdivide the interval at the point of discontinuity. The integrand is never evaluated exactly at the endpoints of the intervals, so it is possible to integrate functions that diverge at the endpoints as long as the singularity is integrable (for example, a log(x) or 1/sqrt(x) singularity).

For real-valued endpoints, the starting and/or ending points may be infinite. (A coordinate transformation is performed internally to map the infinite interval to a finite one.)

In normal usage, quadgk(...) will allocate a buffer for segments. You can instead pass a preallocated buffer allocated using alloc_segbuf(...) as the segbuf argument. This buffer can be used across multiple calls to avoid repeated allocation.

roots,rvals = few_poly_roots(c::Vector, vals::Vector, nodes::Vector, n::Int; verb=0)

Return nr polynomial roots roots given by coefficients c, and rvals corresponding polynomial values.

Speed goal is 1 us for nr about 3 and degree-14. May use vals function values at nodes which should fill out [-1,1]. Alternates Newton for next root, then deflation to factor it out.

No failure reporting yet. User should use rvals as quality check.

roots = find_near_roots(vals, nodes, rho, fac, meth)

Returns complex-valued roots of unique polynomial approximant g(z) matching the vector of vals at the vector nodes. The nodes are assumed to be well-chosen for interpolation on [-1,1]. 'roots' are returned in order of increasing (Bernstein) distance from the interval [-1,1].

rho > 0.0 sets the Bernstein ellipse parameter within which to keep roots. Recall that the ellipse for the standard segment [-1,1] has semiaxes cosh(rho) horizontally and sinh(rho) vertically.

fac allows user to pass in a pre-factorized (eg LU) object for the Vandermonde matrix. This accelerates things by 3us for 15 nodes.

meth controls method for polynomial root-finding: "PR" - PolynomialRoots.roots() "PR5" - PolynomialRoots.roots5() degree-5 only (worse perf) "F" - fewpolyroots local attempt

To do:

1. template so compiles for known n (speed up roots? poly eval?)
2. compare Boyd version using Cheby points (needs twice the degree)

Alex Barnett 6/29/23 - 7/4/23 edits by LXVM 8/4/23

compute the vandermonde matrix for a symmetric Kronrod rule

roots_companion(a)

find all complex roots of polynomial a[1]*z^n + a[2]*z^(n-1) + ... + a[n+1]
via companion matrix EVP in O(n^3) time. Similar to MATLAB roots.
Note poly coeffs are in reverse order that in many Julia pkgs.
If the entire C plane is a root, returns [complex(NaN)].

Local reference implementation; superceded by other pkgs.