Compute the preimage of the set Y under function F.
Parameter Y is optional and without it solve
F(x) = 0
for x ∈ X0. Without a starting box
X0 the function is assumed to be univariate and the solution is
searched among all real numbers.
The function must be an interval arithmetic function and may be multivariate, that is, X0 and Y may be vectors or arrays of intervals. The computational complexity largely depends on the dimension of X0.
Return value X is an interval enclosure of
{x ∈ X0 |F(x) ∈ Y}
. The optional
return value X_PAVING is a matrix of non-overlapping interval values
for x in each column, which form a more detailed enclosure for the
preimage of Y. An index vector X_INNER_IDX indicates the
columns of X_PAVING, which are guaranteed to be subsets of the
preimage of Y.
This function uses the set inversion via interval arithmetic (SIVIA) algorithm. That is, X is bisected until F(X) is either a subset of Y or until they are disjoint.
fsolve (@cos, infsup(0, "pi")) ⇒ ans ⊂ [1.5646, 1.5708]
It is possible to use the following optimization options: MaxFunEvals, MaxIter, TolFun, TolX, Vectorize, Contract.
If Vectorize is true
, the function F will be called
with input arguments x(1)
, x(2)
, …,
x(numel (X0))
and each input argument will carry a vector
of different values which shall be computed simultaneously. If Y is a
scalar or vector, Vectorize defaults to true
. If
Vectorize is false
, the function F will receive only
one input argument x at a time, which has the size of X0.
# Solve x1 ^ 2 + x2 ^ 2 = 1 for -3 ≤ x1, x2 ≤ 3, # the exact solution is a unit circle x = fsolve (@(x1, x2) x1.^2 + x2.^2, infsup ([-3; -3], [3; 3]), 1) ⇒ x ⊂ 2×1 interval vector [-1.002, +1.002] [-1.0079, +1.0079]
If Contract is true
, the function F will be called
with Y as an additional leading input argument and, in addition to the
function value, must return a contraction of its input argument(s).
A contraction for input argument x is a subset of x which
contains all possible solutions for the equation
F (x) = Y
. Contractions can be computed using
interval reverse operations, for example with @infsup/absrev
which
contracts the input argument for the absolute value function.
# Solve x1 ^ 2 + x2 ^ 2 = 1 for -3 ≤ x1, x2 ≤ 3 again, # but now contractions speed up the algorithm. function [fval, cx1, cx2] = f (y, x1, x2) # Forward evaluation x1_sqr = x1 .^ 2; x2_sqr = x2 .^ 2; fval = x1_sqr + x2_sqr; # Reverse evaluation and contraction y = intersect (y, fval); # Contract the squares x1_sqr = intersect (x1_sqr, y - x2_sqr); x2_sqr = intersect (x2_sqr, y - x1_sqr); # Contract the parameters cx1 = sqrrev (x1_sqr, x1); cx2 = sqrrev (x2_sqr, x2); endfunction x = fsolve (@f, infsup ([-3; -3], [3; 3]), 1, ... struct ('Contract', true)) ⇒ x = 2×1 interval vector [-1, +1] [-1, +1]
It is possible to combine options Vectorize and Contract. Depending on the combination, function F should have one of the following signatures.
function fval = f (x)
Vectorize = false
and Contract = false
.
function fval = f (x1, x2, …, xN)
Vectorize = true
and Contract = false
.
function [fval, cx] = f (y, x)
Vectorize = false
and Contract = true
.
cx
is a contraction of x
.
function [fval, cx1, cx2, …, cxN] = f (y, x1, x2, …, xN)
Vectorize = true
and Contract = true
.
cx1
is a contraction of x1
, cx2
is a contraction of
x2
, and so on.
Note on performance: The bisection method is a brute-force approach to exhaust the function’s domain and requires a lot of function evaluations. It is highly recommended to use a function F which allows vectorization. For higher dimensions of X0 it is also necessary to use a contraction function.
Accuracy: The result is a valid enclosure.
See also: @infsup/fzero, ctc_union, ctc_intersect, optimset.
The following code
clf hold on grid on axis equal shade = [238 232 213] / 255; blue = [38 139 210] / 255; cyan = [42 161 152] / 255; red = [220 50 47] / 255; # 2D ring f = @(x, y) hypot (x, y); [outer, paving, inner] = fsolve (f, infsup ([-3; -3], [3; 3]), ... infsup (0.5, 2), ... optimset ('TolX', 0.1)); # Plot the outer interval enclosure plot (outer(1), outer(2), shade) # Plot the guaranteed inner interval enclosures of the preimage plot (paving(1, inner), paving(2, inner), blue, cyan); # Plot the boundary of the preimage plot (paving(1, not (inner)), paving(2, not (inner)), red);
Produces the following figure
Figure 1 |
---|
The following code
clf hold on grid on shade = [238 232 213] / 255; blue = [38 139 210] / 255; # This 3D ring is difficult to approximate with interval boxes f = @(x, y, z) hypot (hypot (x, y) - 2, z); [~, paving, inner] = fsolve (f, infsup ([-4; -4; -2], [4; 4; 2]), ... infsup (0, 0.5), ... optimset ('TolX', 0.2)); plot3 (paving(1, not (inner)), ... paving(2, not (inner)), ... paving(3, not (inner)), shade, blue); view (50, 60)
Produces the following figure
Figure 1 |
---|
Package: interval