MathGroup Archive 2013

[Date Index] [Thread Index] [Author Index]

Search the Archive

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[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



  • 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