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:





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

  • 1
    $\begingroup$ 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! $\endgroup$ – Michael E2 Feb 9 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}]

enter image description here


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}]

enter image description here

  • $\begingroup$ 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. $\endgroup$ – M M Feb 10 at 7:50
  • $\begingroup$ 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. $\endgroup$ – Dominic Feb 10 at 12:40
  • $\begingroup$ Ah I see what you mean, ofcourse that is possible. Thanks a lot! $\endgroup$ – M M Feb 10 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.