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 $t[\tau]$?

M = 500;
r0 = 3 M;
uppertime = 1996.733812103724;

s = NDSolve[{t'[τ] == Sqrt[1 - 2 M/r0]/(1 - 2 M/r[τ]),
r''[τ] - (M/r[τ]^2)/(1 - 2 M/r[τ])*
  r'[τ]^2 + (M/
    r[τ]^2)*(1 - 2 M/r0)/(1 - 2 M/r[τ]) == 0, t[0] == 0,
r[0] == r0, r'[0] == 0}, {t, r}, {τ, 0, uppertime}]

Plot[Evaluate[{t[τ], r[τ]} /. s], {τ, 0, uppertime},
PlotStyle -> Automatic]

(Cross-posted on the Wolfram forum.)


4 Answers 4


(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)


M = 500;
r0 = 3 M;
uppertime = 1996.733812103724;

sys = {t'[τ] == Sqrt[1 - 2 M/r0]/(1 - 2 M/r[τ]), 
   r''[τ] - (M/r[τ]^2)/(1 - 2 M/r[τ])*
      r'[τ]^2 + (M/r[τ]^2)*(1 - 2 M/r0)/(1 - 2 M/r[τ]) == 0, 
   t[0] == 0, r[0] == r0, r'[0] == 0};

s = NDSolve[sys, {t, r}, {τ, 0, uppertime}] // Quiet;

plt = Plot[Evaluate[{t[τ], r[τ]} /. s], {τ, 0, uppertime},
  PlotRange -> {0, 10000},
  PlotLegends -> Placed[
    {HoldForm[t[τ]], HoldForm[r[τ]]},
    {.4, .6}]]

enter image description here

Extract the points for the first curve

pts = Cases[plt, Line[pts_] :> pts, Infinity][[1]];

The inverse function is the interpolation of the reversed points.

inv = Interpolation[Reverse /@ pts];

Comparing the inverse function with the inverse using ParametricPlot

 ParametricPlot[Evaluate[{t[τ], τ} /. s], {τ, 0, uppertime},
  PlotRange -> {{0, 10000}, Automatic}, 
  AspectRatio -> GoldenRatio ],
 Plot[inv[τ], {τ, 0, 10000},
  PlotStyle -> Directive[
    ColorData[97][2], Dashed]]]

enter image description here

  • 1
    $\begingroup$ The // Quiet makes it look like you're hiding something... $\endgroup$
    – Michael E2
    Commented Jul 19, 2022 at 3:28
  • $\begingroup$ @MichaelE2 - the uppertime value of 1996.733812103724 goes slightly beyond what NDSolve can handle. If you prefer, truncating to 1996.733812 will avoid the warning. $\endgroup$
    – Bob Hanlon
    Commented Jul 19, 2022 at 3:37

InverseFunction also work in interpolating function. InverseFunction[t] /. s[[1]] is what you want.

{a, b} = {t[0] /. s[[1]], t[1996] /. s[[1]]};
inverse1 = InverseFunction[t] /. s[[1]];
Plot[inverse1@x, {x, a, b}, AspectRatio -> 1]

enter image description here

To faster the drawing, we can use FunctionInterpolation.

inverse2 = 
  FunctionInterpolation[InverseFunction[t]@x /. s[[1]], {x, a, b}];
Plot[inverse2@x, {x, a, b}, AspectRatio -> 1]
  • 1
    $\begingroup$ Try inverse1[b - 10]. Why chop off b/uppertime? $\endgroup$
    – Michael E2
    Commented Jul 19, 2022 at 3:36

An other way with no loss of precision is to take the coordinates and values generated by NDSolve, reverse coordinates and values and generate the inverse interpolating function tau[t].

M = 500;
r0 = 3 M;
uppertime = 1996.733812103724;

s = NDSolve[{t'[\[Tau]] == Sqrt[1 - 2 M/r0]/(1 - 2 M/r[\[Tau]]), 
   r''[\[Tau]] - (M/r[\[Tau]]^2)/(1 - 2 M/r[\[Tau]])*
      r'[\[Tau]]^2 + (M/
        r[\[Tau]]^2)*(1 - 2 M/r0)/(1 - 2 M/r[\[Tau]]) == 0, t[0] == 0,
    r[0] == r0, r'[0] == 0}, {t, r}, {\[Tau], 0, uppertime}]

tsol[tau_] = t[tau] /. s[[1]];

taucoordinates = First@tsol["Coordinates"];

tvalues = tsol["ValuesOnGrid"];

tp = Transpose[{tvalues, taucoordinates}];

ListPlot[tp, PlotRange -> All]

tau = Interpolation@tp

Plot[tau[t], {t, 0, tvalues[[-1]]}, PlotRange -> All, 
             GridLines -> Automatic]

enter image description here


We demonstrate how to get a fast accurate inverse. InverseFunction is highly accurate when it works, but it tends to be slow. It depends on FindRoot which is not guaranteed to succeed. In fact, it does not succeed near the singularity of the OP's solution, but that is no great loss, imho. Simple interpolation methods that switch $(x,y)$ to $(y,x)$ are fast and easy, but they may or may not be accurate enough. The ones currently present, Bob Hanlon's and Akku14's, are accurate to about 7 digits on average. Often that will be perfectly fine. I'll give one that's better for cases in which greater accuracy is needed.

Inverting the interpolation table is an easy way to compute the inverse function of an interpolation. It will give exact entries at the interpolation nodes. Between the nodes there may be significant interpolation error, since the inverse function of a polynomial is not always well represented by another polynomial of low degree. It is certainly the case in the OP's example in which the interpolating order is 3. To get a more accurate result, a higher degree interpolant may be used based on the derivatives of the given interpolation; however, this approach meets with limited success. Adding the first and second derivatives yields an improvement of only a order of magnitude or a little more in relative accuracy. We'll look at this method, simply because I had a function for computing the inverse this way, and maybe it would be of interest.

A better way is to construct an interpolant of the inverse function between each step, by picking points between the nodes instead of adding derivative information at the nodes. Adding just the midpoints improves the interpolation by an order of magnitude, similar to adding derivative values at the nodes. One might be tempted, say, to use the Chebyshev points. Indeed, this theoretically should be able to construct an arbitrarily accurate inverse. However, one has to solve y == tIFN[tau] for many, many values of y. This should not be hard to program, but it felt like a lot of work compared to, say, letting NDSolve do almost all the work for you, which we demonstrate below.

NDSolve method

An important strategy is to advance NDSolve successively from one step in the OP's solution to the next with NDSolve`Iterate[]. This will insure that the result is accurate. A solution produced by NDSolve such as the OP's solution has weak singularities (e.g. the third derivative is discontinuous at each step/node). By stopping and starting the integration at each step, we insure that NDSolve does not try to step across a singularity. For instance error estimation is based on the ODE and solution being continuously differentiable through a certain order. Another strategy is to use the "Projection" method to improve each step. Note that the use of InterpolationOrder -> All reduces the interpolation error in NDSolve solutions. Note also that the code gives warnings because of the ridiculously small "IterationSafetyFactor". This small factor is what makes the answer have machine-precision accuracy. We're basically setting an impossibly small error tolerance for the solution. There is no need to increase MaxIterations; in fact, I reduced it from the default 32 to 16, because it is good enough in this case.

tIFN = t /. First@s;  (* OP's solution s *)
{state} = 
   {tau'[t] == 1/tIFN'[tau[t]], tau[tIFN[0]] == 0}, 
   tau, {t, tIFN[0], tIFN[uppertime]}, InterpolationOrder -> All,
   Method -> {"Projection", "Invariants" -> {tIFN@tau[t] == t}, 
     "IterationSafetyFactor" -> 1/10^10, MaxIterations -> 16
   PrecisionGoal -> 6, MaxStepFraction -> 1];
NDSolve`Iterate[state, #] & /@ tIFN@"ValuesOnGrid";
sINV = NDSolve`ProcessSolutions[state]
ndINV = tau /. First[sINV];
{tmin, tmax} = First@ndINV@"Domain";

Inverting the interpolation data including derivatives

If you take all the checks out, the code is fairly simple. (Some like to see careful code; others think it obscures the main idea.) There are a few lines to calculate derivative formulas for inverse functions, which demonstrates one of the cool things about Mathematica: letting it generate symbolically generate the code you need to use in another function. In this case, they are used to create the table for yy. And then the Interpolation@Transpose@... line that creates the answer. The checks for invertibility may be of interest. Like most numerical procedures, they are approximate and not foolproof. For instance, Interpolation[{{{0.},0.,4.},{{1.},1.,4.}}] will pass the checks, but it is not an invertible function. However, NDSolve is unlikely to produce such a pathological case.

(* formulas for the derivatives of the inverse function
 * in terms of the values of the derivatives of the function *)
invD // ClearAll;
mem : invD[n_Integer?Positive] := mem = Evaluate[
     Block[{f, y},
      Derivative[n][InverseFunction[f]][y] /.
        {Derivative[k_][_][_] :> Slot[k]} //
     ] &;

(* inverts an interpolating function if it is (probably)
 * monotonic; otherwise returns InverseFunction[f]  *)
invIF // ClearAll;
invIF::univ = "invIF defined for univariate interpolation functions only.";
invIF::ninv = 
  "Interpolating function is not monotonic. InverseFunction is returned.";
invIF[f_InterpolatingFunction] := With[{x = f@"ValuesOnGrid",
    y = f@"Coordinates",
    yy = Table[
      {k, First@f@"InterpolationOrder" - 1}]},
   (* monotonicity data *)
   With[{monotone1 = DeleteDuplicates@Sign@Differences@x,
     monotone2 = DeleteDuplicates@Sign[f'@"ValuesOnGrid"]},
    (* interpolation of f^(-1) with derivative values *)
    Interpolation@Transpose@{List /@ x,
        Sequence @@ Table[invD[k] @@ yy, 
          {k, First@f@"InterpolationOrder" - 1}]} /;
     (* univariate check *)
     Length@y == 1 &&
      (* monotonicity checks *)
      monotone1 === monotone2 && MatchQ[monotone1, {1} | {-1}]
invIF[f_InterpolatingFunction] := 
    InverseFunction[f]) /; Length@f@"Domain" == 1;
invIF[f_InterpolatingFunction] := 
  (Message[invIF::univ]; Null /; False);

myTau = invIF[tIFN];

Error plots

The main take-away is that the NDSolve method achieves machine-precision accuracy, or thereabouts. We use ndINV and myTau (from invIF[]) from above and the solutions from Bob Hanlon invBH (= inv) and Akku14 invAk (= tau). One should not be too concerned about the increase in error near the singularity of t[τ]. We include the interpolation nodes in the plot to show that the errors in invAk and the invIF solution are zero at those points. The spikes go down way below the plot range, but that cannot be shown without making the comparison of the methods impossible. (I tend to use Akku14's method most often.) We separate the plots to make the comparison easier.

  (invBH@tIFN@tau - tau)/($MachineEpsilon + tau) // RealExponent,
  (invAk@tIFN@tau - tau)/($MachineEpsilon + tau) // RealExponent,
  (myTau@tIFN@tau - tau)/($MachineEpsilon + tau) // RealExponent,
  (ndINV@tIFN@tau - tau)/($MachineEpsilon + tau) // RealExponent},
 {tau, 0, uppertime}, PlotRange -> {-20.1, 0.1},
 PlotPoints -> {200, First@tIFN@"Coordinates"}, 
 PlotStyle -> 
  Thread@Directive[ColorData[97] /@ {3, 2, 5, 1}, Opacity[0.5]],
 PlotLegends -> {"BobHanlon", "Akku14", "invIF", "NDSolve"},
 GridLines -> {None, {-8, -6}}]


Side note on the OP's uppertime: The number looks like it might have been copied from a NDSolve::ndsz message. When you change the limits a or b of {τ, a, b} in NDSolve, you can affect the steps ever so slightly. When I run the OP's code, I get an error/warning:

NDSolve::ndsz: At τ == 1996.7338120977709, step size is effectively zero; singularity or stiff system suspected.

Τhe value for τ is very close to uppertime. If I add a little to uppertime, I will get different stopping time that differs only in the last few digits. Somehow NDSolve is adjusting the step sizes, but I cannot explain exactly what is going on.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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