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