# Finding an approximate solution to this integral: $\alpha(\phi,r,p,d)=\int_0^\infty w(z,r,p,d)Q(z,r,\phi)dz$

I'm working on a physics problem and encountered a rather complex integral for which I'm trying to find an approximate solution. The integral is of the following form: $$\alpha(\phi,r,p,d)=\int_0^\infty w(z,r,p,d)Q(z,r,\phi)dz$$.

Where $$w(z,r,p,d)=\frac{r(+(r+z)\text{sech }^2(\frac{z+p}{d})-d\tanh{(\frac{z+p}{d})})(z+r\tanh{(\frac{z+p}{d})})}{(r+z)^3d\text{ sech}^2(\frac{p}{d})}$$

And $$Q(z,r,\phi)=-2\exp{\frac{z\phi(6r^2(\phi-1)^2-3rz(\phi^2+\phi-2)+2z^2(\phi^2+\phi+1))}{2r^3(\phi-1)^3}}$$

Numerical integration does just fine and finds a solution for given values of $$r,\phi,d,p$$ but for further work I need an approximate function for $$\alpha(\phi,r,p,d)$$. I'm currently trying to find a solution using AsymptoticIntegrate but this doesn't seem to yield any results:

AsymptoticIntegrate[-((2 E^((z \[Phi] (6 r^2 (-1 + \[Phi])^2 - 3 r z (-2 + \[Phi] + \[Phi]^2) +
2 z^2 (1 + \[Phi] + \[Phi]^2)))/(2 r^3 (-1 + \[Phi])^3))r (-1 + \[Phi]) (d + (r + z) Sech[(P + z)/d]^2 - d Tanh[(P + z)/d]) (z + r Tanh[(P + z)/d]))/(d *Sech[P/d]^2*(r + z)^3)), {z, 0, \[Infinity]},{\[Phi], 0, 3}, Assumptions -> { Re[d] > 0, Re[P] >= 0, Re[r] > 0, 1 >= Re[\[Phi]] >= 0}]


The boundary conditions are:

$$r>0$$

$$d>0$$

$$p\geq0$$

$$0\leq\phi\leq1$$.

Any help is very much appreciated, I'm quite new to mathematica. Thanks a lot.

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 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! Feb 9, 2021 at 15:51

**** See update as per comments below ****

Lots of functions to fit numerical data (I use FindFormula below). Below I start with setting up the integral function with integration limits of $$1\leq z\leq 10$$ for starters and vary just $$\phi$$ generating a table of {$$\phi$$,integralFunction[r,p,d,$$\phi$$]}. Then use FindFormula to fit this data.

Then I would try to vary two variables like $$\phi$$ and $$r$$, generate a table of {$$\phi$$,r,integralFunction[r,p,d,$$\phi$$]} and then try to fit this 3D data. Then try to generate a table of another variable and so on and also increase the integration limits. Kinda brute force but it's something you may wish to work with.

Here is a start just varying $$\phi$$ from 0.01 to 0.99. The red fitted function fits the points nicely.

(*
define functions
*)
wFun[z_, r_, p_, d_] := (
r (d + (r + z) Sech[(z + p)/d]^2 - d Tanh[(z + p)/d]) (z +
r Tanh[(z + p)/d]))/(d (r + z)^3 Sech[p/d]^2);
qFun[z_, r_, \[Phi]_] := -2 Exp[(
z \[Phi] (6 r^2 (\[Phi] - 1)^2 - 3 r z (\[Phi]^2 + \[Phi] - 2) +
2 z^2 (\[Phi]^2 + \[Phi] + 1)))/(2 r^2 (\[Phi] - 1)^3)];
(*
define integral function
*)
myIntFun[r_?NumericQ, p_?NumericQ, d_?NumericQ, \[Phi]_?NumericQ] :=
NIntegrate[wFun[z, r, p, d] qFun[z, r, \[Phi]], {z, 0, 10}];

(*
for now, create table varying just phi from 0.01 to 0.99 and fit a \
formula to this 1D data.
*)
phiTable = Table[
{phi, myIntFun[1, 1, 2, phi]},
{phi, 0.01, 0.99, 0.01}];

(*
find a formula for data
*)

theF[x_] = FindFormula[phiTable, x]
(*
superimpose ListPlot of points with fitted function
*)
lp = ListPlot[phiTable, PlotStyle -> {Black, PointSize[0.01]}]
p1 = Plot[theF[x], {x, 0.01, 0.99}, PlotStyle -> Red]
Show[{lp, p1}]


Update:

As per comments below, FindFormula only works with one variable. So the following uses FindFit and applies data of myIntFun[r,1,2,$$\phi$$] to: $$a+b\phi+c\phi^2+dr+er^2+fr\phi$$ in the range of $$0.1\leq r\leq 0.9$$ and $$0.1\leq \phi\leq 0.9$$ and compares the fitted formula to a ListPointPlot3D of the data points:

phiTable = Table[
{phi, r, myIntFun[r, 1, 2, phi]},
{phi, 0.1, 0.9, 0.05}, {r, 0.1, 0.9, 0.05}];

myParms =
FindFit[Flatten[phiTable, 1],
a + b x + c x^2 + d y + e y^2 + f x y, {a, b, c, d, e, f}, {x, y}];

myFit[x_, y_] = (a + b x + c x^2 + d y + e y^2 + f x y) /. myParms;

my2DPlot =
Plot3D[myFit[x, y], {x, 0, 1}, {y, 0, 1},
PlotStyle -> {Opacity[0.2], Blue}];

lp = ListPointPlot3D[phiTable, BoxRatios -> {1, 1, 1}];
Show[{lp, my2DPlot}]


• Hi! Thanks for your suggestions. FindFormula doesn't work with multiple variables right? I tried it with two and it doesn't have any output. Feb 10, 2021 at 7:50
• Ok. Sorry about that. I updated the code above using FindFit. You'll of course have to develop it further according to your needs. But it's a start. Feb 10, 2021 at 12:40
• Ah I see what you mean, ofcourse that is possible. Thanks a lot! Feb 10, 2021 at 12:59