Double integral over the region enclosed by several lines

Question

The question is:

Calculate the double integral $$\iint_D (18x^2+51xy+30y^2)^2\,dx\,dy$$ where $D$ is an area enclosed by the following lines: $$6x+5y=9,\quad 6x+5y=-9,\quad 3x+6y=7, \quad 3x+6y=-7$$

I have tried doing it the following way but it's wrong and I can't seem to get it working.

NIntegrate[
 Boole[-9 <= 6 x + 5 y <= 9 && -7 <= 3 x + 6 y <= 7], 
 {x, -Infinity, Infinity}, {y, -Infinity, Infinity}
]

Answer 1

Since Mathematica 10 there's been ImplicitRegion. We define that region this way:

IR = ImplicitRegion[-9 <= 6 x + 5 y <= 9 && -7 <= 3 x + 6 y <= 7, {x, y}];

then we can calculate this integral directly:

Integrate[(18 x^2 + 51 x*y + 30 y^2)^2, {x, y} ∈ IR]

5292

This integral can be calculated without the newest version of the system, to enlighten its structure we proceed in the following way:

18 x^2 + 51 x*y + 30 y^2 // Factor

3 (x + 2 y) (6 x + 5 y)

We change linearly variables: w == 3x + 6y; z == 6x + 5y;. The Jacobian of this transformation we calculate directly:

Det @ D[{3 x + 6 y, 6 x + 5 y}, {{x, y}}]

-21

We take the absolute value of its inverse and the integral can be factorized, thus :

int = 1/21 Integrate[w^2 z^2, {w, -9, 9}, {z, -7, 7}]

 5292

For completness we provide another two straightforward ways which work also in earlier versions of the system:

Integrate[ 
  Boole[-9 <= 6 x + 5 y <= 9 && -7 <= 3 x + 6 y <= 7] (18 x^2 + 51 x y + 30 y^2)^2,
  {x, -∞, ∞}, {y, -∞, ∞}]

Integrate[ 
  HeavisideTheta[6 x + 5 y + 9] HeavisideTheta[9 - 6 x - 5 y] 
  HeavisideTheta[3 x + 6 y + 7] HeavisideTheta[7 - 3 x - 6 y] 
   (18 x^2 + 51 x y + 30 y^2)^2, {x, -∞, ∞}, {y, -∞, ∞}]

Answer 2

This post is cosmetic. It adds nothing to Artes excellent and instructive answer (which I have voted for). I post it to illustrate visualization of the region (a parallelogram in this instance) and use same change of variables as Artes. I think SliceContourPlot3D is also a a nice way to visualize.

m = {{6, 5}, {3, 6}};
mi = Inverse[m];
mi3 = ArrayFlatten[{{mi, 0}, {0, 1}}];
f[x_, y_] := Times @@ (m.{x, y})^2;
j = Abs[Det@D[m.{x, y}, {{x, y}}]];
cp = ContourPlot[
   Evaluate@Table[m.{x, y} == k {9, 7}, {k, {1, -1}}], {x, -5, 
    5}, {y, -5, 5}, PlotLegends -> "Expressions"];
pp = ParametricPlot[mi.{u, v}, {u, -9, 9}, {v, -7, 7}, 
   PlotStyle -> Yellow];
sh1 = Legended[Show[{cp, pp}], SwatchLegend[{Yellow}, {"Region D"}]];
reg = ParametricRegion[mi.{u, v}, {{u, -9, 9}, {v, -7, 7}}];
i1 = HoldForm[Integrate[a^2 b^2, {a, -9, 9}, {b, -7, 7}]/j];
i2 = HoldForm[Integrate[f[x, y], {x, y} \[Element] reg]];
p3 = Plot3D[f[x, y], {x, y} \[Element] reg, Mesh -> None, 
   PlotStyle -> Opacity[0.2]];
reg3 = ParametricRegion[
   mi3.{u, v, w}, {{u, -9, 9}, {v, -7, 7}, {w, 0, 1}}];
scp = SliceContourPlot3D[
   z - f[x, y], {z == 0}, {x, y, z} \[Element] reg3, 
   PlotPoints -> 100, Contours -> {0., -10, -100, -1000}];
sh2 = Show[p3, scp, ImageSize -> 300];
vr = ParametricRegion[
   mi3.{u, v, w}, {{u, -9, 9}, {v, -7, 7}, {w, 0, f @@ (mi.{u, v})}}];
vol = Volume[vr];
TraditionalForm@
 Grid[{{Row[{HoldForm@f[x, y], " = ", f[x, y], " = ", 
      Expand[Times @@ (m.{x, y})]^2}], SpanFromLeft}, {sh1, sh2}, 
   Row[{#, " = ", ReleaseHold@#}, BaseStyle -> {20}] & /@ {i1, 
     i2}, {Row[{"Parametric region D: ", vr[[1]]}], 
    Row[{"Volume of D = ", vol}]}}, Alignment -> Left, Frame -> All]