Re: Calculation of a not so simple integral
- To: mathgroup at smc.vnet.net
- Subject: [mg131239] Re: Calculation of a not so simple integral
- From: "Dr. Wolfgang Hintze" <weh at snafu.de>
- Date: Thu, 20 Jun 2013 04:45:52 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-outx@smc.vnet.net
- Delivered-to: mathgroup-newsendx@smc.vnet.net
- References: <kpmns8$h52$1@smc.vnet.net>
On 17 Jun., 12:19, Alexei Boulbitch <Alexei.Boulbi... at iee.lu> wrote:
> (*Partial Fractions decomposition, Fourier Integrals
>
> The problem diagnostics:
>
> Calculate
>
> Integrate[ Sin[x/2]^2 / x^2 / (x^2-4*Pi^2 )^2 / (x^2 + a^2)^2 ,
> {x,-oo,oo}, Assumptions-> a>0]
>
> The integrand is nonnegative, has no poles on the real line and decays
> rapidly ~ x^-10 as x->+-oo
>
> Hi, Roland,
>
> I checked your expression, it has a pole on the x axis. Check this
>
> Sin[x/2]^2 / x^2 / (x^2-4*Pi^2 )^2 / (x^2 + a^2)^2//StandardForm
>
> Could it be the case, that you have just written it down in a wrong way?
>
> Alexei BOULBITCH, Dr., habil.
> IEE S.A.
> ZAE Weiergewan,
> 11, rue Edmond Reuter,
> L-5326 Contern, LUXEMBOURG
>
> Office phone : +352-2454-2566
> Office fax: +352-2454-3566
> mobile phone: +49 151 52 40 66 44
>
> e-mail: alexei.boulbi... at iee.lu
Let me explain what I meant with my short remark.
The integral can be calculated completely in Mathematica with the
result
fwh[a_] = (3*a^7*E^a + 28*a^5*E^a*Pi^2 - 112*a^2*(-1 + E^a)*Pi^4 +
16*a^3*(1 + 6*E^a)*Pi^4 - 192*(-1 + E^a)*Pi^6 +
64*a*(1 + 2*E^a)*Pi^6)/(E^a*(64*a^5*Pi^3*(a^2 + 4*Pi^2)^3))
Numerical example
In[62]:= fwh[0.1]
Out[62]= 0.25188490534423524
As described in my contribution the aforementioned thread, the main
idea is decomposition into partial fractions.
The point is, however, that other methods lead to wrong results
Method 1: Calculate the integral directly using Integrate -> fi[a]
Without Assumptions on a it gives a complicated expression involving
functions MeierG.
With the Assumptions->a>0 the result contains only elementary
functions and "looks good".
But it is wrong, as the numerical value at a=0.1 is -141.075.
Method 2: Calculate the indefinite integral and taking limits -> fl[a]
Leads to another well-looking results containing only elementary
functions.
But it is also wrong giving fl[0.1] = - 161.232 (which is diferent
from the wrong result fi[a]).
Method 3: Make decomposition into partial fractions and integrate term
by term
This is the method I proposed and decribed in detail earlier (I don't
repat it here to save space).
This leads to the correct result above.
Why do I consider this solution correct?
Well, mainly because it is in agreement with the numerical integration
using NIntegrate .
Actually, defining
In[27]:= fn[aa_] := NIntegrate[in /. a -> aa, {x, -Infinity,
Infinity}]
we find
In[28]:= fn[0.1]
Out[28]= 0.2518853035259861
in good agreement with fwh[0.1].
The Plot shows the agreement over a broad range of values:
In[32]:= Plot[{fwh[a]/fn[a]}, {a, 0.1, 5}, PlotRange -> {1 - 0.1, 1 +
0.1}]
What do we learn from this study?
First of all, that Mathematica's Integrate-Procedure is very far from
being perfect (as was pointed out already by many users on many
occasions).
It can have the very unpleasant feature to give "good-looking" results
which are competely wrong.
What I strongly suggest is to perform an "obligatory" numerical check
of the result using NIntegrate.
Best regards,
Wolfgang