solve

solve(funcs, a, b; verbose=false, returnBoundingBoxes=false, exact=false,
      minBoundingIntervalSize=1e-5, roundoff=53)

Find the roots of a system of functions on the search interval [a, b].

Generates a Chebyshev approximation for each function on the given interval, then uses properties of those approximations to shrink the search interval. When the information in the approximation is insufficient to shrink it further, the interval is subdivided and the search recurses until it zeros in on each root. One point – and optionally a bounding box – is returned per root found.

Arguments

• funcs: Vector of functions to solve simultaneously. Each element is either a callable taking one argument per dimension, or a MultiCheb / MultiPower polynomial.
• a: Vector holding the lower bound of the search interval in each dimension, in dimension order, or a single number to use in every dimension.
• b: The upper bound, likewise. It must be above a in every dimension, and both must be finite; solve throws an ArgumentError otherwise, and when a function cannot be called with one argument per dimension or a polynomial's dimension is not the number of functions.

Keyword arguments

• verbose::Bool = false: Print progress of the approximation and rootfinding to the terminal. Useful for systems that take a long time to solve.
• returnBoundingBoxes::Bool = false: Also return a bounding box for each root.
• exact::Bool = false: Run the transformations on the approximation in higher precision, minimising error at some cost in speed.
• minBoundingIntervalSize::Real = 1e-5: If a root is found whose bounding interval is larger than this in every dimension, the system is solved again on the smaller interval. Smaller values give more accurate roots but increase solve time, and can cause trouble if the functions cannot be evaluated accurately at points close together. The value is absolute while the interval in question lies within [-1, 1], and relative otherwise: for an interval with an endpoint of magnitude greater than 1, it is multiplied by that magnitude in that dimension.
• roundoff::Int = 53: Bits of precision to solve at. 53 selects Float64 and the fast solver; <= 24 selects Float32, <= 11 selects Float16, and anything above 53 switches to BigFloat at that precision.

Returns

A vector of roots, each a point in the search interval. With returnBoundingBoxes = true, returns (roots, boundingBoxes) instead, where each box is a 2 x dim array whose first row holds the lower bound in each dimension and whose second row holds the upper bound.

Examples

julia> using YRoots

julia> solve([(x, y) -> x + y - 0.3, (x, y) -> x - y - 0.1], [-1.0, -1.0], [1.0, 1.0])
1-element Vector{Any}:
 [0.2, 0.09999999999999994]

julia> M1 = MultiPower([0 3 0 2; 1.5 0 7 0; 0 0 4 -2; 0 0 0 1]);

julia> M2 = MultiCheb([0.02 0.31; -0.43 0.19; 0.06 0]);

julia> solve([M1, M2], [-5.0, -5.0], [5.0, 5.0]);