Solving an equation with no analytical form (NIntegrate)

Question

I have the next doble integral:

$$ \int_{0}^{\infty}\int_{0}^{\infty}\frac{130000 E^{-0.36k^2}rSin[kr]Sin[kR]}{(72000+E^{2r})k^2R}drdk$$

In which the integration is made first in "r" and then is made in "k", finally giving a function of "R". The integral cannot be done analytically , so I used the next commands to find a table and the plot:

 fc[k_?NumericQ, r_?NumericQ] =  Simplify[130000*Exp[-0,36k^2]*r*Sin[kr]*Sin[kR]/((72000+Exp[2r])*k^2*R)] //
Expand]

and then:

  tb=Table[{R,NIntegrate[ fc[k,r],{k,0,Infinity},{r,0,Infinity},Method -> {"LevinRule", "Points" -> 5}, PrecisionGoal -> 2, MaxRecursion -> 50]}, {R, 0.1, 20, 0.1}]

So finally I got the next plot:

I Have to find the x-coordinate for a given y-coordinate , namely, solving and equation, say, for a y=20 I should find aproximately x=10, however I don't know how to solve and equation involving a table or this kind of integral, I have to solve this later adding another horrific integral like the one here and do the same thing. Moreover, the integration is taking to much just with one integral !!. Can anybody help with these issues? Thanks in advance :D

Answer 1

One strategy is to make an interpolating function from your tabular data and then invert that. From your code, with some typos removed:

fc[k_?NumericQ, r_?NumericQ] = Simplify[130000*Exp[-0.36 k^2]*r*Sin[k r]*
     Sin[k R]/((72000 + Exp[2 r])*k^2*R)] // Expand;
tb = Table[{R, NIntegrate[fc[k, r], {k, 0, Infinity}, {r, 0, Infinity}, 
     Method -> {"LevinRule", "Points" -> 5}, PrecisionGoal -> 2, 
     MaxRecursion -> 50]}, {R, 0.1, 20, 0.1}];

The interpolating function:

tbInt = Interpolation[tb, InterpolationOrder -> 1]

Visualizing this to make sure everything is OK:

Plot[tbInt[z], {z, 0.1, 20}]

Solve[20 == tbInt[x], x]
{{x -> 8.9264}}

(Note that you get the warning that inverse functions are being used. In this case, the function is invertible, so there is no problem).

Answer 2

fc[k_, r_, R_] = 
  130000*Exp[-36/100 k^2]*r*Sin[k r]*Sin[k R]/((72000 + Exp[2 r])*k^2*R);

The order of the integrations can be exchanged since the integration variables are independent. Either, single integration can be done analytically; however, integrating by k first provides a more manageable form.

f[r_, R_] = Assuming[r >= 0 && R > 0,
  Integrate[fc[k, r, R], {k, 0, Infinity}] //
   Simplify]

(*  (1/((72000 + E^(
   2 r)) R))130000 r (-(3/10)
      E^(-(25/36) (r + R)^2) (-1 + E^((25 r R)/9)) Sqrt[π] - 
   1/4 π (r - R) Erf[(5 (r - R))/6] + 
   1/4 π (r + R) Erf[(5 (r + R))/6])  *)

Plot3D[f[r, R], {r, 0, 20}, {R, 1/10, 20}, PlotRange -> All,
 PlotPoints -> 50,
 AxesLabel -> Automatic,
 LabelStyle -> {Bold, 12}]

The second integration must then be done numerically.

fN[r_?NumericQ, R_?NumericQ] := f[r, R]

tab = Table[
   {R, NIntegrate[fN[r, R], {r, 0, Infinity}]},
   {R, 1/10, 20, 1/10}];

Use Interpolation on the table to define an InterpolatingFunction

if = Interpolation[tab];

Plot[if[R], {R, 1/10, 20},
 PlotStyle -> AbsoluteThickness[3],
 Epilog -> {Yellow, AbsolutePointSize[1], Point[tab]}]

The inverse can be done with InverseFunction

invf1[y_?NumericQ] := InverseFunction[if][y]

invf1[20]

(*  8.91728  *)

However, there is an anomaly in an interval around R == 28

Plot[invf1[y], {y, if[20], if[.1]}]

Another way to invert the function is with FindRoot

invf2[y_?NumericQ] := x /. FindRoot[if[x] == y, {x, 10}]

invf2[20]

(*  8.91728  *)

However, this also has anomalous behavior beginning below 36.

Plot[invf2[y], {y, if[20], if[.1]}]

These methods can be combined to avoid the anomalies.

invf3[y_?NumericQ] := Piecewise[{
   {x /. FindRoot[if[x] == y, {x, 10}], y < 35}},
  InverseFunction[if][y]]

Plot[invf3[y], {y, if[20], if[.1]}]

EDIT: As suggested by @aardvark2012 in a comment to answer by @bills you can also just interpolate the reverse of the entries in the data table (tab)

invf = Interpolation[Reverse /@ tab];

Plot[invf[y], {y, if[20], if[.1]}]