Below is a problem in an degree exam in real analysis: $$\int_{\Bbb R^2}e^{-(4x^2+4xy+5y^2)}d(x,y)=?$$
The original idea by me is to diagonalize such quadratic form. The result of diagonalization of the quadratic form $4x^2+4xy+5y^2$ is: 
However, the eigenvalues are very ugly. How to do the rest step to compute the integral? Is there an easy way?
Using the substitution I hinted at in a comment to this answer: $u=2x+y$ and $v=2y$, we get $$ \begin{align} \int_{\mathbb{R}}\int_{\mathbb{R}}e^{-4x^2-4xy-5y^2}\,\mathrm{d}x\,\mathrm{d}y &=\frac14\int_{\mathbb{R}}\int_{\mathbb{R}}e^{-u^2-v^2}\,\mathrm{d}u\,\mathrm{d}v\\ &=\frac14\int_0^{2\pi}\int_0^\infty e^{-r^2}r\,\mathrm{d}r\,\mathrm{d}\theta\\ &=\frac14\cdot2\pi\cdot\frac12 \end{align} $$ The polar substitution is $u=r\cos(\theta)$ and $v=r\sin(\theta)$.
Hopefully the test mentioned in your comment to that answer is an old test and not one you are currently taking.
You do not need to actually perform a diagonalization. If $q(x,y)=a x^2+2bxy+cy^2$ is associated to a positive definite matrix $Q=\begin{pmatrix}a & b \\ b & c \end{pmatrix}$, then $Q=J^{-1} D J$ where $J^{-1}=J^T$ and $D$ is a diagonal matrix containing the eigenvalues $\lambda_1,\lambda_2>0$ of $Q$. By performing the substitution $J(x,y)^T = (X,Y)^T$ we have
$$\begin{eqnarray*} \iint_{\mathbb{R}^2}e^{-q(x,y)}\,dx\,dy &=&\iint_{\mathbb{R}^2}e^{-(x,y)Q(x,y)^T}\,dx\,dy\\&=&\iint_{\mathbb{R}^2}e^{-(x,y)J^{-1} D J(x,y)^T}\,dx\,dy\\&=&\iint_{\mathbb{R}^2}e^{-(X,Y)Q(X,Y)^T}\,\left|\det J^{-1}\right|dX\,dY\\&=&\iint_{\mathbb{R}^2}e^{-(\lambda_1 X^2+\lambda_2 Y^2)}\,dX\,dY\\\small{\left(X=\tfrac{u}{\sqrt{\lambda_1}},Y=\tfrac{v}{\sqrt{\lambda_2}}\right)}\quad&=&\frac{1}{\sqrt{\lambda_1 \lambda_2}}\iint_{\mathbb{R}^2}e^{-u^2-v^2}\,du\,dv\\(\text{Fubini})\quad&=&\frac{1}{\sqrt{\det Q}}\left(\int_{\mathbb{R}}e^{-u^2}\,du\right)^2\\&=&\frac{\pi}{\sqrt{\det Q}}=\color{red}{\frac{\pi}{\sqrt{ac-b^2}}}\end{eqnarray*}$$ and a similar identity holds for $\iiint_{\mathbb{R}^3}e^{-q(x,y,z)}\,dx\,dy\,dz$, for instance. The elements of $J$ or the exact values of $\lambda_j$ do not really matter, just the positive definiteness is needed to ensure convergence.
Your quadratic form is given as:
$$\boldsymbol{x}^T\boldsymbol{Px}$$
Now use the substitution (for a justification see below)
$$\boldsymbol{x} = \boldsymbol{V}\boldsymbol{\Lambda}^{-1/2} \boldsymbol{w} ,$$
in which $\boldsymbol{\Lambda}$ is a diagonal matrix containing the eigenvalues of $\boldsymbol{P}$ and $\boldsymbol{V}$ is the matrix containing the normalized eigenvectors associated with $\boldsymbol{\Lambda}$. Note, that in your case you can choose the eigenvectors in such a fashion that $\boldsymbol{V}$ is orthonormal. Then do the substitution of the integral by using the determinant of the Jacobian for this substitution.
Edit: In order to solve the problem, you will need to determine the Jacobian of the substitution
$$\boldsymbol{x} = \boldsymbol{V}\boldsymbol{\Lambda}^{-1/2} \boldsymbol{w}$$
as this is a linear expression the Jacobi determinant is given by
$$\det \left[\boldsymbol{V}\boldsymbol{\Lambda}^{-1/2}\right].$$
The integral is then given by
$$\int_{\Bbb R^2}e^{-\boldsymbol{w}^T\boldsymbol{w}}\det \left[\boldsymbol{V}\boldsymbol{\Lambda}^{-1/2}\right]d\boldsymbol{w}=\det \left[\boldsymbol{V}\boldsymbol{\Lambda}^{-1/2}\right]\int_{\Bbb R^2}e^{-\boldsymbol{w}^T\boldsymbol{w}}d\boldsymbol{w}.$$
The last expression is a well-known result obtained by Laplace by applying a more general formula of Euler (it is often not correctly referred as the Gauss integral). It can be solved by the trigonometric substitution
$$\boldsymbol{w} = \begin{bmatrix}r\cos \varphi\\ r \sin \varphi\end{bmatrix} \implies d\boldsymbol{w} = rdrd\varphi.$$
In order to motivate my answer, we will do this step by step.
From the eigenvalue equation in matrix form and the orthonormality of $\boldsymbol{V}$
$$\boldsymbol{PV} = \boldsymbol{V\Lambda}$$
we can obtain
$$\boldsymbol{\Lambda} = \boldsymbol{V}^{-1}\boldsymbol{PV}=\boldsymbol{V}^T\boldsymbol{PV}.$$
If we transform our quadratic form by the substitution $$\boldsymbol{x} = \boldsymbol{Vz}$$
we obtain
$$\boldsymbol{z}^T\boldsymbol{V}^T\boldsymbol{PVz}=\boldsymbol{z}^T\boldsymbol{\Lambda z}.$$
Now to transform the diagonal matrix $\boldsymbol{\Lambda}$ into the identity matrix we use
$$\boldsymbol{z}=\boldsymbol{\Lambda}^{-1/2}\boldsymbol{w}$$
if we apply this to the transformed quadratic form we can obtain
$$\boldsymbol{z}^T\boldsymbol{V}^T\boldsymbol{PVz}=\boldsymbol{w}^T\boldsymbol{\Lambda}^{-T/2}\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{-1/2}\boldsymbol{w}.$$
As $\boldsymbol{\Lambda}$ is a diagonal matrix
$$\boldsymbol{\Lambda}^{-T/2}=\boldsymbol{\Lambda}^{-1/2}.$$
$$\boldsymbol{w}^T\boldsymbol{\Lambda}^{-T/2}\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{-1/2}\boldsymbol{w}=\boldsymbol{w}^T\boldsymbol{\Lambda}^{-1/2}\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{-1/2}\boldsymbol{w}=\boldsymbol{w}^T\boldsymbol{\Lambda}^{-1/2}\boldsymbol{\Lambda}^{1/2}\boldsymbol{\Lambda}^{1/2}\boldsymbol{\Lambda}^{-1/2}\boldsymbol{w}=\boldsymbol{w}^T\boldsymbol{w}.$$
So we really obtain a simple quadratic form. Chaining both substitutions gives
$$\boldsymbol{x} = \boldsymbol{V\Lambda}^{-1/2}\boldsymbol{w}$$
as was proposed previously.
Let $A = \pmatrix{2 & 1 \\0 & 2}$, so that $A(x, y) = (2x + y, 2y)$. Then making the substitution $(u, v) = A(x, y)$, \begin{align*} \int_{\mathbb{R}^2}\, dx\, dy\; e^{-(4x^2 + 4xy + 5y^2)} = \frac{1}{\det A}\int_{\mathbb{R}^2}\, du\, dv\; e^{-(u^2 + v^2)}, \end{align*} which should be clear.
