1
$\begingroup$

I have this complex function $wig2$ which I am trying to double integrate with respect to 2 variables $q$ and $p$. Using the code below, I have tried to run it several times; however, it takes forever to compute and doesn't produce a final answer. It just keeps on running. I have tried leaving it for hours but it didn't work.

I have also tried using NIntegrate by fixing limits of $q$ from $-5$ to $5$ but it says the integrand has evaluated to a non-numerical value at the sampling points.

Please guide me why this is happening and how can I fix this?

wig2 = (1/(81*Pi))*4*Abs[(97 - 648*I*Sqrt[2]*p - 2992*p^2 + 2624*I*Sqrt[2]*p^3 + 1088*p^4 - 648*Sqrt[2]*q + 5984*I*p*q + 7872*Sqrt[2]*p^2*q - 4352*I*p^3*q + 2992*q^2 - 7872*I*Sqrt[2]*p*q^2 - 6528*p^2*q^2 - 2624*Sqrt[2]*q^3 + 4352*I*p*q^3 + 1088*q^4 + 64*(17 + 64*p^4 + 64*I*p^3*(3*Sqrt[2] - 4*q) - 88*Sqrt[2]*q + 304*q^2 - 192*Sqrt[2]*q^3 + 64*q^4 + 16*p^2*(-19 + 36*Sqrt[2]*q - 24*q^2) - 8*I*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3))*Conjugate[p]^4 - 8*(81*Sqrt[2] + 704*Sqrt[2]*p^4 - 1000*q + 2144*Sqrt[2]*q^2 - 3520*q^3 + 704*Sqrt[2]*q^4 - 704*I*p^3*(-5 + 4*Sqrt[2]*q) - 32*p^2*(67*Sqrt[2] - 330*q + 132*Sqrt[2]*q^2) + 8*I*p*(-125 + 536*Sqrt[2]*q - 1320*q^2 + 352*Sqrt[2]*q^3))*Conjugate[q] + 16*(187 + 1216*p^4 + 128*I*p^3*(25*Sqrt[2] - 38*q) - 1072*Sqrt[2]*q + 4240*q^2 - 3200*Sqrt[2]*q^3 + 1216*q^4 + 16*p^2*(-265 + 600*Sqrt[2]*q - 456*q^2) - 16*I*p*(67*Sqrt[2] - 530*q + 600*Sqrt[2]*q^2 - 304*q^3))*Conjugate[q]^2 - 64*(41*Sqrt[2] + 192*Sqrt[2]*p^4 - 440*q + 800*Sqrt[2]*q^2 - 1088*q^3 + 192*Sqrt[2]*q^4 - 64*I*p^3*(-17 + 12*Sqrt[2]*q) - 32*p^2*(25*Sqrt[2] - 102*q + 36*Sqrt[2]*q^2) + 8*I*p*(-55 + 200*Sqrt[2]*q - 408*q^2 + 96*Sqrt[2]*q^3))*Conjugate[q]^3 + 64*(17 + 64*p^4 + 64*I*p^3*(3*Sqrt[2] - 4*q) - 88*Sqrt[2]*q + 304*q^2 - 192*Sqrt[2]*q^3 + 64*q^4 + 16*p^2*(-19 + 36*Sqrt[2]*q - 24*q^2) - 8*I*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3))*Conjugate[q]^4 + 64*Conjugate[p]^3*(-192*I*Sqrt[2]*p^4 - 64*p^3*(-17 + 12*Sqrt[2]*q) + 32*I*p^2*(25*Sqrt[2] - 102*q + 36*Sqrt[2]*q^2) + 8*p*(-55 + 200*Sqrt[2]*q - 408*q^2 + 96*Sqrt[2]*q^3) - I*(41*Sqrt[2] - 440*q + 800*Sqrt[2]*q^2 - 1088*q^3 + 192*Sqrt[2]*q^4) - 4*(-64*I*p^4 + 64*p^3*(3*Sqrt[2] - 4*q) - 16*I*p^2*(-19 + 36*Sqrt[2]*q - 24*q^2) - 8*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3) + I*(-17 + 88*Sqrt[2]*q - 304*q^2 + 192*Sqrt[2]*q^3 - 64*q^4))*Conjugate[q]) + 16*Conjugate[p]^2*(-187 - 1216*p^4 - 128*I*p^3*(25*Sqrt[2] - 38*q) + 1072*Sqrt[2]*q - 4240*q^2 + 3200*Sqrt[2]*q^3 - 1216*q^4 + p^2*(4240 - 9600*Sqrt[2]*q + 7296*q^2) + 16*I*p*(67*Sqrt[2] - 530*q + 600*Sqrt[2]*q^2 - 304*q^3) + 12*(41*Sqrt[2] + 192*Sqrt[2]*p^4 - 440*q + 800*Sqrt[2]*q^2 - 1088*q^3 + 192*Sqrt[2]*q^4 - 64*I*p^3*(-17 + 12*Sqrt[2]*q) - 32*p^2*(25*Sqrt[2] - 102*q + 36*Sqrt[2]*q^2) + 8*I*p*(-55 + 200*Sqrt[2]*q - 408*q^2 + 96*Sqrt[2]*q^3))*Conjugate[q] + 24*(-17 - 64*p^4 - 64*I*p^3*(3*Sqrt[2] - 4*q) + 88*Sqrt[2]*q - 304*q^2 + 192*Sqrt[2]*q^3 - 64*q^4 + p^2*(304 - 576*Sqrt[2]*q + 384*q^2) + 8*I*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3))*Conjugate[q]^2) + 8*Conjugate[p]*(704*I*Sqrt[2]*p^4 + 704*p^3*(-5 + 4*Sqrt[2]*q) - 32*I*p^2*(67*Sqrt[2] - 330*q + 132*Sqrt[2]*q^2) - 8*p*(-125 + 536*Sqrt[2]*q - 1320*q^2 + 352*Sqrt[2]*q^3) + I*(81*Sqrt[2] - 1000*q + 2144*Sqrt[2]*q^2 - 3520*q^3 + 704*Sqrt[2]*q^4) + 4*(-1216*I*p^4 + 128*p^3*(25*Sqrt[2] - 38*q) - 16*I*p^2*(-265 + 600*Sqrt[2]*q - 456*q^2) - 16*p*(67*Sqrt[2] - 530*q + 600*Sqrt[2]*q^2 - 304*q^3) + I*(-187 + 1072*Sqrt[2]*q - 4240*q^2 + 3200*Sqrt[2]*q^3 - 1216*q^4))*Conjugate[q] + 24*(192*I*Sqrt[2]*p^4 + 64*p^3*(-17 + 12*Sqrt[2]*q) - 32*I*p^2*(25*Sqrt[2] - 102*q + 36*Sqrt[2]*q^2) - 8*p*(-55 + 200*Sqrt[2]*q - 408*q^2 + 96*Sqrt[2]*q^3) + I*(41*Sqrt[2] - 440*q + 800*Sqrt[2]*q^2 - 1088*q^3 + 192*Sqrt[2]*q^4))*Conjugate[q]^2 + 32*(-64*I*p^4 + 64*p^3*(3*Sqrt[2] - 4*q) - 16*I*p^2*(-19 + 36*Sqrt[2]*q - 24*q^2) - 8*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3) + I*(-17 + 88*Sqrt[2]*q - 304*q^2 + 192*Sqrt[2]*q^3 - 64*q^4))*Conjugate[q]^3))/E^(2*Abs[-(1/Sqrt[2]) + I*p + q]^2)/ ((Sqrt[2] - 4*I*p - 4*q)^2*(Sqrt[2] + 4*I*Conjugate[p] - 4*Conjugate[q])^2)];

abs = ComplexExpand[Abs[wig2]]];
intofwig = Integrate[abs, q];
doubint = Integrate[intofwig, p];
delta = Simplify[doubint] - 1;
Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ You write Integrate[abs, q]; but where is abs defined? May be you forgot the code for that? $\endgroup$
    – Nasser
    Commented Aug 20, 2023 at 8:09
  • 1
    $\begingroup$ Better is use: Conjugate[q] == Abs[q]^2/q and Conjugate[p] == Abs[p]^2/p for analytical solution if closed-form exist !. I Try with: NIntegrate[wig2, {q, -5, 5}, {p, -5, 5}, Method -> "LocalAdaptive"] and give me: 0.6804 with no errors. $\endgroup$
    – Mariusz Iwaniuk
    Commented Aug 20, 2023 at 8:25
  • $\begingroup$ @Nasser I edited it now, I wrote it in the code but there was problem with it displaying here $\endgroup$
    – Anaya
    Commented Aug 20, 2023 at 8:35
  • 1
    $\begingroup$ NIntegrate[abs, {q, -5, 5}, {p, -5, 5}, Method -> "LocalAdaptive"] $\endgroup$
    – Daniel Huber
    Commented Aug 20, 2023 at 9:12
  • 1
    $\begingroup$ The wig2 expression is to complicated to computed analyticaly. Probably only way is NIntegrate . $\endgroup$
    – Mariusz Iwaniuk
    Commented Aug 20, 2023 at 9:38

1 Answer 1

Reset to default
4
$\begingroup$

By going to polar coordinates the integrand becomes a lot simpler and can potentially be studied in more detail.

First, define wig2 without the Abs function:

wig2[p_, q_] = (1/(81*Pi))*4*((97 - 648*I*Sqrt[2]*p - 2992*p^2 + 2624*I*Sqrt[2]*p^3 + 1088*p^4 - 648*Sqrt[2]*q + 5984*I*p*q + 7872*Sqrt[2]*p^2*q - 4352*I*p^3*q + 2992*q^2 - 7872*I*Sqrt[2]*p*q^2 - 6528*p^2*q^2 - 2624*Sqrt[2]*q^3 + 4352*I*p*q^3 + 1088*q^4 + 64*(17 + 64*p^4 + 64*I*p^3*(3*Sqrt[2] - 4*q) - 88*Sqrt[2]*q + 304*q^2 - 192*Sqrt[2]*q^3 + 64*q^4 + 16*p^2*(-19 + 36*Sqrt[2]*q - 24*q^2) - 8*I*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3))*Conjugate[p]^4 - 8*(81*Sqrt[2] + 704*Sqrt[2]*p^4 - 1000*q + 2144*Sqrt[2]*q^2 - 3520*q^3 + 704*Sqrt[2]*q^4 - 704*I*p^3*(-5 + 4*Sqrt[2]*q) - 32*p^2*(67*Sqrt[2] - 330*q + 132*Sqrt[2]*q^2) + 8*I*p*(-125 + 536*Sqrt[2]*q - 1320*q^2 + 352*Sqrt[2]*q^3))*Conjugate[q] + 16*(187 + 1216*p^4 + 128*I*p^3*(25*Sqrt[2] - 38*q) - 1072*Sqrt[2]*q + 4240*q^2 - 3200*Sqrt[2]*q^3 + 1216*q^4 + 16*p^2*(-265 + 600*Sqrt[2]*q - 456*q^2) - 16*I*p*(67*Sqrt[2] - 530*q + 600*Sqrt[2]*q^2 - 304*q^3))*Conjugate[q]^2 - 64*(41*Sqrt[2] + 192*Sqrt[2]*p^4 - 440*q + 800*Sqrt[2]*q^2 - 1088*q^3 + 192*Sqrt[2]*q^4 - 64*I*p^3*(-17 + 12*Sqrt[2]*q) - 32*p^2*(25*Sqrt[2] - 102*q + 36*Sqrt[2]*q^2) + 8*I*p*(-55 + 200*Sqrt[2]*q - 408*q^2 + 96*Sqrt[2]*q^3))*Conjugate[q]^3 + 64*(17 + 64*p^4 + 64*I*p^3*(3*Sqrt[2] - 4*q) - 88*Sqrt[2]*q + 304*q^2 - 192*Sqrt[2]*q^3 + 64*q^4 + 16*p^2*(-19 + 36*Sqrt[2]*q - 24*q^2) - 8*I*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3))*Conjugate[q]^4 + 64*Conjugate[p]^3*(-192*I*Sqrt[2]*p^4 - 64*p^3*(-17 + 12*Sqrt[2]*q) + 32*I*p^2*(25*Sqrt[2] - 102*q + 36*Sqrt[2]*q^2) + 8*p*(-55 + 200*Sqrt[2]*q - 408*q^2 + 96*Sqrt[2]*q^3) - I*(41*Sqrt[2] - 440*q + 800*Sqrt[2]*q^2 - 1088*q^3 + 192*Sqrt[2]*q^4) - 4*(-64*I*p^4 + 64*p^3*(3*Sqrt[2] - 4*q) - 16*I*p^2*(-19 + 36*Sqrt[2]*q - 24*q^2) - 8*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3) + I*(-17 + 88*Sqrt[2]*q - 304*q^2 + 192*Sqrt[2]*q^3 - 64*q^4))*Conjugate[q]) + 16*Conjugate[p]^2*(-187 - 1216*p^4 - 128*I*p^3*(25*Sqrt[2] - 38*q) + 1072*Sqrt[2]*q - 4240*q^2 + 3200*Sqrt[2]*q^3 - 1216*q^4 + p^2*(4240 - 9600*Sqrt[2]*q + 7296*q^2) + 16*I*p*(67*Sqrt[2] - 530*q + 600*Sqrt[2]*q^2 - 304*q^3) + 12*(41*Sqrt[2] + 192*Sqrt[2]*p^4 - 440*q + 800*Sqrt[2]*q^2 - 1088*q^3 + 192*Sqrt[2]*q^4 - 64*I*p^3*(-17 + 12*Sqrt[2]*q) - 32*p^2*(25*Sqrt[2] - 102*q + 36*Sqrt[2]*q^2) + 8*I*p*(-55 + 200*Sqrt[2]*q - 408*q^2 + 96*Sqrt[2]*q^3))*Conjugate[q] + 24*(-17 - 64*p^4 - 64*I*p^3*(3*Sqrt[2] - 4*q) + 88*Sqrt[2]*q - 304*q^2 + 192*Sqrt[2]*q^3 - 64*q^4 + p^2*(304 - 576*Sqrt[2]*q + 384*q^2) + 8*I*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3))*Conjugate[q]^2) + 8*Conjugate[p]*(704*I*Sqrt[2]*p^4 + 704*p^3*(-5 + 4*Sqrt[2]*q) - 32*I*p^2*(67*Sqrt[2] - 330*q + 132*Sqrt[2]*q^2) - 8*p*(-125 + 536*Sqrt[2]*q - 1320*q^2 + 352*Sqrt[2]*q^3) + I*(81*Sqrt[2] - 1000*q + 2144*Sqrt[2]*q^2 - 3520*q^3 + 704*Sqrt[2]*q^4) + 4*(-1216*I*p^4 + 128*p^3*(25*Sqrt[2] - 38*q) - 16*I*p^2*(-265 + 600*Sqrt[2]*q - 456*q^2) - 16*p*(67*Sqrt[2] - 530*q + 600*Sqrt[2]*q^2 - 304*q^3) + I*(-187 + 1072*Sqrt[2]*q - 4240*q^2 + 3200*Sqrt[2]*q^3 - 1216*q^4))*Conjugate[q] + 24*(192*I*Sqrt[2]*p^4 + 64*p^3*(-17 + 12*Sqrt[2]*q) - 32*I*p^2*(25*Sqrt[2] - 102*q + 36*Sqrt[2]*q^2) - 8*p*(-55 + 200*Sqrt[2]*q - 408*q^2 + 96*Sqrt[2]*q^3) + I*(41*Sqrt[2] - 440*q + 800*Sqrt[2]*q^2 - 1088*q^3 + 192*Sqrt[2]*q^4))*Conjugate[q]^2 + 32*(-64*I*p^4 + 64*p^3*(3*Sqrt[2] - 4*q) - 16*I*p^2*(-19 + 36*Sqrt[2]*q - 24*q^2) - 8*p*(11*Sqrt[2] - 76*q + 72*Sqrt[2]*q^2 - 32*q^3) + I*(-17 + 88*Sqrt[2]*q - 304*q^2 + 192*Sqrt[2]*q^3 - 64*q^4))*Conjugate[q]^3))/E^(2*Abs[-(1/Sqrt[2]) + I*p + q]^2)/((Sqrt[2] - 4*I*p - 4*q)^2*(Sqrt[2] + 4*I*Conjugate[p] - 4*Conjugate[q])^2));

Go to polar coordinates by setting $p=r\cos\phi$ and $q=r\sin\phi$ and simplifying:

f[r_, ϕ_] = wig2[r Cos[ϕ], r Sin[ϕ]] // ComplexExpand // FullSimplify
(*    E^(-1 - 2 r^2 + 2 Sqrt[2] r Sin[ϕ]) *
      (97 + 736 r^2 + 64 r^4 + 8 r (-34 r Cos[2 ϕ] - 5 Sqrt[2] (13 + 8 r^2) Sin[ϕ])) /
      (81 π)    *)

Integrate numerically (without forgetting the Jacobian):

NIntegrate[Abs[f[r, ϕ]] r, {r, 0, ∞}, {ϕ, 0, 2 π}, Method -> "LocalAdaptive"]
(*    0.6804    *)

confirming @MariuszIwaniuk's result.

Improve this answer
$\endgroup$
3
  • $\begingroup$ This worked! Thanks a lot! Just wanted to know if there is any way to compute this integral analytically? $\endgroup$
    – Anaya
    Commented Aug 21, 2023 at 7:32
  • 1
    $\begingroup$ @Anaya it's easy if you get rid of the Abs function, or use Abs[...]^2. But simply with Abs, as you have it, it requires quite a few case distinctions. I'll give it a try. $\endgroup$
    – Roman
    Commented Aug 21, 2023 at 12:44
  • $\begingroup$ The $Abs$ function is part of the formula that I am using to compute my answer, so it's not possible to remove it. Thanks a lot for your help $\endgroup$
    – Anaya
    Commented Aug 22, 2023 at 6:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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