# How to correctly calculate an improper integral with Mathematica?

Bug introduced in 13.0 or earlier and persisting through 13.2.0 or later.

I try in version 13 on Windows 10

j = Integrate[Exp[a*x^2 + b*x*y + c*y^2], {x, -Infinity, Infinity}, {y, -Infinity,
Infinity}, Assumptions -> {a, b, c} ∈ Reals]


ConditionalExpression[(2 π)/(Sqrt[-4 a + b^2/c] Sqrt[-c]), 4 a c^2 < b^2 c]

The above result is not correct, as

j /. {a -> 1, b -> 3, c -> 1}


-((2 I π)/Sqrt[5])

and

Integrate[Exp[a*x^2 + b*x*y + c*y^2] /. {a -> 1, b -> 3, c -> 1}, {x,-Infinity, Infinity}, {y, -Infinity, Infinity}]


Integrate::idiv:Integral of E^(x^2+3 x y+y^2) does not converge on {-∞,∞}.

demonstrate.

Is it possible to correctly calculate that integral with Mathematica? I tried to use the definition. Because the integrand is positive, it is enough to calculate

Integrate[Exp[a*x^2 + b*x*y + c*y^2], {x, -d, d}, {y, -d, d},
Assumptions -> {a, b, c} ∈ Reals && d > 0]


exhausting the plane over the squares and then to find the limit as d->Infinity. However, the latest integral almost crashes my comp.

• Integrate[Exp[a*x^2 + b*x*y + c*y^2], {x, -d, d}, {y, -d, d}, Assumptions -> {a, b, c} \[Element] Reals && d > 0, GenerateConditions -> False] returns a one-dimensional integral. Mar 14, 2022 at 17:42
• The same issue with Integrate[Exp[a*x^2 + b*x*y + c*y^2], {x, -d, d}, {y, -d, d}, Assumptions -> {a, b, c} > 0 && d > 0]. Mar 14, 2022 at 17:53

This can be done with help of Mathematica as follows. Let us exhaust the plane by disks centered at the origin. Switching to the polar coordinates, one obtains

Integrate[r*Exp[a*x^2 + b*x*y + c*y^2] /. {x -> r*Cos[ϕ],y -> r*Sin[ϕ]}, {r, 0, R}]


(-1 + E^(R^2 (a Cos[ϕ]^2 + b Cos[ϕ] Sin[ϕ] + c Sin[ϕ]^2)))/(a + c + (a - c) Cos[2 ϕ] + b Sin[2 ϕ])

It is clear the finite limit of the function defined by the above expression as R->Infinity exists only if ForAll[ϕ, ϕ > -Pi && ϕ <= Pi, a Cos[ϕ]^2 + b Cos[ϕ] Sin[ϕ] + c Sin[ϕ]^2 < 0] is valid as the result of

Integrate[r*Exp[a*x^2 + b*x*y + c*y^2] /. {x -> r*Cos[ϕ], y -> r*Sin[ϕ]},
{r, 0,Infinity},Assumptions -> {a, b, c} ∈ Reals && ϕ > -Pi && ϕ <= Pi]


ConditionalExpression[-(1/( 2 (a Cos[ϕ]^2 + b Cos[ϕ] Sin[ϕ] + c Sin[ϕ]^2))), a Cos[ϕ]^2 + b Cos[ϕ] Sin[ϕ] + c Sin[ϕ]^2 < 0]

demonstrates.

Unfortunately, Mathematica fails with

Resolve[ForAll[ϕ, ϕ > -Pi && ϕ <= Pi,
a Cos[ϕ]^2 + b Cos[ϕ] Sin[ϕ] + c Sin[ϕ]^2 <  0], Reals]


, returning the input. However, Mathematica is able to derive the condition in such a way

Resolve[ForAll[y, y >= -1 && y <= 1,a y^2 + b *y*Sqrt[1 - y^2] + c (1 - y^2) < 0], Reals]


c < 0 && a < b^2/(4 c)

Now

Integrate[(-1)/(a + c + (a - c) Cos[2 ϕ] + b Sin[2 ϕ]), {ϕ, -Pi, Pi},
Assumptions -> c < 0 && a < b^2/(4 c) && a < 0 && b ∈ Reals]


results in

ConditionalExpression[ -Piecewise[{{0, (Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) - Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] < 1 && Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) + Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] < 1) || (Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) - Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] >= 1 && Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) + Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] >= 1)}, {(2*Pi)/Sqrt[-b^2 + 4*a*c], Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) + Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] < 1}}, (-2*Pi)/Sqrt[-b^2 + 4*a*c]], !(Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) - Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] >= 1 && Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) - Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] <= 1) && !(Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) + Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] >= 1 && Sqrt[Abs[-(a/(a - I*b - c)) - c/(a - I*b - c) + Sqrt[-b^2 + 4*a*c]/(a - I*b - c)]] <= 1)]

, thinking a few hours. In fact, this is $$\frac{2 \pi }{\sqrt{4 a c-b^2}}$$ as Simplify under the conditions c < 0 && a < b^2/(4 c) && a < 0 && b ∈ Reals] shows.