I am trying to calculate the following integral which contains a parameter. enter image description here

I have used the Integrate and FullSimplify using assumptions but Mathematica fails to produce an analytical solution.

Integrate[((Sin[u]^1.82 + (parameter^(-1))^0.63*Sin[u]^2.45)*
    Sin[parameter + u]^1.82)/(parameter + Sin[u])^2, {u, 0, Pi}, 
 Assumptions -> Inequality[0, Less, parameter, Less, 1]]

Is there another function I can use? If not, what function would you recommend in order to estimate the integral? My end goal is to replace π (integral upper limit) with a second parameter.

  • 1
    $\begingroup$ What do you mean by "estimate the integral"? Mathematica has no problem doing the integral numerically for any value of parameter not equal to zero, but you can't be sure that an analytic solution exists. I tried the Rubi package (apmaths.uwo.ca/~arich) and it didn't give a solution. $\endgroup$ – Jason B. Oct 5 '15 at 9:43

It really depends on the level at which you want to estimate this function. Do you want to end up with a nice closed expression? Do you simply need an expression to model the data? You can't be sure that an analytic solution exists. I tried the Rubi package (apmaths.uwo.ca/~arich) and it didn't give a solution. But you can always fit it to some curve.

{redata, imdata} = 
  Transpose[{{#1, Re[#2]}, {#1, Im[#2]}} & @@@ 
      {parameter, NIntegrate[((Sin[u]^1.82 + (parameter^(-1))^0.63*Sin[u]^2.45)*
          Sin[parameter + u]^1.82)/(parameter + Sin[u])^2, {u, 0,Pi}]}
        , {parameter, 0.005, 1, .005}]]; 

ListLinePlot /@ {redata, imdata}

enter image description here

This suggests to me that we could model both the real and imaginary parts with a multi-exponential decay.

  func = Sum[A[n] Exp[-B[n] x], {n, 1, nexp}];
  params = Flatten[Table[{A[n], B[n]}, {n, nexp}]];
  refit = NonlinearModelFit[redata, func, params, x];
  imfit = NonlinearModelFit[imdata, func, params, x];
  {Show[ListLinePlot[redata, ImageSize -> 300], 
    Plot[refit[x], {x, 0, 1}, PlotStyle -> {Dashed, Red}]],
   ListPlot[Transpose[{redata[[All, 1]], refit["FitResiduals"]}],ImageSize -> 300],
   Show[ListLinePlot[imdata, ImageSize -> 300], 
    Plot[imfit[x], {x, 0, 1}, PlotStyle -> {Dashed, Red}]],
   ListPlot[Transpose[{imdata[[All, 1]], imfit["FitResiduals"]}],ImageSize -> 300]}
  , {nexp, 2, 5}]]

enter image description here

So you could decide to ignore errors smaller than, say 0.5, which would allow you to ignore the imaginary part altogether, and take an answer with three exponential decay terms

$$51.6583 e^{-164.728 x}+15.9181 e^{-23.0781 x}+6.3396 e^{-2.9963 x}$$

A similar strategy would work for the upper integration limit.


Well you might think about a semi-analytical approach:

    lst = Table[{p, 
   NIntegrate[((Sin[u]^1.82 + (p^(-1))^0.63*Sin[u]^2.45)*
       Sin[p + u]^1.82)/(p + Sin[u])^2, {u, 0, Pi}]}, {p, 0.1, 1, 

giving a list of complex values (right? you expect it?) Then you can plot separately the real and imaginary parts and fit them to some reasonable functions. For the real part this gives:


odel1 = a/x^0.5 + b;
ff1 = FindFit[ Re[lst], model1, {a, b}, x]
  Plot[model1 /. ff1, {x, 0, 1}, PlotStyle -> Red]

(*  {a -> 2.76556, b -> -2.47875}  *)

and for the imaginary one - the following:

    model2 = a - b*x^2.5;
ff2 = FindFit[ 
  Transpose[{Transpose[lst][[1]], Im[Transpose[lst][[2]]]}], 
  model2, {a, b}, x]
  ListPlot[Transpose[{Transpose[lst][[1]], Im[Transpose[lst][[2]]]}]],
  Plot[model2 /. ff2, {x, 0, 1}, PlotStyle -> Red]

(*   {a -> -0.0000222685, b -> 0.0152534}   *)

The images which will appear enable you to visually check the fitting quality, and vary it, if necessary. I do not publish them, since I have today a poor internet connection.

Have fun!


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.