# Solving an equation with no analytical form (NIntegrate)

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

• There are syntax errors in your code and undefined symbols, such as kr and kR.... Aug 29, 2017 at 17:36

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).

• If the inverse of tbInt is the main point, you could just reverse R and NIntegrate in the definition of tb, interpolate the points to get tbIntinverse (and extract the domain by tbIntinverse["Domain"]). Aug 29, 2017 at 23:01
• Thank you very much Bill s, you helped me a lot. I have another question related. How can I find the solution if there are more than one?. In the interpolation plot we see an one-to-one relation but in the main plot I have to do ( Following your steps) mi plot in fact is not like that. I used solve but it gives me just one solution and sometimes there are two or three solution. I will apprecitate any help . Thanks :D Aug 31, 2017 at 19:34
• Try using Findroot instead of Solve. This lets you specify a starting value, so if you search over many starting values, you can find all the roots. Aug 31, 2017 at 20:53
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]}]
`

• Thank You very much for your help Bob Hanlon, it was very useful for me. As I posted I have to do the same procedure with two more integrals , then add them , plot the result and find for a given y-coordinate its corresponding x. How can I do that if there are several solution?. When I used this , it gives me just un solution . Thank you in advance Aug 31, 2017 at 19:36
• @JhoanPerez - I can not guess what problem you are having. Post a new question which shows what you have tried using the information provided here and explain what difficulty you are having in continuing. Aug 31, 2017 at 19:41
• Okay , I will do it :D. Thanks :D Aug 31, 2017 at 21:00