1
$\begingroup$

I have a very ugly function with 157 terms that I want to integrate over three variables. Each term is a Gaussian integral in those variables multiplied by a polynomial in the variables. However, in each term the argument of the exponential is different (still Gaussian) and each polynomial is different. I want to I will give an example term at the bottom (I've also attached a picture in case this is too tough to parse).

Is there a way to extract the coefficient in front of the exponent in such an expression, and to extract the argument of the exponent? I'm asking because I saw a post that deals with improving the speed of these types of integrals. That post provides a module to do the integration given the exponent's argument, the coefficient expression, and the variables. So if I can extract the relevant argument and exponent, I should in principle be able to use that module.

(I E^((3 I m \[Omega] (-2 q q1 + (q^2 + q1^2) Cos[t1 \[Omega]]) Csc[
    t1 \[Omega]])/(2 \[HBar]) + (
  I m \[Omega] (-2 q1 q2 + (q1^2 + 
        q2^2) Cos[(-t1 + t2) \[Omega]]) Csc[(-t1 + t2) \[Omega]])/(
  2 \[HBar]) + (
  I m \[Omega] (-2 q q2 + (q^2 + 
        q2^2) Cos[(-t2 + t3 - I \[Beta]) \[Omega]]) Csc[(-t2 + t3 - 
       I \[Beta]) \[Omega]])/(2 \[HBar])) m q1^5 q2 \[Lambda] Cot[
  t1 \[Omega]] Csc[t1 \[Omega]]^3 Sqrt[-((
  I m \[Omega] Csc[t1 \[Omega]])/\[HBar])] Sqrt[-((
  I m \[Omega] Csc[(-t1 + t2) \[Omega]])/\[HBar])] Sqrt[-((
  I m \[Omega] Csc[(-t2 + t3 - 
       I \[Beta]) \[Omega]])/\[HBar])])/(64 Sqrt[2] \[Pi]^(
 5/2) \[HBar])

expression

$\endgroup$
  • $\begingroup$ I suspect that your integration variables are q1, q2 and q. Is that correct? $\endgroup$ – Jens Jul 19 '17 at 0:00
1
$\begingroup$

I think you're referring to my answer in the linked thread. Using that, you could do the following:

gaussMoment[fPre_, fExp_, vars_] := 
 Module[{coeff, dist, ai, μ, norm}, 
  coeff = CoefficientArrays[fExp, vars, "Symmetric" -> True];
  ai = Inverse[2 coeff[[3]]];
  μ = -ai.coeff[[2]];
  dist = MultinormalDistribution[μ, -ai];
  norm = 1/PDF[dist, vars] /. Thread[vars -> μ];
  Simplify[
   norm Exp[1/2 coeff[[2]].μ + coeff[[1]]] Distribute@
     Expectation[fPre, vars \[Distributed] dist]]]


expr = (I E^((3 I m ω (-2 q q1 + (q^2 + q1^2) Cos[
              t1 ω]) Csc[t1 ω])/(2 ℏ) + (I m ω (-2 q1 q2 + \
(q1^2 + q2^2) Cos[(-t1 + t2) ω]) Csc[(-t1 + 
              t2) ω])/(2 ℏ) + (I m ω (-2 q q2 + \
(q^2 + q2^2) Cos[(-t2 + t3 - I β) ω]) Csc[(-t2 + t3 - 
              I β) ω])/(2 ℏ)) m q1^5 q2 λ \
Cot[t1 ω] Csc[t1 ω]^3 Sqrt[-((I m ω Csc[
           t1 ω])/ℏ)] Sqrt[-((I m ω Csc[(-t1 + 
              t2) ω])/ℏ)] Sqrt[-((I m ω Csc[(-t2 +
               t3 - I β) ω])/ℏ)])/(64 Sqrt[2] Pi^(5/2) ℏ);

expr /. Times[f___, Power[E, g_]] :> gaussMoment[Times[f], g, {q, q2, q1}]

$$\frac{15 \lambda \hbar ^2 \cot (\text{t1} \omega ) \csc ^3(\text{t1} \omega ) \sqrt{-\frac{i m \omega \csc (\text{t1} \omega )}{\hbar }} (\sin (\omega (i \beta +2 \text{t1}-\text{t3}))-2 \sin (\omega (\text{t3}-i \beta )))^2 \sqrt{\frac{i m \omega \csc (\omega (\text{t1}-\text{t2}))}{\hbar }} (-2 \sin (\omega (-i \beta +\text{t1}-\text{t2}+\text{t3}))+\sin (\omega (i \beta +\text{t1}+\text{t2}-\text{t3}))+3 \sin (\omega (\text{t1}-\text{t2}))) \sqrt{\frac{i m \omega \csc (\omega (i \beta +\text{t2}-\text{t3}))}{\hbar }} \sqrt{\frac{i \hbar ^3 \sin (\text{t1} \omega ) \sin (\omega (\text{t2}-\text{t1})) \sin (\omega (-i \beta -\text{t2}+\text{t3}))}{m^3 \omega ^3 (\cos (\omega (i \beta +2 \text{t1}-\text{t3}))-4 \cos (\omega (\text{t3}-i \beta ))+3)}}}{256 \sqrt{2} \pi m^2 \omega ^3 (\cos (\omega (i \beta +2 \text{t1}-\text{t3}))-4 \cos (\omega (\text{t3}-i \beta ))+3)^3}$$

What I did here is to replace the combination of E^(something) times any prefactor by the Gaussian integral, as calculated with the function gaussMoment from the answer linked here. This assumes that the expression is simply a sum of such exponentials, so that the integration over q, q1 and q2 can be distributed over it. That's what the replacement rule essentially does.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.