       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:= fwh[0.1]
Out= 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:= fn[aa_] := NIntegrate[in /. a -> aa, {x, -Infinity,
Infinity}]

we find

In:= fn[0.1]
Out= 0.2518853035259861

in good agreement with fwh[0.1].

The Plot shows the agreement over a broad range of values:

In:= 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

```

• Prev by Date: Global Fitting multiple nonlinear equations to multiple datasets.
• Next by Date: Re: How does one get data out of a TemporalData object?
• Previous by thread: Re: Calculation of a not so simple integral
• Next by thread: Re: Calculation of a not so simple integral