ClearAll[chebMinimize];
Options[chebMinimize] = {MinRecursion -> 4, MaxRecursion -> 10,
Tolerance -> 1.*^-8,
"SuppressMessages" :> {FindMinimum::lstol, FindMinimum::reged},
"MaxMinima" -> Infinity};
chebMinimize::ncv =
"chebMinimize failed to converge to prescribed accuracy after `1` \
recursive bisections of {`2`, `3`}. chebMinimize obtained `4` for the \
relative error estimate.";
Clear[cheb`minpos, cheb`approx];
(* Returns the indices of the search starting point of the bracketing sample points.
* Note that the Chebyshev abscissae are reversed *)
cheb`minpos[v_] :=
Flatten[SparseArray[PeakDetect[-v, 0, 0, Infinity]]["NonzeroPositions"]];
(* Returns the abscissae and function values of a
* Chebyshev approximation of f of absolute error less than tolerance *)
cheb`$keys = {
"X", (* x-values/abscissae *)
"Y", (* y-values/function values *)
"F", (* the function f *)
"I", (* the interval {x, a, b} *)
"E", (* error *)
"C", (* Chebyshev coefficients *)
"M", (* positions of the local minima *)
"FM" (* FindMinimum results *)
};
cheb`data = AssociationMap[None &, cheb`$keys];cheb`$keys = {
"X", (* x-values/abscissae *)
"Y", (* y-values/function values *)
"F", (* the function f *)
"I", (* the interval {x, a, b} *)
"E", (* error *)
"C", (* Chebyshev coefficients *)
"M", (* positions of the local minima *)
"FM" (* FindMinimum results *)
};
cheb`data = AssociationMap[None &, cheb`$keys]; (* not needed; used for diagnostics *)
cheb`approx[f_, {x_, a_, b_}, tolerance_, minR_: 4, maxR_: 10] :=
Module[{vf, xvals, fvals, error, r = minR,
coeffs}, (* coeffs are not needed; saved in cheb`data *)
If[Quiet@ListQ[f /. x -> {N@a}], (* vf = Listable-ized f *)
vf = Evaluate[f /. x -> #] &,
vf = Map[Evaluate[f /. x -> #] &, #] &;
];
{xvals, fvals} = NestWhile[
Function[{xx, ff},
(* recursive bisection of abscissae *)
With[{new = Rescale[Cos[Pi Range[1., 2^r, 2]/2^r], {-1, 1}, {a, b}]},
r++;
{Riffle[xx, new], Riffle[ff, vf[new]]} (* riffle new values *)
]
] @@ # &,
(r++; {#, vf[#]}) &@ (* initial subdivision of {a,b} *)
Rescale[Cos[Pi Range[0., 2^r]/2^r], {-1, 1}, {a, b}],
(error = (* calculate error from tail of Chebyshev series *)
Total@Abs@#[[-Max[2, r] ;;]]/Max@Abs[#] &[
coeffs = FourierDCT[Last@#, 1]]) > tolerance &,
1,
1 + maxR - minR
];
If[error > tolerance, (* check error *)
Message[chebMinimize::ncv, maxR, a, b, error]
];
cheb`data =
AssociationThread[cheb`$keys -> {xvals, fvals, f, {x, a, b}, error, coeffs, {}, {}}];
{xvals, fvals}
];
chebMinimize[f_, {x_, a_, b_}, opts : OptionsPattern[]] :=
Module[{xvals, fvals, minima, res},
(* get sample points & values *)
{xvals, fvals} =
cheb`approx[f, {x, a, b},
OptionValue[Tolerance],
OptionValue[MinRecursion],
OptionValue[MaxRecursion]];
(* get positions of minima in fvals *)
minima = cheb`minpos[fvals];
If[NumericQ[OptionValue["MaxMinima"]] &&
Length@minima > OptionValue["MaxMinima"],
minima = minima[[
Ordering[fvals[[minima]], OptionValue["MaxMinima"]]
]]
];
cheb`data["M"] = minima;
(* polish minima with FindMinimum *)
res = With[{msgs = OptionValue[chebMinimize, {opts}, "SuppressMessages", Hold] /.
{All | Automatic :> Sequence[], None -> {}}},
msgs /. Hold[m___] :> MinimalBy[First]@
Quiet[
cheb`data["FM"] = FindMinimum[f,
{x,
xvals[[#]],
xvals[[Clip[# + 1, {1, Length@xvals}]]],
xvals[[Clip[# - 1, {1, Length@xvals}]]]}
] & /@ minima,
m]];
MapAt[First, Transpose@res, 1]
];
A similar thing happens with the other term Sin[50000 x + 1]/1000, except that convergence does not start until we have about 50000 points. Note that the level of the pre-convergent phase is about the amplitude of the function, 1/1000. In general, one should expect that for the algorithm to be guaranteed to be successful, the approximation should have an error much less than this minimum amplitude of the oscillations. When this is not known, one can see that a machine-precision approximation is possible. If the oscillations of this order, then it will hopeless anyway, unless the working precision is increased.
The full example is a combination of these two processes. Note that oddness of the second term is reflected in the plateau at 10^-3, which is magnitude of the first term (lower left).
In general, one should expect that for the algorithm to be guaranteed to be successful, the approximation should have an error much less than this minimum amplitude of the oscillations.
When this is not known, one can see that a machine-precision approximation is possible. If the oscillations of this order, then it will hopeless anyway, unless the working precision is increased.
Here is a close=upclose-up of the region where the absolute minimum was found. We can see that there are at least three sample points per oscillation, just enough to identify a good seedseeds for FindMinimum, and that every local minimum has been found.
