Transform an InterpolatingFunction

Update-- It turns out that a direct construction from an InterpolatingFunction does not satisfy all cases and that what is needed is a more general function interpolation. It would be nice to have a method that had some sort of automatic precision tracking. FunctionInterpolation is an apparently orphaned system function that seems to make no attempt at tracking accuracy or precision. The best way I know, which I've discussed with some of the numerics people at WRI, is to use NDSolve.

Function interpolation

Basically there are two way you can proceed to interpolate a function f[x], either solve the equation y[x] == f[x] as a differential-algebraic equation with a dummy ODE, or differentiate the equation and solve it as an ODE:

NDSolveValue[{y[x] == f[x], t'[x] == 1, t[a] == a}, y, {x, a, b}] (* DAE *)
NDSolveValue[{y'[x] == D[f[x], x], y[a] == f[a]}, y, {x, a, b}]   (* ODE *)

The DAE method can be highly accurate but it is currently limited to MachinePrecision. The ODE method may be done at any finite precision, but it also has the error from the integration method. In the OP's examples, the function being integrated contains a machine-precision interpolating function, and so one would not expect accuracy or precision better than one half of machine precision; perhaps worse in the ODE method, which would contain the derivative of the interpolating function.

I have used both methods several times on the site:

  • Interpolating a ParametricFunction

  • How to find an approximate solution in a region when Solve or NSolve does not work?

  • How to use `FindRoot` to solve an equation containing a parameter?

  • How do I feed data points into an equation to solve NUMERICALLY?

  • Interpolating an Antiderivative

  • Solve trigonometric equation

  • automatic processing of numerical results in `Plot`

Since FunctionInterpolation does not work well, here is my version of it. It does the DAE method for machine precision and the ODE otherwise.

ClearAll[functionInterpolation];
SetAttributes[functionInterpolation, HoldAll];
functionInterpolation[f_, {x_, a_, b_}, opts : OptionsPattern[NDSolve]] /; 
   MatchQ[WorkingPrecision /. {opts}, WorkingPrecision | MachinePrecision] :=
  Block[{x},
   NDSolveValue[
    {\[FormalY][x] == f, \[FormalT]'[x] == 1, \[FormalT][a] == a},
    \[FormalY], {x, a, b}, opts]
   ];
functionInterpolation[f_, {x_, a_, b_}, opts : OptionsPattern[NDSolve]] :=
  Block[{x},
   NDSolveValue[
    {\[FormalY]'[x] == D[f, x], \[FormalY][a] == f /. x -> a},
    \[FormalY], {x, a, b}, opts]
   ];

OP's second example:

sol3nd = functionInterpolation[Exp[ln[t] /. First@sol4], {t, 0, 100}]
Plot[sol3nd[t] - n[t] /. sol3b, {t, 0, 100}, PlotRange -> All]

Mathematica graphics

If high precision is needed, then, as in the OP's solution sol3b, one needs a higher AccuracyGoal (the relative error of the one above exceeds 1, sometimes by a lot, when the value n[t] is close to zero):

sol3nd = functionInterpolation[Exp[ln[t] /. First@sol4], {t, 0, 100}, AccuracyGoal -> 30];
Plot[(sol3nd[t] - n[t])/n[t] /. sol3b, {t, 0, 100; 61}, PlotRange -> All]

Mathematica graphics

Direct construction (original solution)

A direct way to construct "Exp[sol2]" is to exponentiate the values stored in the interpolating function:

sol3IFN = Interpolation[Transpose[{ln["Grid"], Exp[ln["ValuesOnGrid"]]} /. First@sol2]]

As for accuracy, the interpolation between the grid points will be cubic, whereas for Exp[ln[t]] /. sol2, the interpolation will be the exponential of a cubic. So this method will yield exactly Exp[ln[t]] on the grid points but vary from it in between.

Plot[sol3IFN[t] - Exp[ln[t]] /. sol2, {t, 0, 20}, PlotPoints -> 100]

Mathematica graphics

Comparing with the exact solution shows it is quite accurate:

sol1ex = DSolve[{n'[t] == n[t] (1 - n[t]), n[0] == 10^-5}, n, t];
Plot[sol3IFN[t] - n[t] /. sol1ex, {t, 0, 20}]

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

Mathematica graphics

I did something similar in my answer to Interpolating a ParametricFunction. In that case it was a rescaling. It seems a nonlinear transformation can have problems (see OP's update). Perhaps it is due to changes in convexity.


To obtain the solution to each equation as an InterpolatingFunction, use

sol1 = n /. NDSolve[{n'[t] == n[t] (1 - n[t]), n[0] == 10^-5}, n, {t, 0, 20}][[1]];
sol2 = ln /. NDSolve[{ln'[t] == (1 - E^ln[t]), ln[0] == Log[10^-5]}, ln, {t, 0, 20}][[1]]

Unfortunately, simply constructing Exp[sol2] to match sol1 does not work, as explained in Question 67494. A solution is to construct a Table, apply Exp to each element, and construct a new InterpolatingFunction from the result:

sol3 = Interpolation[Table[{0.1 i, Exp[sol2[0.1 i]]}, {i, 0, 200}]]

The relative error is small.

Plot[(sol3[t]/sol1[t] - 1), {t, 0, 20}]

Relative difference between sol1 and sol3

The solution posted as this one was being written takes an alternative approach, creating a function that applies Exp at each evaluation of sol2, rather than creating a new InterpolatingFunction. Both work.

Update based on my last comment

Actually, the plot above shows the relative difference between sol3 and sol1, but not necessarily the accuracy of sol3. To obtain a good estimate of its accuracy, we must have a high accuracy result for sol1, obtained by reducing the NDSolve step size:

sol1 = n /. NDSolve[{n'[t] == n[t] (1 - n[t]), n[0] == 10^-5}, n, {t, 0, 20}, MaxStepFraction -> 1/4000][[1]]

The value of MaxStepFraction is chosen so that sol1 changes imperceptibly for further decreases. Using the new sol1 result yields a comparison plot,

Relative error between sol1 and sol3

Thus, sol3 is very accurate indeed.


I would augment your good logarithmic NDSolve call to also produce your desired normal interpolating function:

{nIF, lnIF} = NDSolveValue[
    {
    ln'[t] == Piecewise[{{(1 - E^ln[t]), Mod[t, 100] < 60}, {-1, Mod[t, 100] >= 60}}],
    ln[0] == Log[10^-15],
    n'[t] == Exp[ln[t]] ln'[t],
    n[0] == 10^-15
    },
    {n, ln},
    {t, 0, 100}
];

This is similar to @MichaelE2's answer, but I think using 1 NDSolve call allows the step size to be adjusted so that both interpolating functions get the same accuracy and precision goals.

Here is a plot of nIF:

Plot[nIF[x], {x, 0, 100}]

enter image description here

And nIF is indeed an interpolating function:

nIF //OutputForm

InterpolatingFunction[{{0., 100.}}, <>]