Re: Numerically inverting functions?
- To: mathgroup at smc.vnet.net
- Subject: [mg6714] Re: [mg6659] Numerically inverting functions?
- From: Daniel Lichtblau <danl>
- Date: Sat, 12 Apr 1997 02:15:53 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
GB wrote:
>
> Hello. I have a family of functions vv[m]:[0,1] --> R, parametrized
> by m, where m is any real number in [1, Infty]. For each m, vv[m] is
> defined by:
>
> / 1-b-Exp[-b*m]
> | -------------, 0 < b <= 1
> vv[m][b] := < b
> |
> \ m - 1, b == 0
>
> I am interested in the inverse bb[m] of vv[m], namely, the function
> bb[m] such that bb[m][vv[m][b]] == b and vv[m][bb[m][v]] == v. For
> each m, the function vv[m] is strictly monotonic in b, so the desired
> functions bb[m] are all well-defined. Is there a way in Mathematica
> of producing these function bb[m] numerically? I have attempted
> NSolve and InverseFunction with no success.
>
> [N.B. The strange object after the "=" sign above is an ASCII
> rendition of a large "{"]
>
> Thanks in advance,
>
> Gabriel
One way to proceed is by direct solving. This can be done in version 3
of Mathematica.
In[1]:= eq = vv == (1-b-Exp[m*b])/b;
In[2]:= InputForm[inv = Solve[eq, b]]
InverseFunction::ifun:
Warning: Inverse functions are being used. Values may be lost for
multivalued inverses.
Solve::ifun: Inverse functions are being used by Solve, so some
solutions may
not be found.
Out[2]//InputForm=
{{b -> (m/(1 + vv) - ProductLog[(E^(m/(1 + vv))*m)/(1 + vv)])/m}}
The ProductLog is the multi-valued function that inverts y==x*Exp[x] (it
is sometimes called Lambert's W function in the literature). We actually
have a family of solutions of the form
{{b -> (m/(1 + vv) - ProductLog[k, (E^(m/(1 + vv))*m)/(1 + vv)])/m}}
indexed by integers k. Graphing the ProductLog function for the zero
(default) and -1 branches shows real values near the origin. One can see
this using e.g.
Plot[ProductLog[x], {x,-1,1}] and Plot[ProductLog[-1,x], {x,-1,1}].
Warning messages and the plots themselves will show where each branch
runs into trouble. Alternatively, play around with the graph of
x*Exp[x]. Anyway, for many purposes of functional inversion, it is
better to proceed analytically. We show this approach below.
In[3]:= invert[m_] := FindRoot[#==(1-x-Exp[m*x])/x, {x,.7}]&
In[4]:= v[m_] := (1-b-Exp[m*b])/b
Now let's try to see what might be a reasonable value to invert.
In[5]:= Plot[v[3.], {b,.1, 1}]
Out[5]= -Graphics-
A look at the plot reveals that in the relevant range 0<b<=1 we take on
values in the range around -4 to -20. We will invert at -10.
In[6]:= invert[3.][-10]
Out[6]= {x -> 0.634605}
Now let us see if this is consistant with the result from Solve.
In[7]:= inv /. {m->3., vv->-10.}
-17
Out[7]= {{b -> 3.70074 10 }}
Not even close. We'll change the Solve result by hand, using the -1
branch of ProductLog.
In[8]:= inv2 = {b -> (m/(1 + vv) - ProductLog[-1,(E^(m/(1 + vv))*m)/(1 +
vv)])/m};
In[10]:= inv2 /. {m->3., vv->-10.}
Out[10]= {b -> 0.634605}
Much better.
Daniel Lichtblau
Wolfram Research
danl at wolfram.com