       Calculation of a not so simple integral

• To: mathgroup at smc.vnet.net
• Subject: [mg131169] Calculation of a not so simple integral
• From: Roland Franzius <roland.franzius at uos.de>
• Date: Sat, 15 Jun 2013 04:23:08 -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

```(*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

Mathematica  version 6+8 are failing to give a result in a reasonable
time of calculation.
Time constraint does not seem to work

The indefinite integral is calculated readily with an obscure complex
result, the primitive function giving completely useless limits for
x->0,+-oo

The FourierTransform with Limit k->0 seems to give a correct result
compared to numerical integration for values of
a -> 1/20

*)

(* Every rational function mod polynomials can be decomposed in a sum of
fractions with at most quadratic numerators *)

(*Denominator of the integrand*)

In:=f[x_,a_,b_]:= x^-2 (x^2-a^2)^-2 (x^2+b^2)^-2

(*Coefficient list for the quadratic numerators of the partial fractions *)

In:= coefficients={a1,a2,a3,a4,a5};

(* the primitve fractions occurring in the sum decomposition,
linear denominators not allowed for integrability*)

In:= primitives[x_,u_,v_]:=
{ x^-2,
(x^2-u^2)^-2,  x^2 (x^2-u^2)^-2,
(x^2+v^2)^-2,  x^2 (x^2+v^2)^-2};
(* Single first order terms in x cannot occur because of symmetry *)

(* The general form with coefficients inserted)

In:= g[x_,u_,v_]=coefficients.primitives[x,u,v]

Out= a1/x^2 + a2/(-u^2+x^2)^2 + (a3 x^2)/(-u^2+x^2)^2 +
a4/(v^2+x^2)^2 + (a5 x^2)/(v^2+x^2)^2

(* Compare expressions for the given values of the zeros *)

In:= { g[x,2\[Pi],a],f[x,2\[Pi],a]}

Out= {a1/x^2 + a4/(a^2+x^2)^2 +
(a5 x^2)/(a^2+x^2)^2 +
a2/(-4 \[Pi]^2+x^2)^2 +
(a3 x^2)/(-4 \[Pi]^2+x^2)^2,
1/(x^2 (a^2+x^2)^2 (-4 \[Pi]^2+x^2)^2)}

(* Equating to zero all coefficients of the common numerator polynomial
in x*)

In:= equation[u_,v_]=(0==#&)/@
Cases[
FullSimplifyq@CoefficientList[  Numerator[ Together[
g[x,u,v]-f[x,u,v]]],x],
Except]

Out= {0==-1+a1 u^4 v^4,
0==a4 u^4+2 a1 u^4 v^2+(a2-2 a1 u^2) v^4,
0==-2 a4 u^2+(a1+a5) u^4+2 (a2-2 a1 u^2) v^2+(a1+a3) v^4,
0==a2+a4-2 (a1+a5) u^2+2 (a1+a3) v^2,
0==a1+a3+a5}

(* solving the equation for the coefficients*)

In:= rules[u_,v_] = Solve[equation[u,v],coefficients][]

Out= {a1 -> 1/(u^4 v^4),
a2 -> (2 (2 u^2+v^2))/(u^2 (u^2+v^2)^3),
a3 -> -((3 u^2+v^2)/(u^4 (u^2+v^2)^3)),
a4 -> -((2 (u^2+2 v^2))/(v^2 (u^2+v^2)^3)),a5->-((u^2+3
v^2)/(v^4 (u^2+v^2)^3))}

(* Confirming the solution *)

In:= equation[u,v]/.rules[u,v]//Simplify

Out= {True,True,True,True,True}

(* Confirming the form of the fraction *)

In:= g[x,u,v]/.rules[u,v]//FullSimplify

Out= 1/(x^2 (u^2-x^2)^2 (v^2+x^2)^2)

(*Integrating all the primitives with the common factor  Sin[x/2]^2 *)

In:= primInts=(Integrate[ Sin[x/2]^2
#,{x,-\[Infinity],\[Infinity]},Assumptions->a>0]&/@
primitives[x,2\[Pi],a])//TrigToExp//ExpandAll

Out= {\[Pi]/2,
1/(16 \[Pi]),
\[Pi]/4,
\[Pi]/(4 a^3)-(E^-a \[Pi])/(4 a^3)-(E^-a \[Pi])/(4 a^2),
\[Pi]/(4 a)+(E^-a \[Pi])/4-(E^-a \[Pi])/(4 a)}

(* Exapand in terms of Exp[-a] with csoeffients that are polynomials in
the constants a and \[Pi}]*)

In:= res=Together@CoefficientList[
Expand[FullSimplify[primInts.coefficients/.rules[2\[Pi],a]]],Exp[-a]].{1,Exp[-a]}

Out= (E^-a \[Pi] (7 a^2+a^3+12 \[Pi]^2+4 a \[Pi]^2))/4 a^5 (a^2+4
\[Pi]^2)^3)
+
(3 a^7+28 a^5 \[Pi]^2-112 a^2 \[Pi]^4+96 a^3 \[Pi]^4-192
\[Pi]^6+128 a \[Pi]^6)/(64 a^5 \[Pi]^3 (a^2+4 \[Pi]^2)^3)

In:= f[x,2\[Pi],a]

Out= 1/(x^2 (a^2+x^2)^2 (-4 \[Pi]^2+x^2)^2)

(*Compare to result of FourierTransform at k->0*)

In:= fres{1,Exp[-a]}.(
Together@CoefficientList[
Assuming[a>0,
Limit[ Sqrt[2\[Pi]]*FourierTransform[Sin[x/2]^2
f[x,2\[Pi],a],x,k],    k->0]]//
TrigToExp//ExpandAll,Exp[-a]]
)

Out= (E^-a \[Pi] (7 a^2+a^3+12 \[Pi]^2+4 a \[Pi]^2))/(4 a^5 (a^2+4
\[Pi]^2)^3)
+
(3 a^7+28 a^5 \[Pi]^2-112 a^2 \[Pi]^4+96 a^3 \[Pi]^4-192
\[Pi]^6+128 a \[Pi]^6)/(64 a^5 \[Pi]^3 (a^2+4 \[Pi]^2)^3)

(*

My personal conclusion:

Mathematica would profit of a redesign of its integration complex.
It needs much more user friendly prepraring  routines.

More preprocessing by linear decomposition could be easily implemented
to give partial results.

We'd like to have at hands a smooth working application of domain
substitutions, something like:

Integrate[f, {x,a,b}]  -> Integrate[Transform[f], Transform[{x,a,b}],
Transform->function]

The simplest form could be one dimensional form for monotone g

Transform[Integrate[f[x],{x,a,b}],Transformations->{ g, g-1}] :>
Integrate[f[ g[x] ]*D[ g[x], x],{ x, g-1[a],g-1[b]]}]

Decompositions into sums of integrals over decompositions of the domain
into monotone intervalls are easily implemented, too.

*)

--

Roland Franzius

```

• Prev by Date: Kalman Filter for Finance
• Next by Date: Re: What is wrong with 2Pi?
• Previous by thread: Kalman Filter for Finance
• Next by thread: Re: Calculation of a not so simple integral