Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Calculation of a not so simple integral
Replies: 4   Last Post: Jun 21, 2013 5:45 AM

 Messages: [ Previous | Next ]
 Roland Franzius Posts: 586 Registered: 12/7/04
Calculation of a not so simple integral
Posted: Jun 15, 2013 4:15 AM

(*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[1]:=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[4]:= coefficients={a1,a2,a3,a4,a5};

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

In[2]:= 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[6]:= g[x_,u_,v_]=coefficients.primitives[x,u,v]

Out[6]= 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[7]:= { g[x,2\[Pi],a],f[x,2\[Pi],a]}

Out[7]= {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[11]:= equation[u_,v_]=(0==#&)/@
Cases[
FullSimplifyq@CoefficientList[ Numerator[ Together[
g[x,u,v]-f[x,u,v]]],x],
Except[0]]

Out[11]= {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[12]:= rules[u_,v_] = Solve[equation[u,v],coefficients][[1]]

Out[12]= {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[13]:= equation[u,v]/.rules[u,v]//Simplify

Out[13]= {True,True,True,True,True}

(* Confirming the form of the fraction *)

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

Out[14]= 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[17]:= primInts=(Integrate[ Sin[x/2]^2
#,{x,-\[Infinity],\[Infinity]},Assumptions->a>0]&/@
primitives[x,2\[Pi],a])//TrigToExp//ExpandAll

Out[17]= {\[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[25]:= res=Together@CoefficientList[
Expand[FullSimplify[primInts.coefficients/.rules[2\[Pi],a]]],Exp[-a]].{1,Exp[-a]}

Out[25]= (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[28]:= f[x,2\[Pi],a]

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

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

In[32]:= 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[32]= (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

Date Subject Author
6/15/13 Roland Franzius
6/17/13 Alexei Boulbitch
6/19/13 Dr. Wolfgang Hintze
6/20/13 Dr. Wolfgang Hintze
6/21/13 Roland Franzius