# Principal value of a double integral numerically

I have several integrals of the form : where a,b,wc and T are constants and PV denotes the Cauchy principal value I tried to integrate one of this form using NIntegrate however when trying to compute the principal value. I tried:

HLS[nu_, phi_, t_] :=
t^2 (Sinc[(omega - phi)*t/2]^2)*(alpha*
Exp[-nu/wc]*
nu)*( (1/(Exp[nu/T] - 1)) /(phi + nu))

NIntegrate[
HLS[nu, phi, 5], {nu, -Infinity, Infinity}, {phi, 0, Infinity},
Method -> "PrincipalValue"]


But I get NIntegrate::pvdim: PrincipalValue can be used for one-dimensional integrals only.

The principal value should make the integral convergent, I have spent a few days on it by now. I don't usually use mathematica but python and there the function is also only supported for one dimensional integrals so I gather this is a difficult problem. However I'd appreciate any advice or work arounds.

• The integral under consideration is not any principal value, but an iterated integral, where the inner integral is taken as its ptrincipal value. – user64494 Mar 10 at 8:24
• Perhaps I did not phrase it correctly still, the principal value needs to be taken in the inner integral, do you have any clue as to how to code it ? – Gerardo Suarez Mar 10 at 8:37
• Evaluate parameters, e.g. t = 1; omega = 4; phi = 5; alpha = 2; wc = 1; T = 2;. Unfortunately, then Integrate[ t^2 (Sinc[(omega - phi)*t/2]^2)*(alpha*Exp[-nu/wc]* nu)*((1/(Exp[nu/T] - 1))/(phi + nu)), {nu, -Infinity, Infinity}, PrincipalValue -> True] produces an error "Integrate::idiv: Integral of (2 E^-nu nu Sinc[1/2]^2)/((-1+E^(nu/2)) (5+nu)) does not converge on {-[Infinity],[Infinity]}.". – user64494 Mar 10 at 8:42

## 1 Answer

If you wish to compute it numerically for a specific set of constants then you could possibly create a table of principal-values of the inner integral over a range of $$\phi$$. Note in the code below I specify the two singular points of 0 and $$\phi$$ in the integration limits as {0,s1,s2,infty} then run the inner integration as principal-value as $$\phi$$ ranges from 0.1 to 15 and note the sum does appear to be converging for at least $$\phi\geq 0.1$$ for delta phi=1/4. Not sure what it would look like if you start decreasing delta phi to 1/100 or 1/1000. However your notation was a bit unclear. Perhaps you can further develop the code below to your liking.

Update: I think to integrate it this way, would have to multiply the value of each inner principal-value integral by delta phi like a Riemann sum which I did not do in the code below.

(*
Let a=1, t=1, b=1, w_c=1, T=2
*)
myf2[v_, \[Phi]_] := (v Exp[-v])/((Exp[v/2] - 1) (\[Phi] - v))
sum = 0;
myDoubleSum = Table[
sum +=
NIntegrate[
Sinc[1 + thePhi]^2 myf2[v, thePhi], {v, 0, 0,
thePhi, \[Infinity]}, Method -> "PrincipalValue"];
{thePhi, sum},
{thePhi, 0.1, 15, 1/4}];
ListPlot[myDoubleSum, Joined -> True] • Thanks I'll check it out and get back to you if I have any doubts. Sorry about the notation I mainly wanted to see a way to do it, summing over the second integral seems like a great way I did not think about it, it looks like what I expected! – Gerardo Suarez Mar 10 at 17:26
• Not quite sure about the legitimacy of doing that. Something you can work on perhaps. – Dominic Mar 10 at 17:35