I have an expression of the form:

expr = TerribleNumericalFunction[U[p]] == SimpleFunction[p]

The Simple function is a polynomial, and the Terrible function involves elliptic integrals. I want to invert this, solving for U as a function of p as an interpolating function .

Naively, I would plug in p=0 and solve for U[0]. Then I could repeat this for many values of p and interpolate. Is there a better method? I feel like there should be a way to vectorize this process.

As a test case, this expression is way faster to evaluate than mine:

expr = Sqrt[0.04 + U[p] + U[p]^2] == 2.0 + 3.0*p;
soln = Table[0, {i, 0, 99}];
For[p = 0, p < 100, p++,
  soln[[p]] = NSolve[expr, U[p]]

This gets me a list of values (as substitution rules). This also gives 2 roots per p, which is true for my actual case.

How do I best turn my list of substitutions into a list of points for Interpolation? I guess there should really be two interpolating functions... one for each branch of U.

  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – bbgodfrey
    Apr 14, 2015 at 13:49

3 Answers 3


I might do the simple example like this, with a differential-algebraic equation, assuming initial points on each branch are given (U10, U20) -- had the wrong code before:

Clear[U, U1, U2, x, p];
{U10, U20} = U[0] /. NSolve[expr /. p -> 0, U[0]]; (* initial values *)

{sol} = NDSolve[{
  expr /. {{U -> U1}, {U -> U2}},  (* implicit equation*)
  {U1[0], U2[0]} == {U10, U20},    (* initial condition *)
  x'[p] == 1, x[0] == 0},          (*dummy ODE,to force interpolation of U1,U2*)
 {U1, U2}, {p, 0, 100}]
  {{U1 -> InterpolatingFunction[{{0., 100.}}, <>], 
    U2 -> InterpolatingFunction[{{0., 100.}}, <>]}}

Not a particularly interesting equation to plot:

  Plot @@ {Through[{U1, U2}[p]] /. sol, 
    Flatten@{p, U1["Domain"] /. sol},
    AspectRatio -> 1, Frame -> True, Axes -> False},
   Evaluate[expr /. U[p] -> U],
   {p, 0, 100}, {U, -300, 300}]

Mathematica graphics


What you want to do can be accomplished pretty simply by rewriting your expression as a function.

eqn[p_] := Sqrt[0.04 + u + u^2] == 2.0 + 3.0 p

uPts = Table[Flatten[{p, NSolve[eqn[p], u][[All, 1, 2]]}], {p, 0, 99}];
u1 = Interpolation[uPts[[All, {1, 2}]]];
u2 = Interpolation[uPts[[All, {1, 3}]]]

Plot @@ {Through[{u1, u2}[t]], Flatten[{t, u1["Domain"]}], 
           AspectRatio -> 1, PlotLegends -> {"u1", "u2"}}



To set up the points for interpolation, I used two Tables.

pts1 = Table[{p, U[p] /. soln[[p, 1]]}, {p, 0, 99}];
pts2 = Table[{p, U[p] /. soln[[p, 2]]}, {p, 0, 99}];
U1 = Interpolation[pts1];
U2 = Interpolation[pts2];

I still think there should be a better way than a For loop. It is terribly slow.


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.