$$I = \iint_D dx dy, \quad D: \frac{x^2}{a^2} +\frac{y^2}{b^2} \leq 1$$
\(\displaystyle I = \iint_D dx dy\) は領域 \(D\) の面積を表すのであった。 \(\displaystyle \frac{x^2}{a^2} +\frac{y^2}{b^2} \leq 1\) つまり，領域 \(D\) は長半径・短半径 \(a, b\) の楕円の面積であるから，公式より \(I = \pi a b\) である。では実際に \(\displaystyle \iint_{D} dx dy\) をどうやって計算するかという話。
「円の面積を2重積分で求める」で，円の面積は求めているので，参考までに楕円の面積についてもまとめておく。楕円の面積はケプラーの第2法則（面積速度一定則）の際に使うので，念のために Maxima で求めている（「Maxima-Jupyter で楕円の面積を求める」）が，Maxima に頼らず，人力でも積分したくなるでしょ？
以下で示すように，素直に楕円の中心を原点としたデカルト座標で累次積分してください。間違って，楕円の焦点を原点とした極座標を使って2重積分しようとするものなら，痛い目にあいます。
累次積分を使う
領域 \(D\) の条件式から，
$$D: \frac{x^2}{a^2} +\frac{y^2}{b^2} \leq 1 \quad \Rightarrow \quad y^2 \leq b^2 \left(1 – \frac{x^2}{a^2} \right)$$
$$\therefore -b\sqrt{1 – \frac{x^2}{a^2}} \leq y \leq b \sqrt{1 – \frac{x^2}{a^2}} $$
\begin{eqnarray}
I &=& \iint_{D} dx dy \\
&=& \int_{-a}^{a} \left\{\int_{-b\sqrt{1 – \frac{x^2}{a^2}}}^{b \sqrt{1 – \frac{x^2}{a^2}}} dy\right\} dx\\
&=& 2 \int_{-a}^a b \sqrt{1 – \frac{x^2}{a^2}} dx
\end{eqnarray}
高校数学では，$y = f(x)$ と $y = g(x)$ が区間 $x_1 \leq x \leq x_2$ で $f(x) \geq g(x)$ のとき，$y = f(x), \ y = g(x), \ x = x_1, \ x = x_2$ で囲まれる部分の面積 $S$ が
$$S = \int_{x_1}^{x_2} \left( f(x) – g(x) \right) \, dx$$
となると習ったと思うが，この式が上記の累次積分で導かれたことになる。
ここで，\(x = a \sin\theta\) と変数変換すると，\(\displaystyle \sqrt{1 – \frac{x^2}{a^2}} = \cos\theta, dx = a \cos\theta d\theta\)。
\begin{eqnarray}
I &=& 2 \int_{-a}^a b \sqrt{1 – \frac{x^2}{a^2}} dx \\
&=& 2 a b \int_{-\pi/2}^{\pi/2} \cos^2\theta d\theta\\
&=& a b \int_{-\pi/2}^{\pi/2} (1 + \cos 2\theta) d\theta \\
&=&a b \left[ \theta + \frac{1}{2}\sin 2\theta\right]_{-\pi/2}^{\pi/2} = \pi a b
\end{eqnarray}
$$\therefore I = \pi a b$$
焦点を原点とした極座標による累次積分
焦点を原点とした極座標 $r, \phi$ で表した楕円の式は
$$r(\phi) = \frac{a (1-e^2)}{1 + e \cos\phi}, \quad b = a \sqrt{1-e^2}$$
であった。この極座標で表すと，領域 \(D\) は
$$ 0 \leq r \leq r(\phi) , \quad 0 \leq \phi \leq 2 \pi$$
となる。したがって楕円の面積は以下のようになるが…
\begin{eqnarray} I &=& \iint_{D} r \,dr \,d\phi \\
&=& \int_0^{2\pi} d\phi \int_0^{r(\phi)} r \, dr \\
&=& \frac{1}{2} \int_0^{2\pi} r^2(\phi) \,d\phi \\
&=& a^2 (1-e^2)^2 \int_0^{2\pi} \frac{1}{2 (1+e\cos\phi)^2}\, d\phi
\end{eqnarray}
$\cos\phi$ の有理関数 $\displaystyle \frac{1}{(1+e\cos\phi)^2}$ の積分になり，ハマります。以前やってた天文学の授業で，「天文学では楕円の中心を原点としたデカルト座標系など使わない。焦点を原点とした極座標だけが出てくるのだ。」などと力説した手前，楕円の面積も $r(\phi)$ だけで求めようとして，大ハマりしました。これはやってはいけない例として。
答えは既に $ I = \pi a b$ と求まっているので，こんな煩わしい積分に時間を潰したくないと思う人は，以下は読み飛ばしてください。
まず，授業でやった「sin 𝑥, cos 𝑥 の有理関数の積分」のセオリーにそって，以下のような変換を行なって置換積分の格好にする。
\begin{eqnarray}
\tan\frac{\phi}{2}&\equiv& t\\
d\phi &=& \frac{2}{1 + t^2} dt \\
\cos\phi &=& \frac{1-t^2}{1+t^2}
\end{eqnarray}
まずは不定積分の形で，変数変換後，部分分数に分解し，さらに部分積分も駆使して…（ってこんなの私の腕力ではできません。ここだけこっそり Maxima 使いました。）
\begin{eqnarray}&&
\int \frac{1}{2 (1+e\cos\phi)^2}\, d\phi \\
&=& \int \frac{1}{\left(1+e\frac{1-t^2}{1+t^2} \right)^2}\frac{dt}{1+t^2}\\
&=& \int \frac{1+t^2}{\left\{(1-e) t^2 + (1+e) \right\}^2}\,dt \\
&=& \int \left\{\frac{1+e}{ (1-e^2)^{-\frac{3}{2}}} \frac{\sqrt{\frac{1-e}{1+e}}}{1 + \frac{1-e}{1+e} t^2}
– \frac{2 e}{ (1-e^2)^{-\frac{3}{2}}}\frac{\sqrt{\frac{1-e}{1+e}}}{\left(1 + \frac{1-e}{1+e} t^2\right)^2}\right\} dt \\
&&\qquad\qquad \left( x \equiv \sqrt{\frac{1-e}{1+e}} t \right)\\&=& \frac{1+e}{ (1-e^2)^{-\frac{3}{2}}}\int \frac{dx}{1 + x^2}
– \frac{e}{ (1-e^2)^{-\frac{3}{2}}}\int \frac{2 dx}{\left(1 + x^2\right)^2} \\
&=& \frac{1}{ (1-e^2)^{-\frac{3}{2}}}\int \frac{dx}{1 + x^2} – \frac{e}{ (1-e^2)^{-\frac{3}{2}}}\frac{ x}{1+x^2} \\
&=& (1-e^2)^{-\frac{3}{2}} \tan^{-1} \left(\sqrt{\frac{1-e}{1+e}} t \right)
– \frac{e t}{(1-e^2) \left\{ (1-e)t^2 + (1+e)\right\}}
\end{eqnarray}
(ここで $ \int \frac{2 dx}{\left(1 + x^2\right)^2} = \int \frac{dx}{1+x^2} + \frac{x}{1+x^2}$ を使った。)
上記第2項は，$t = \tan\frac{\phi}{2} $ が $0$ のときも $\pm\infty$ のときもゼロになるので定積分には寄与しない。したがって（積分区間内で $\tan\frac{\phi}{2}$ の発散があるから慎重に）
\begin{eqnarray}&&
\int_0^{2\pi} \frac{1}{2 (1+e\cos\phi)^2}\, d\phi \\
&=& (1-e^2)^{-\frac{3}{2}} \Biggl[\tan^{-1} \left(\sqrt{\frac{1-e}{1+e} }\tan\frac{\phi}{2} \right)
\Biggr]_0^{\pi-0}\\
&& + (1-e^2)^{-\frac{3}{2}} \Biggl[\tan^{-1} \left(\sqrt{\frac{1-e}{1+e} }\tan\frac{\phi}{2} \right)
\Biggr]_{\pi+0}^{2\pi} \\
&=& (1-e^2)^{-\frac{3}{2}} \left\{ \left(\frac{\pi}{2} – 0\right) + \left(0 – \left(-\frac{\pi}{2}\right) \right)\right\}\\
&=& (1-e^2)^{-\frac{3}{2}} \,\pi
\end{eqnarray}
最終的に
\begin{eqnarray}
I &=& a^2 (1-e^2)^2 \int_0^{2\pi} \frac{1}{2 (1+e\cos\phi)^2}\, d\phi\\
&=& a^2 (1-e^2)^2 \times (1-e^2)^{-\frac{3}{2}} \,\pi \\
&=& \pi a^2 \sqrt{1-e^2} \\
&=& \pi a b
\end{eqnarray}
追記：離心近点離角を使った置換積分
Memo「真近点離角と離心近点離角との関係についてもう少し」にまとめたように，以下のように変数 $u$ を使った置換積分にすると，簡単。$u$ は業界用語で「離心近点離角」または「離心近点角 (Wikipedia)」と呼ぶ。ここでは木下宙著「天体と軌道の力学」にならい，「離心近点離角」で。
離心近点離角の幾何学的意味づけも興味深いが，ここではあくまで置換積分のための変数変換であるという数学的道具としての有効性のみを強調するのみにとどめる。角度座標 $\phi$ （これは業界用語で「真近点離角」とも呼ばれる）と $u$ の関係は，
\begin{eqnarray}
\frac{1 -e^2}{1 + e \cos\phi} &\equiv& 1 -e \cos u \\
d\phi &=& \frac{\sqrt{1 -e^2}}{1 -e \cos u} \,du
\end{eqnarray}
$u$ を使って置換積分すると，（積分範囲が $[0, 2\pi]$ のままなのも便利なところ）
\begin{eqnarray}
\int_0^{2\pi} \frac{1}{2 (1+e\cos\phi)^2}\, d\phi
&=& \frac{1}{2 (1 -e^2)^2} \int_0^{2 \pi} (1 -e \cos u)^2 \cdot \frac{\sqrt{1 -e^2}}{1 -e \cos u} \,du \\
&=& \frac{\sqrt{1 -e^2}}{2 (1 -e^2)^2} \int_0^{2 \pi} (1 -e \cos u)\, du \\
&=& \frac{\sqrt{1 -e^2}}{2 (1 -e^2)^2} \Bigl[ u -e \sin u \Bigr]_0^{2 \pi} \\
&=& \frac{ \pi \sqrt{1 -e^2}}{(1 -e^2)^2}
\end{eqnarray}