(* :Title: GaussInt.m *) (* :Author: Frank Kuehnel *) (* :Summary: Mathematicas integration package is very versatile but extremely slow when only the special class of gaussian integrals is to be considered. This package performs the analytic integration of high dimensional gaussian integrals in a managable time frame. *) (* :Package Version: 1.0 *) (* :Copyright: Copyright 2000 Frank Kuehnel. *) (* :History: *) BeginPackage["GaussInt`"] GaussInt::usage = "GaussInt[expr, {vecvars}] integrates analytically over a set of two-dimensional variables {vecvars}. The integration range for each vector variable is the entire infinite plane." Wedge::usage = "Wedge[a, b] gives wedge product of vectors" Begin["`Private`"] Attributes[VecXComp] = {Listable} VecXComp[expr_Symbol] := ToExpression[StringReplace["nx", {"n"->ToString[expr]}]] Attributes[VecYComp] = {Listable} VecYComp[expr_Symbol] := ToExpression[StringReplace["ny", {"n"->ToString[expr]}]] NegIntQ[expr_] := IntegerQ[expr] && TrueQ[expr < 0] NegQ[expr_] := TrueQ[expr < 0] || NegIntQ[expr] (* an algebraic vector qualifier *) AlgVectorQ[expr_] := VectorQ[expr] || MatchQ[expr, a_ + b_ /; AlgVectorQ[a] || AlgVectorQ[b]] || MatchQ[expr, a_ b_ /; AlgVectorQ[a] && !AlgVectorQ[b]] (* wedge product rules *) $WedgeRules = Dispatch[{ Wedge[a_, a_] -> 0, Wedge[a_ + b_, c_] :> Wedge[a, c] + Wedge[b, c], Wedge[a_, b_ + c_] :> Wedge[a, b] + Wedge[a, c], Wedge[a_ c_ /; !AlgVectorQ[c], b_] :> c Wedge[a, b], Wedge[a_, b_ c_ /; !AlgVectorQ[c]] :> c Wedge[a, b], Wedge[a_Symbol, b_Symbol] /; !OrderedQ[{a,b}] :> -Wedge[b, a]} ] (* convert the wedge product to its explicit form *) $ConvWedgeProd = {Wedge[a_Symbol /; VectorQ[a], b_Symbol /; VectorQ[b]] :>\ VecXComp[a] VecYComp[b] - VecYComp[a] VecXComp[b]} (* scalar product rules *) $ScalarProdRules = Dispatch[{ Dot[a_ c_ /; !VectorQ[c], b_] :> c Dot[a, b], Dot[a_, b_ c_ /; !VectorQ[c]] :> c Dot[a, b]} ] (* convert the scalar product into component form *) $ConvScalarProd = {Dot[a_Symbol /; VectorQ[a], b_Symbol /; VectorQ[b]] :>\ VecXComp[a] VecXComp[b] + VecYComp[a] VecYComp[b]} $Convindexmp = {a_Symbol[m] :> VecXComp[a] Sqrt[1/2] + VecYComp[a] I Sqrt[1/2], a_Symbol[p] :> VecXComp[a] Sqrt[1/2] - I Sqrt[1/2] VecYComp[a]} $TrigRules = Dispatch[{ a_ (b_ + c_) :> a b + a c, Sin[a_]^m_Integer/; m>1 :> (1/2-Cos[2 a]/2) Sin[a]^(m-2), Cos[a_]^m_Integer/; m>1 :> (1/2+Cos[2 a]/2) Cos[a]^(m-2), Sin[a_] C[a_] :> Sin[2 a]/2, Sin[a_]Cos[b_] :> Sin[a + b]/2 + Sin[a - b]/2, Cos[a_]Cos[b_] :> Cos[a + b]/2 + Cos[a - b]/2, Sin[a_]Sin[b_] :> Cos[a - b]/2 - Cos[a + b]/2} ] $TrigToExp = { Sin[expr_] :> Exp[I expr]/(2I) - Exp[-I expr]/(2I), Cos[expr_] :> Exp[I expr]/2 + Exp[-I expr]/2} $CollectVar[var_] = {a_. var^2 + b_. var^2 :> (a+b) var^2, a_. var + b_. var :> (a+b) var} $ExpandExpRules = {Exp[expr_] :> Exp[Expand[expr]]} GaussInt[expr_, {vecvars__}] := Block[{newvars, rename, modexpr, result, count=1}, newvars = Map[Unique[q]&, {vecvars}]; Map[(VectorQ[#] ^= True)&, newvars]; rename = Dispatch[MapThread[Rule[#1,#2] &, {{vecvars}, newvars}]]; modexp = expr //. rename //. $WedgeRules //. $ScalarProdRules; (* integrate out easy gaussian integrals drop imaginary parts *) result = Fold[UserInt0[(#1 /. $ExpandExpRules), #2]&,modexp,newvars]; (* Print[result]; *) modexp = Map[( (* Print["integrate next part ", count ]; \ *) (* count+=1; \ *) # //. {UnInt[a_ /; FreeQ[a, UnInt], var_] :> ( \ (* Print["integrate: ", var]; \ *) (* Print[(a //. $ConvScalarProd //. $ConvWedgeProd //. $TrigRules //. $TrigToExp)]; \ *) UserInt1[(a //. $ConvScalarProd //. \ $ConvWedgeProd //. $TrigRules //. $TrigToExp), Join[VecXComp[var],VecYComp[var]]])} \ )&, If[Head[result]===Plus, result, Hold[0]+result]] // ReleaseHold; Map[ClearAll[#]&, newvars]; Return[modexp // Simplify]; ] (* NegQ[expr_] := Negative[expr] *) (* speed up the gaussian integration by testing simple cases (step 1 integration) *) UserInt0[0, var_] := 0 UserInt0[1, var_] := Infinity UserInt0[a_ b_., var_] := a UserInt0[b, var] /; FreeQ[a, var] UserInt0[expr_Plus, var_ ] := Map[UserInt0[#, var]&, expr] UserInt0[UnInt[a_, b_],var_] := UnInt[UserInt0[a, var], b] UserInt0[c_ UnInt[a_, b_],var_] := UnInt[UserInt0[a c, var], b] UserInt0[(Dot[a_,a_]) Exp[c_ Dot[a_,a_]], a_] := If[NegQ[c], Pi/(-c)^2, Infinity] UserInt0[(Dot[a_,a_])^m_Integer Exp[c_ Dot[a_,a_]], a_] := If[NegQ[c], Pi/(-c)^((2m+2)/2) m!, Infinity] UserInt0[(Dot[a_,a_]) Exp[c_ Dot[a_,a_] + d_], a_] := If[NegQ[c], Pi/(-c)^2 Exp[d], Infinity] /; FreeQ[d, a] UserInt0[(Dot[a_,a_])^m_Integer Exp[c_ Dot[a_,a_] + d_], a_] := If[NegQ[c], Pi/(-c)^((2m+2)/2) m! Exp[d], Infinity] /; FreeQ[d,a] UserInt0[dummy_, a_] := UnInt[dummy, {a}] /; FreeQ[dummy, UnInt] (* UserInt0[UnInt[dummy_, a_], b_] := UnInt[dummy, Join[a,{b}]] /; FreeQ[dummy, UnInt] *) (* UserInt0[c_ UnInt[dummy_, a_], b_] := UnInt[c dummy, Join[a,{b}]] /; FreeQ[dummy, UnInt] *) (* more intricate guassian integrals (step 2 integration) *) (* higher dimensional integrals are integrated recursively *) UserInt1[a_, var_List] := Fold[UserInt1[#1,#2]&, a, var] (* linearity rules *) UserInt1[0, var_] := 0 UserInt1[1, var_] := Infinity UserInt1[a_ b_., var_] := a UserInt1[b, var] /; FreeQ[a, var] UserInt1[expr_Plus, var_ ] := Map[UserInt1[#, var]&, expr] UserInt1[a_ (b_ + c_), var_] := UserInt1[a b + a c, var] UserInt1[a_. (b_ + c_)^m_Integer, var_] := UserInt1[a Expand[(b + c)^m], var] (* simple gauss integral *) UserInt1[Exp[a_. var_^2 + b_.], var_] := If[NegQ[a], (Pi/-a)^(1/2) Exp[b], Infinity] /; FreeQ[b, var] UserInt1[Exp[a_. var_^2 + b_. var_ + c_.], var_] := If[NegQ[a], (Pi/-a)^(1/2) Exp[ b^2/(-4a) + c // Expand], Infinity] /; FreeQ[{b,c}, var] (* gauss integrals with powers *) UserInt1[Exp[a_. var_^2 + b_.] var_ , var_] := 0 /; FreeQ[b, var] (* 3.462 - 2 special *) UserInt1[Exp[a_. var_^2 + b_. var_ + c_.] var_, var_] := If[NegQ[a], b/(-2a) (Pi/-a)^(1/2) Exp[b^2/(-4a) + c // Expand], Indetermined] /; FreeQ[{b,c}, var] (* 3.461 - 2 *) UserInt1[Exp[a_. var_^2 + b_.] var_^m_Integer /; m>0 , var_] := If[OddQ[m], 0, If[NegQ[a], Gamma[(m+1)/2]/(-a)^((m+1)/2) Exp[b], Infinity]] /; FreeQ[b, var] (* 3.462 - 2 general *) UserInt1[Exp[a_. var_^2 + b_. var_ + c_.] var_^m_Integer, var_] := If[NegQ[a], (Pi/-a)^(1/2) Exp[c] (D[Exp[qprot^2/(4 pprot)], {qprot, m}] //. {qprot->b, pprot->(-a)}), Indetermined]\ /; FreeQ[{b,c}, var] (* nothing did apply, handle argument of exponential function *) UserInt1[Exp[expr_] term_., var_] := UserInt1[Exp[(expr // Expand) //. $CollectVar[var]] term, var] /; !FreeQ[expr, var] End[] Protect[GaussInt] EndPackage[]