1
$\begingroup$

Folks,

I want to numerically evaluate a InverseFunction h[x] close to the edges x5 and x2 of its domain, with no success up to now. As you can see below, h[x] can only be evaluated sufficiently far from these points, which will not be suitable in my case.

I suspect I can use some hidden option of InverseFunction like MaxIterations. How can I properly perform this evaluation (in particular the NIntegrate)?

The code is

Block[{f, if, g, ig, h, ih, x0, y0, x2, x3, y3, x4, x5, y5, ph2, ph3, 

  cont, norm}, With[{d = 10^-4, r = 10^-8, s = .005, L = 16, n = 10^6},

  if[x_] := -((L s)/(d Log[s/(# r)]))
    ProductLog[-((
    d # Log[s/(# r)] n^(1/L) (s/(# r))^-((d #)/(L s)) )/(L s))] &[x];

  f[x_] := InverseFunction[-((L s)/(d Log[s/(# r)]))
     ProductLog[-((
     d # Log[s/(# r)] n^(1/L) (s/(# r))^-((d #)/(L s)) )/(
     L s))] &][x];

  h[x_] := InverseFunction[-((L s)/(d Log[s/(# r)]))
     ProductLog[-1, -((
     d # Log[s/(# r)] n^(1/L) (s/(# r))^-((d #)/(L s)) )/(
     L s))] &][x];

  ih[x_] := -((L s)/(d Log[s/(# r)]))
    ProductLog[-1, -((
    d # Log[s/(# r)] n^(1/L) (s/(# r))^-((d #)/(L s)) )/(L s))] &[x];

  ig[x_] := InverseFunction[(s Log[1/n] + d # Log[s/(# r)])/(
  d Log[s/(# r)]) &][x];

  g[x_] := (s Log[1/n] + d # Log[s/(# r)])/(d Log[s/(# r)]) &[x];

  y3 = Re[y3] /. FindRoot[if[y3] == ih[y3], {y3, 25}];
  x3 = if[y3];

  x4 = ig[.001];
  x2 = ih[.001];

  y5 = Re[y5] /. FindRoot[ih[y5] == ig[y5], {y5, 10}, MaxIterations -> 1000];
  x5 = ih[y5];

  (* Here my problem is evident *)
  Print[{y5, h[x5], h[1.02 x5], h[1.019 x5]}];
  Print[{x2, h[x2], h[.94 x2], h[.98 x2]}]

  (* Here is the operation I want to perform *)
NIntegrate[1, {a, .001, f[y5]}, {c, x5, ih[a]}, {b, Max[a, if[a]], 
Min[c, h[c]]}]

]]

Thanks!

Improve this question
$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) 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! $\endgroup$ – Dr. belisarius Nov 23 '15 at 19:20
  • $\begingroup$ you should try to provide a more minimal example. The issue here is only that h[x] does not evaluate for 77~<x~<97. Everything else is just distracting. $\endgroup$ – george2079 Nov 23 '15 at 20:11
0
$\begingroup$

The problem is the function you want to invert becomes complex out of range. Here are a couple of approaches:

With[{d = 10^-4, r = 10^-8, s = .005, L = 16, n = 10^6},
 f[y_] := -((L s)/(d Log[s/(y r)])) 
        ProductLog[-1, -((d y
            Log[s/(y r)] n^(1/L) (s/(y r))^-((d y)/(L s)))/(L s))];
 myh[x_] := y /. FindMinimum[Abs[f[y] - x], {y, 10}][[2, 1]];
 myh2[x_] := 
  y /. FindMinimum[Re[#] + 1000 Im[#] &@ (f[y]  - x)^2, {y, 10}][[2,1]]]

 myh[80]

14.2824

 myh2[80]

14.2824

Of course for robustness you should add a check that the min is actually near zero.

I think you'll find these same tricks will speed up finding your other roots as well.

Improve this answer
$\endgroup$
  • $\begingroup$ Thanks, george2079. The function itself works great. The problem is that it does not work inside NIntegrate. The integration I want to perform is NIntegrate[1, {a, .001, f[y5]}, {c, x5, ih[if[a]]}, {b, if[a], myh[c]}]. If I use c or Exp[c] instead of myh[c] an integration is performed, but not with your function yet. I've been trying to use Hold in its definition with ReleaseHold at NIntegrate, with no avail. $\endgroup$ – asoares Nov 25 '15 at 2:38
  • $\begingroup$ ReleaseHold[ NIntegrate[ 1, {a, .001, myf[y5]}, {c, x5, ih[if[a]]}, {b, if[a], Hold[myh[c]]}]] finally worked out, with lots of messages though. Perhaps to ask how to do the same but circumventing them may be worthwhile. Thanks again, george2079. $\endgroup$ – asoares Nov 25 '15 at 12:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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