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