Transform an InterpolatingFunction

you can work on the results as follows:

f1 = n /. sol1[[1]];
f2 = ln /. sol2[[1]];
f3[x_] := Exp[f2[x]]

To check the match between f1 and f3 we can plot:

Plot[{f1[x], f2[x], f3[x]}, {x, 0, 20}, 
 PlotStyle -> {Directive[Red, Dashed, Thickness[0.02]], Automatic, 
   Black}, PlotRange -> All]

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]

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]

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]

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. >>

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.

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}]

And nIF is indeed an interpolating function:

nIF //OutputForm

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