# Obtaining minimax polynomial from an interpolation function

I have an interpolation function f:

rdata = {{0., 1.225}, {2941.1764705882356, 0.9146722941176471}, {5882.35294117647,
0.6683875882352942}, {8823.529411764704, 0.47635299999999997}, {11764.70588235294,
0.32260076470588234}, {14705.882352941177, 0.20289264705882354}, {17647.058823529413,
0.12759417647058824}, {20588.235294117647, 0.08003145882352941}, {23529.411764705885,
0.04988052352941177}, {26470.58823529412, 0.03128139411764706}, {29411.764705882353,
0.019735988235294117}, {32352.941176470587, 0.012494405882352939}, {35294.117647058825,
0.007847365294117647}, {38235.294117647056, 0.0050069258823529405}, {41176.470588235294,
0.00324233}, {44117.64705882353, 0.0021290770588235294}, {47058.82352941176,
0.0014170652941176472}, {50000., 0.000977525}};
f = Interpolation[rdata, InterpolationOrder -> 2, Method -> "Spline"];
hMax = f[[1, 1, 2]];


And I'm trying to get a 5th order polynomial approximation, which would minimize relative error. So I found that FunctionApproximationsMiniMaxApproximation is exactly for this purpose. But when I try to use it I get

Needs["FunctionApproximations"]
MiniMaxApproximation[f[h], {h, {0, hMax}, 5, 0}]


MiniMaxApproximation::extalt: The extrema of the error do not alternate in sign. It may be that MiniMaxApproximation has lost track of the extrema by going too fast. If so try increasing the values in the option Brake. It may be that the WorkingPrecision is insufficient. Otherwise there is an extra extreme value of the error, and MiniMaxApproximation cannot deal with this problem.

Looking at the relative error plot and comparing with the output I see that MiniMaxApproximation missed the local maximum near $h=10000$, thus it got two consecutive errors with the same sign.

I've tried to use Brake option with various values, increasing working precision, but the result is still the same. So what else can I do to finally obtain the minimax polynomial?

• Could it be that a 5th order polynomial simply can't fit the constraints of alternating signs on your curve? I'd experiment with the interpolation order of both the spline and the polynomial, to gain some insight. Jun 8, 2017 at 10:53
• @MikeY I did try to make the order higher and lower, without any results. Jun 8, 2017 at 11:23

Update: Found a simpler way to get a first approximation to minimax polynomial.

## On MiniMaxApproximation[] and the OP's f[h]

The file

FileNameJoin[{\$InstallationDirectory, "AddOns", "Packages",
"FunctionApproximations", "Approximations.m"}]


contains several hundred lines of documentation, discussion and examples. You may wish to consult it for alternatives to what I shall present.

The file warns that functions with singular points can be difficult, if not impossible. I'm not sure if infinite singularities are what they have in mind, or any sort of discontinuities in the function or its first two derivatives (which are used in the algorithm!). The suggestion is to modify the function so that it is more amenable to the algorithm. The example they give is subtracting out a logarithmic singularity, approximating the difference, and adding the singularity back in.

The order-2 interpolation has a discontinuous second derivative:

Plot[{f[h], 10000 f'[h], 10000^2 f''[h]}, {h, 0, 50000},
PlotLabel -> Row[{"Interpolation Order ", f["InterpolationOrder"]}]]


One can bump up the interpolation order to 3 or 4 and get smooth derivatives, but MiniMaxApproximation[] still fails. So I suspected the problem might be arising from a combination of a few issues:

• the oscillation around 14000-16000 (most apparent in f''[h]) or perhaps other less obvious oscillations from the interpolation.
• the exponential decrease in the value of f[h] over the interval (over a factor 1000), which makes minimizing the relative error more difficult.
• the potential rounding error over such a long interval (h^5 can be as big as 3.125*^23).

Since Exp[-c h] can be approximated (for c small enough to avoid underflow but large enough that the value changes by a factor of 1000 or more), it cannot simply be the last two. It seems more likely that the little oscillation introduces extra extrema once the degree of the polynomial is sufficiently high and the polynomial is somewhat close to f[h]. MiniMaxApproximation[] cannot figure out which extrema should be ignored and the function fails. In the OP's case, there is one too many extrema (by the Equioscillation Theorem), which is a pickle because no matter which one is dropped -- and it cannot be one of the end points -- the signs will not alternate. This should have suggested that the starting point was bad, but at first, I did not know you could set the starting point of MiniMaxApproximation[].

mmaOP = MiniMaxApproximation[f[h], {h, {0, hMax}, 5, 0}];

Plot[(f[h] - mmaOP[[2]])/(1 + 0 f[h]), {h, 0, 50000},
GridLines -> {mmaOP[[1]], None},
PlotLabel -> "Error of failed minimax"]


The file also suggests that in difficult cases, it may be possible to approach the minimax polynomial by a sequence of approximations, such as by approximating f[h] to get good starting values for MiniMaxApproximation[] or by starting with a smaller interval and extending it.

## Obtaining the minimax approximant

To get the minimax approximant, I started with an exponential model, since the data looks like a decaying exponential, by forming the interpolation logf of the logarithm of the data rdata. The Chebyshev series approximant mma0 to the interpolation logf was found à la Boyd/Chebfun (see this answer). I then used the MiniMaxApproximation[] to mma0 as the starting point for finding the minimax approximant to f[h].

logf = Interpolation[
MapAt[Log, rdata, {All, 2}],      (* take the Log of the function values *)
InterpolationOrder -> 2, Method -> "Spline"];

Block[{n = 8},  (* Chebyshev interpolation on 2n+1 points *)
cnodes = Rescale[Sin[Pi/2 N@Range[-n, n]/n], {1, -1}, {0, hMax}];
yy = logf@cnodes;
cc = Sqrt[1/n] FourierDCT[yy, 1];
cc[[{1, -1}]] /= 2;
mma0 = Exp[cc[[;; 6]].ChebyshevT[Range[0, 5], 2 h/hMax - 1]]
] // Simplify
(*
1.20576 E^(h (-0.0000759961 - 4.3012*10^-9 h + 9.65032*10^-14 h^2 -
9.43409*10^-19 h^3 + 4.03776*10^-24 h^4))
*)

era = MiniMaxApproximation[mma0, {h, {0, 50000}, 5, 0}];     (* initial approximation *)
mma = MiniMaxApproximation[f[h], era, {h, {0, hMax}, 5, 0}]  (* fed as second argument *)
mma = %[[2, 1]];                                             (* extract the approximant *)
(*
{{0., 10288.8, 18883.7, 31749.7, 41032.1, 47359., 50000.},
{1.27756 - 0.000135637 h + 5.86599*10^-9 h^2 - 1.28016*10^-13 h^3 +
1.40093*10^-18 h^4 - 6.12305*10^-24 h^5,
-0.0429097}}
*)

Plot[{(f[h] - era[[2, 1]])/f[h], (f[h] - mma0)/f[h], (f[h] - mma)/f[h]},
{h, 0, 50000},
PlotStyle -> {Automatic, Automatic, Thick},
PlotLabel -> "Relative error of minimax",
PlotLegends -> {Row[{"Minimax of ", HoldForm@Log[f[h]]}],
Row[{"Exp@Cheby@", HoldForm@Log[f[h]]}],
Row[{"Minimax of ", HoldForm@f[h]}]}]


Note the exponential/Chebyshev/logarithmic approximant is better than the minimax. Of course, you would need an Exp[] function in your actual application to use it. One can also see how the wobbliness of the interpolation is magnified and might create extra extrema; a higher interpolation order produces a somewhat smoother error plot.

I concur with @Anton Antonov that dropping the degree to 3 lets it run. But another issue is that the MiniMaxApproximation[ ] minimizes relative error. Plotting the solution polynomial versus f gives

Which I doubt is acceptable. This distortion is due to working near zero.

If you displace the original data upward, then the 3rd degree fit looks better. you can extract out the bias on the backside.

sdata = rdata + Table[{0, 100}, {rdata // Length}];
g = Interpolation[sdata, InterpolationOrder -> 2, Method -> "Spline"]
hMax = g[[1, 1, 2]];

{pts, {poly, err}} = MiniMaxApproximation[g[h], {h, {0, hMax}, 3, 0}]

Show[Plot[poly, {h, 0, 50000}], Plot[g[x], {x, 0, 50000}]]


Finally, throwing in some denominator rational fraction gets you closer still.

{pts, {poly, err}} = MiniMaxApproximation[g[h], {h, {0, hMax}, 3, 1}]
Show[Plot[poly, {h, 0, 50000}], Plot[g[x], {x, 0, 50000}]]


Unless you have a really good reason to argue from physics that a 5th order polynomial is the right fit, I'd play with the numerator and denominator order to get the best fit, bearing in mind again that the algorithm minimizes relative error, so values near zero will heavily influence the outcome.

EDIT

Here is the non-displaced plot with approximant a rational function {3,2}

Interesting problem, thanks!

• Well, biasing is going against my objective: I was expecting to minimize relative error, see the OP. And using the 5th order polynomial was supposed to be a tradeoff between precision and order. I remember that I had some success in manually fitting a 5th order polynomial to look more or less like a minimax one, so I suspect the solution does exist. Jun 8, 2017 at 14:28
• OK, but if you change the order of the approximant to {3,2}, you will get a decent fit. I'll plot it... Jun 8, 2017 at 14:30

The code below produces the computations without the error MiniMaxApproximation::extalt. Note that I reduced the degree from 5 to 3.

In[231]:=
f = Interpolation[SetPrecision[rdata, 40], InterpolationOrder -> 2,
Method -> "Spline"];
hMax = f[[1, 1, 2]]

Out[232]= 50000.00000000000000000000000000000000000

In[233]:= Needs["FunctionApproximations"]
MiniMaxApproximation[f[h], {h, {0, hMax}, 3, 0},
WorkingPrecision -> 100]

Out[234]= {{0, 24946.78815041953881716742766320715591030,
39806.61641341562991362119284976533120552,
47496.30647414305131656113866599892816232,
50000.00000000000000000000000000000000000}, \
{0.7009223053430526054261990569526498941310 -
0.00004797994292104659156732072897181168384775 h +
1.096193914962884908036115270152706202480*10^-9 h^2 -
8.334805000022467079939309994782606314364*10^-15 h^3,
0.4278185262505693431785507401576743863946}}
`