Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

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


Math Forum » Discussions » Software » comp.soft-sys.math.mathematica

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

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Roland Franzius

Posts: 458
Registered: 12/7/04
Calculation of a not so simple integral
Posted: Jun 15, 2013 4:15 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

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




Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.