# Solving an equation with two parameters for many pairs of the parameters

I have the equation:

tau == (8.183*10^-4) * Sqrt[Pi] * Sqrt[T]*U^-1.5 * E^(U/T)


where I am trying to solve for U but I have multiple pairs of tau and T. I was wondering how the best way to numerically solve for all the Us I get from the lists of tau and T so I don't have to enter the taus and Ts one at a time.

• This will throw errors about using InverseFunction, which implies that some of the solutions might be thrown away, but you can use Solve on this equation to get a closed form expression for U in terms of ProductLog. After that, just make a Table over the values of tau and T. Dec 21, 2017 at 21:02
• Yeah, it is 8.183 * 10^-4 Dec 21, 2017 at 21:11

Inverse functions of differentiable function satisfy a well-known differential equation. We can exploit that by solving this ODE with NDSolve:

This is the function that maps a $T$ to a mapping $f_T \colon U \mapsto \tau$).

τFun = T \[Function] (U \[Function] (8.183*10^-4)*Sqrt[Pi]*Sqrt[T]*U^-1.5*E^(U/T));


And this is the function that maps each $T$ to $f_T^{-1} \colon \tau \mapsto U$ (or rather its restriction to the interval from τa to τb:

UFun = T \[Function] Block[{S = T, u, τ, U0 = 1.},
NDSolveValue[
Evaluate[{D[u[τ], τ] == 1/τFun[S]'[u[τ]], u[τFun[S][U0]] == U0}],
u,
{τ, τa, τb}],
PrecisionGoal -> 15
];


Sanity check:

f = τFun[1/10];
finv = UFun[1/10];
Plot[f[finv[τ]] - τ, {τ, τa, τa + 1}]


Admittedly, this method is not very precise...

As an application we can plot the family of inverse functions like this

Manipulate[
With[{U = UFun[T]},
Plot[U[τ], {τ, 1., 20.}, PlotRange -> {-0.1, .1}]], {T,
0.1, 1000}]


Final remark

Originally, I aimed at using ParametricNDSolveValue as this is precisly one of the application it was made for. Unfortunately. I did not get it working ParametricNDSolveValue. Maybe somebody else knows how to do it.

If you use an exact exponent (almost always needed for Solve, it seems), you get three solutions, the first of which is real:

Solve[tau == (8.183*10^-4)*Sqrt[Pi]*Sqrt[T]*U^(-3/2)*E^(U/T), U]
(*
{{U -> -1.5 T ProductLog[-(0.00854215/(T^(2/3) tau^(2/3)))]},
{U -> -1.5 T ProductLog[(0.00427107 - 0.00739772 I)/(T^(2/3) tau^(2/3))]},
{U -> -1.5 T ProductLog[(0.00427107 + 0.00739772 I)/(T^(2/3) tau^(2/3))]}}
*)


So {U -> -1.5 T ProductLog[-(0.00854215/(T^(2/3) tau^(2/3)))]} is the real solution (for position T and tau). You can plug in values like this:

{U -> -1.5 T ProductLog[-(0.00854214749414121/(T^(2/3) tau^(2/3)))]} /.
{T -> 1., tau -> 2.}
(*  {U -> 0.00811561}  *)

f[{tau_, T_}] := Evaluate[U /.
First@Solve[tau == (8.183*10^-4)*Sqrt[Pi]*Sqrt[T]*U^-1.5*E^(U/T), U]]
pts = RandomReal[{0, 1}, {25, 2}]
f /@ pts

• This works so much as it ends up being that there are two solutions and this only gets me the lowest solution. Do you know how I would go about getting the higher one? I think it has to do with a problem with ProductLog only gives the lowest order solution. If I try and do higher order solutions of ProductLog, I get imaginary components though I know my solution should be completely real. Dec 21, 2017 at 22:51
• @JamesD GivenT>0 and anticipating U>0 you can set z=U/T and solve for z: Solve[tau == (8.183*10^-4)*Sqrt[Pi]*1/T z^-(3/2) E^z, z].
– Alan
Dec 21, 2017 at 23:59