6
$\begingroup$

I'd like to transform an InterpolatingFunction from NDSolve but can't figure out how. Here's an example. The equation I want to solve is

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

but for numerical reasons (not in this toy example), the log transformed equation

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

works better. How can I transform the result in sol2 to match sol1 by taking E^sol2, without losing any accuracy?

EDIT: A slightly less minimal example:

This is in response to @MichaelE2's comment below. It's still a toy example, but illustrates the problem with interpolating only on the original grid.

This is the original model in natural units. The output is wrong because extra AccuracyGoal is needed.

sol3 = NDSolve[{n'[t] == 
  Piecewise[{{n[t] (1 - n[t]), Mod[t, 100] < 60}, {-n[t], Mod[t, 100] >= 60}}],
  n[0] == 10^-15}, n, {t, 0, 100}][[1]];
Plot[n[t] /. sol3, {t, 0, 100}]

wrong solution

Here's a better solution of the model:

sol3b = NDSolve[{n'[t] == 
  Piecewise[{{n[t] (1 - n[t]), Mod[t, 100] < 60}, {-n[t], Mod[t, 100] >= 60}}],
  n[0] == 10^-15}, n, {t, 0, 100}, AccuracyGoal -> \[Infinity]][[1]];
Plot[n[t] /. sol3b, {t, 0, 100}]

better solution

Here's a solution of the log-transformed model:

sol4 = NDSolve[{ln'[t] == 
  Piecewise[{{(1 - E^ln[t]), Mod[t, 100] < 60}, {-1, Mod[t, 100] >= 60}}],
  ln[0] == Log[10^-15]}, ln, {t, 0, 100}][[1]];
Plot[E^ln[t] /. sol4, {t, 0, 100}]

good solution

Now to try the two different ways to create an exponentiated InterpolatingFunction. This one uses the original grid and leads to a big artifact.

sol5IFN = Interpolation[Transpose[{ln["Grid"], Exp[ln["ValuesOnGrid"]]} /.First@sol4]];
Plot[sol5IFN[t], {t, 0, 100}]

bad transformed function

Looking at the underlying grid makes it obvious what the problem is:

ListPlot[Transpose[{ln["Grid"][[All, 1]], Exp[ln["ValuesOnGrid"]]} /. First@sol4]]

underlying grid

The alternative approach of making your own grid suggested by @bbgodfrey works better:

sol6 = Interpolation[Table[{t, E^ln[t] /. sol4}, {t, 0, 100, 0.1}], InterpolationOrder -> 1];
Plot[sol6[t], {t, 0, 100}]

better interpolatingfunction

UPDATE 2: Higher InterpolationOrder

sol5IFN = Interpolation[Transpose[{ln["Grid"], Exp[ln["ValuesOnGrid"]]} /. First@sol4], InterpolationOrder -> 8];
Plot[sol5IFN[t], {t, 0, 100}, PlotRange -> {0, 1}]

interpolation disaster

|improve this question
$\endgroup$
  • $\begingroup$ @ChrisK After making my last comment on the solution by @MichaelE2, I tried recomputing his sol3IFN with InterpolationOrder -> 8. With that change the corresponding relative error plot looks essentially identical to the second plot in my answer. Perhaps, you could rerun your comparison above with this change. $\endgroup$ – bbgodfrey Dec 9 '14 at 16:53
  • $\begingroup$ @bbgodfrey I added it above -- very wiggly. $\endgroup$ – Chris K Dec 10 '14 at 3:14
  • $\begingroup$ @ChrisK Too bad. Thanks anyway. Good luck with your project. $\endgroup$ – bbgodfrey Dec 10 '14 at 4:10
2
$\begingroup$

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.

|improve this answer
$\endgroup$
  • $\begingroup$ This is potentially great, but I need to check my real application to see if the inaccuracy introduced is too much. $\endgroup$ – Chris K Dec 8 '14 at 19:56
  • $\begingroup$ I was careless in my terminology. The plot above shows the relative difference between sol3 and sol1, not the error in sol3. In fact, the error appears to be in sol1. I just recomputed sol1 to higher accuracy with MaxStepFraction -> 1/500, and the maximum relative difference in the plot dropped to 0.00002. Thus, as you originally suggested, computing and exponentiating sol2 is rather more accurate than computing sol1 directly. $\endgroup$ – bbgodfrey Dec 8 '14 at 21:23
  • $\begingroup$ This works well in my real application. Two things that are nice: 1) if you need the transformed InterpolatingFunction to be a better approximation, you can increase the number of points that went into it and 2) the end values seem to be exact. $\endgroup$ – Chris K Dec 9 '14 at 13:03
6
$\begingroup$

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:

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.

|improve this answer
$\endgroup$
  • $\begingroup$ Yours is a very elegant approach. It appears to be based on mathematica.stackexchange.com/questions/19042, but what is the source of that information? Does Wolfram have some additional documentation on line somewhere? Thanks. $\endgroup$ – bbgodfrey Dec 9 '14 at 11:58
  • $\begingroup$ This is indeed elegant, but it didn't work so well on my real application. That's because my real function has a sharp corner and the new InterpolatingFunction has big distortions due to the cubic interpolation between widely spaced exact points. bbgodfrey's Table approach fixed that issue. But it might be fine for smooth functions. $\endgroup$ – Chris K Dec 9 '14 at 13:07
  • 1
    $\begingroup$ @bbgodfrey Many data structures take an argument sol3IFN["Properties"] that shows arguments you can use. There's also mathematica.stackexchange.com/q/28337 and the tutorial [InterpolatingFunctionAnatomy](reference.wolfram.com/language /tutorial/NDSolvePackages.html). With the Spelunking package, you can find out InterpolatingFunctionAnatomy is implemented as ifun["Grid"] etc. Some users prefer to use InterpolatingFunctionAnatomy functions. I've grown used to using "Properties" for various data structures. $\endgroup$ – Michael E2 Dec 9 '14 at 13:07
  • $\begingroup$ @ChrisK Interesting. That means the corner comes between the computed data points, since the method above has almost zero error on a data point. It doesn't sound like the usual output of NDSolve. I don't suppose you could post the actual sol2 you want interpolated, since your real problem seems to have issues not present in the question? $\endgroup$ – Michael E2 Dec 9 '14 at 13:31
  • $\begingroup$ @MichaelE2 Sure! This comment space it too small to fully describe it, so I'll add it in an edit to my original question. $\endgroup$ – Chris K Dec 9 '14 at 14:56
1
$\begingroup$

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]

enter image description here

|improve this answer
$\endgroup$
  • $\begingroup$ Thanks for the idea! This is what I've done myself, but I really want it in InterpolatingFunction form. $\endgroup$ – Chris K Dec 8 '14 at 19:53
1
$\begingroup$

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

|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.