Message Boards Message Boards

How to find the inverse of an InterpolatingFunction?

Posted 1 year ago

Hello. I am numerically solving the coupled differential equations (for $t[\tau]$ and $r[\tau]$) below. Mathematica outputs two interpolating functions as solutions. I would like to invert one of these solutions - $t[\tau]$. That is, how can I obtain $\tau[t]$ from the interpolating function? Any help would be much appreciated, thanks!

POSTED BY: Aiden G
8 Replies

How accurate do you want the inverse function to be?

How fast do you want it to be?

POSTED BY: Michael Rogers
Posted 1 year ago

It need not be extremely accurate- I’ll be using the inverse to create conceptual plots. Speed is not crucial either. I don’t plan on increasing or carrying the data much farther than this point.

POSTED BY: Aiden G

Then I think Eric's method (inverting the interpolation table) is probably best suited to your needs. InverseFunction[] is much more accurate, about 700-800 times slower, and fails near uppertime (because of the singularity developing in the solution, I would guess).

Here's another way to do the same thing:

tFN = t /. First@s;
inverseFN = tFN[{"ValuesOnGrid", "Grid"}] //
    MapAt[Flatten, #, 2] & //
   Transpose //
  Interpolation;

(Note t["Methods"] /. First@s will show what string arguments may be passed to an InterpolatingFunction; not all methods are defined for every InterpolatingFunction.)

The following is only 60 times slower than interpolating and highly accurate like InverseFunction[]. We can bracket the solution if the interpolating function is monotonic, an implicit assumption in the code:

fInv[if_InterpolatingFunction][y_] := Module[{$x},
   With[{idx = First@Nearest[if["ValuesOnGrid"] -> "Index", y]},
    With[{xvals = Flatten@if["Grid"]},
     $x /. FindRoot[y == if[$x],
       {$x, xvals[[Max[1, idx - 1]]], 
        xvals[[Min[Length@xvals, idx + 1]]]}, Method -> "Brent"]
     ]]];

Usage:

fInv[tFN][0.1]
(*  0.057735  *)

If you want to check monotonicity, here's a simple way:

tdiff = DeleteDuplicates@Sign@Differences@tFN["ValuesOnGrid"];
MatchQ[tdiff, {1} | {-1}]

A slightly more rigorous check:

dtsign = DeleteDuplicates@Sign[tIFN'["ValuesOnGrid"]];
tdiff == dtsign && MatchQ[tdiff, {1} | {-1}]

Caveat: The following passes both tests but is not monotonic. It is unlikely for NDSolve to produce such a pathological case.

foo = Interpolation[{{0, 0}, {1/100, 49/100}, {99/100, 51/100}, {1, 1}}];
foo'[0.5]
(*  -23.7372  *)
POSTED BY: Michael Rogers

How long has "ValuesOnGrid" and "Grid" been accessible like that, and how much of my life have I wasted loading "InterpolatingFunctionAnatomy"??

POSTED BY: Eric Smith

Hi Eric,

I probably found out about it here (2013): https://mathematica.stackexchange.com/questions/28337/whats-inside-interpolatingfunction1-4

However, as pointed out in that post, you can inspect the code for the package with

NotebookOpen[
 FindFile["DifferentialEquations`InterpolatingFunctionAnatomy`"]]

You'll see some of the string arguments used, so they're probably as old as the package.

Many objects (such as InterpolatingFunction, LinearSolveFunction, MeshRegion...) can be queried for such arguments with either obj["Method"], obj["Properties"] or both. Some, like MeshRegion, disappointingly give a long list of mostly unsupported arguments; but others allow you to query useful information.

POSTED BY: Michael Rogers

Thanks for the stack exchange link.

I'm aware of obj["Properties"] but apparently today I just learned about obj["Methods"] thanks again for the information and your posts.

POSTED BY: Eric Smith

Here's another way:

The correct way to extract the points from an interpolating function is to load "DifferentialEquations`InterpolatingFunctionAnatomy`" and use associated functions. But here's a shortcut using ListPlot's ability to plot an interpolating function directly:

solnOfT=t/.s[[1]];
pts=Cases[ListPlot[solnOfT],Point[x:{{_?NumericQ,_?NumericQ}..}]:>x,5][[1]];
pts=Reverse/@pts;

In that last step we just reverse 'x' and 'y' (create an inverse). Now just interpolate. If you have more than one dimension, you're no longer on a square grid and will need to reduce the InterpolationOrder

invTFun = Interpolation[pts, InterpolationOrder -> 1]
invTFun[1000]
(*567.796*)
invTFun = Interpolation[pts, InterpolationOrder -> 3]
invTFun[1000]
(*567.85*)
POSTED BY: Eric Smith

Hi Aiden,

InverseFunction works with InterpolatingFunction.

tFun = t /. First@s
invTFun = InverseFunction[tFun]

invTFun[1000]
(* 567.85 *)

tFun[567.85]
(* 1000. *)
POSTED BY: Rohit Namjoshi
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract