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} ] 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, -∞, ∞}] 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] 
Integrateinstead ofIntegrateto get exact answer; correct typo iny-domain so it runs from-Infinityat lower end! $\endgroup$Integrateinstead ofNIntegrate(typo was from Mac's "helpful" auto-correct). $\endgroup$