Find the area enclosed by $\sqrt{(x-2)^2+(y-3)^2} + 2\sqrt{(x-3)^2+(y-1)^2} = 4$

This is an incomplete answer I have only covered the case where an oval is "large" in the sense it contains both of its foci in reasonable details.

To simplify the analysis, let start with $(x_1,y_1) = (-1,0), (x_2, y_2) = (1,0)$ and define $r, s, u, v$ such that

$$ \begin{cases}r^2 = (x+1)^2 + y^2\\s^2 = (x-1)^2 + y^2\end{cases} \quad\text{ and }\quad\begin{cases}u = \frac12 (r-s)\\v = \frac12(r+s)\end{cases} $$ We will look at the area of half of an oval in the upper half plane given by

$$\mathscr{O} = \big\{\; (x,y) \in \mathbb{R}^2 : y \ge 0,\, r + a s \le k \;\big\}$$

WOLOG, we can assume $a \ge 1$. Geometrically, there are three 3 possibilities:

1. $k > 2 a$ - The oval is "large" in the sense it encloses both of its foci.
2. $k = 2 a$ - The oval turns into a teardrop with a cusp at the focus $(-1,0)$.
3. $2a > k > 2$ - The oval encloses one and only one of the focus $(1,0)$.

In the following analysis, we will concentrate on the $1^{st}$ case where $k > 2 a$.

Notice $u^2 + v^2 = \frac12 (r^2 + s^2) = x^2 + 1 + y^2$. On the upper half plane, we have $$ \begin{cases} x = \frac14 (r^2 - s^2) = uv\\ \\ y = \sqrt{u^2 + v^2 - u^2 v^2 -1 } = \sqrt{(1-u^2)(v^2-1)} \end{cases}$$

The area element becomes

$$dx \wedge dy = \frac{dx \wedge d y^2}{2y} = \frac{duv \wedge d(u^2 + v^2)}{2y} = \frac{(u dv + v du)\wedge(u du + v dv)}{\sqrt{(1-u^2)(v^2-1)}}\\ = \frac{(v^2 - u^2) du \wedge dv}{\sqrt{(1-u^2)(v^2-1)}} = \frac{( (2v^2-1)+(1-2u^2)) du\wedge dv}{2\sqrt{(1-u^2)(v^2-1)}} $$ Notice $$d u\sqrt{1-u^2} = \frac{1-2u^2}{\sqrt{1-u^2}}\quad\text{ and }\quad d v\sqrt{v^2-1} = \frac{2v^2-1}{\sqrt{v^2-1}} $$ We get $$dx \wedge dy = \frac12 d \left[\frac{u(1-u^2)dv - v(v^2-1)du}{\sqrt{(1-u^2)(v^2-1)}}\right]$$

When $k > 2a$, the half oval $\mathscr{O}$ get mapped to a quadrilateral $\mathscr{Q}$ in the $uv$-plane.

$$\mathscr{Q} = \big\{\; (u,v) \in \mathbb{R}^2 : -1 \le u \le 1; 1 \le v \le \alpha + \beta u \;\}$$ where $ \displaystyle \alpha = \frac{k}{a+1},\;\beta = \frac{a-1}{a+1}$. This gives us

$$2 \text{Area}(\mathscr{O}) = 2\int_{\mathscr{O}} dx \wedge dy = \int_{\partial\mathscr{Q}} \frac{u(1-u^2)dv - v(v^2-1)du}{\sqrt{(1-u^2)(v^2-1)}} $$ It is easy to check on $\partial\mathscr{Q}$, the horizontal and the two vertical edges contribute nothing.
As a result, the area of the full oval can be rewritten as

$$2 \text{Area}(\mathscr{O}) = \int_{-1}^1 \frac{(\alpha+\beta u)((\alpha+\beta u)^2-1) - \beta u(1-u^2)}{\sqrt{(1-u^2)((\alpha+\beta u)^2-1)}} du\tag{*1} $$ Abusing the notation, we will use $v$ as a shorthand for $\alpha + \beta u$ from now on.

Consider the quadratic equation $t = (\alpha+\beta t)(\alpha t + \beta)$. When $k > 2a$, one can verify above equation has two real roots $\eta$ and $\eta^{-1}$ for $t$. We will let $\eta$ be the unique root with $|\eta| < 1$. Define $\lambda, \mu$ and introduce a new variable $z$ such that

$$\begin{cases} \lambda &= \alpha + \beta \eta\\ \mu &= \alpha \eta + \beta \end{cases} \quad\text{ and }\quad u = \frac{z + \eta}{1 + \eta z}\;\;\iff\;\; z = \frac{u - \eta}{1 - \eta u} $$ One can show that $|\mu| < 1$ and we have $$ du = \frac{(1-\eta^2)}{(1+\eta z)^2}dz,\quad \begin{cases} u & = \frac{z+\eta}{1+\eta z }\\ 1 + u & = \frac{(1+\eta)(1+z)}{1+\eta z}\\ 1 - u & = \frac{(1-\eta)(1-z)}{1+\eta z} \end{cases} \quad\text{ and }\quad \begin{cases} v &= \frac{\lambda + \mu z}{1 + \eta z}\\ v + 1 &= \frac{(\lambda+1)(1 + \mu z)}{1+\eta z}\\ v - 1 &= \frac{(\lambda-1)(1 - \mu z)}{1+\eta z} \end{cases}$$ This implies

$$\begin{align} v(v^2 - 1) & = \frac{\lambda^2 - 1}{(1+\eta z)^3}(\lambda + \mu z)(1 - \mu^2z^2)\\ u(1-u^2) &= \frac{1 - \eta^2}{(1+\eta z)^3}( z + \eta )(1 - z^2)\\ \frac{du}{\sqrt{(1-u^2)(v^2-1)}} &= \sqrt{\frac{1-\eta^2}{\lambda^2-1}}\frac{dz}{\sqrt{(1-z^2)(1-\mu^2 z^2)}} \end{align}$$

Notice $\beta( 1 - \eta^2 ) = \mu - \lambda \eta = \mu ( 1 - \lambda^2)$, we can rewrite $(*1)$ as

$$2 \text{Area}(\mathscr{O}) = \sqrt{(1-\eta^2)(\lambda^2-1)} \int_{-1}^1 \frac{P(z) dz}{ (1+\eta z)^3 \sqrt{(1-z^2)(1-\mu^2 z^2)}}\tag{*2}$$

where $P(z) = (\lambda + \mu z)(1 -\mu^2 z^2) + \mu (z + \eta)(1-z^2)$.

By brute force, one can show that

$$\frac{P(z)}{(1+\eta z)^3} = A(\eta,\mu) + \left[ B(\eta,\mu) + C(\eta,\mu)\frac{\partial}{\partial \eta} + D(\eta,\mu)\frac{\partial^2}{\partial \eta^2}\right]\frac{1}{1+\eta z}$$

where $$\begin{cases} A(\eta,\mu) &= -\frac{\mu^3+\mu}{\eta^3},\\ B(\eta,\mu) &= \frac{\mu^4 + (\eta^4+1)\mu^2+\eta^4}{\eta^3\mu},\\ C(\eta,\mu) &= -\frac{{\mu}^{4}+\left( -2\,{\eta}^{4}+2\,{\eta}^{2}+1\right) \,{\mu}^{2}-2\,{\eta}^{4}}{{\eta}^{2}\,\mu}\\ D(\eta,\mu) &= \frac{{\mu}^{4}+\left( {\eta}^{4}-4\,{\eta}^{2}+1\right) \,{\mu}^{2}+{\eta}^{4}}{2\,\eta\,\mu} \end{cases}$$

Substitute this in $(*2)$, then in terms of following complete elliptic integrals of $1^{st}, 2^{nd}$ and $3^{rd}$ kind:

$$\begin{align} K(\mu) &= \int_0^1 \frac{dz}{\sqrt{(1-z^2)(1-\mu^2 z^2)}}\\ E(\mu) &= \int_0^1 \sqrt{\frac{1-\mu^2 z^2}{1-z^2}} dz\\ \Pi(\eta^2,\mu) &= \int_0^1 \frac{dz}{(1 - \eta^2 z^2)\sqrt{(1-z^2)(1-\mu^2 z^2)}} \end{align}$$

We have $$2\text{Area}(\mathscr{O}) = 2\sqrt{(1-\eta^2)(\lambda^2-1)} \left( A K(\mu) + \left[ B + C \frac{\partial}{\partial \eta} + D\frac{\partial^2}{\partial \eta^2}\right]\Pi(\eta^2,\mu) \right)$$ Please note that the partial derivatives of $\Pi(\eta^2,\mu)$ can be expressed in terms of the complete elliptic integrals itself.

$$\frac{\partial\Pi(n,\mu)}{\partial n} = \frac{1}{2n(\mu^2 - n)(n-1)}\left(n E(\mu) + (\mu^2 - n)K(\mu) + (n^2 - \mu^2)\Pi(n,\mu)\right)$$

So in principle, we can get rid of the partial derivatives above and express the area of the oval solely in terms of complete elliptic integrals directly.

Let us switch to the $3^{rd}$ case where $2 < k < 2a$. The oval now intersect the line segment joining the two foci and encloses one of the focus $(1,0)$.

Let $\displaystyle \sigma = \frac{1-\alpha}{\beta}$. On the $uv$-plane, the line $v = \alpha + \beta u$ intersect the line $v = 1$ at $( \sigma, 1 )$ before the line $u = -1$. The half oval $\mathscr{O}$ get mapped to a triangle $\mathscr{T}$:

$$\mathscr{T} = \big\{\; (u,v) \in \mathbb{R}^2 : u \le 1; 1 \le v \le \alpha+\beta u\;\big\}$$

Once again, we can express the area of half oval as a line integral over $\partial \mathscr{T}$ and the horizontal and vertical edges contributes nothing. As a result, the area of the full oval is given by

$$2\text{Area}(\mathscr{O}) = \int_\sigma^1 \frac{(\alpha+\beta u)((\alpha + \beta u)^2 - 1) - \beta u (1-u^2) }{\sqrt{(1-u^2)((\alpha + \beta u)^2-1)}}\tag{*1'} $$

Since the lower limit $\sigma$ now depends on $\alpha, \beta$. We need a different Mobius transform which map between $u$ and a new variable $\tilde{z}$ such that $$\begin{array}{rcr} u& &\tilde{z}\\ \frac{-1-\alpha}{\beta} & \longleftrightarrow & -\frac{1}{\tilde{\mu}}\\ -1 & \longleftrightarrow & \frac{1}{\tilde{\mu}}\\ \sigma = \frac{1-\alpha}{\beta} & \longleftrightarrow & -1\\ 1 & \longleftrightarrow & 1 \end{array}$$ to bring the integrand back to the canonical form. If one perform this substitution, one will obtain something very similar to $(*2)$ above. The derivation will be very complicated and I will stop here. However, it is sort of clear if one do the dirty work along this direction, ultimately we will be able to express the area of oval again in terms of complete elliptic integrals.

Update

Back to original problem. Notice the distance between $(2,3)$ and $(3,1)$ is $\sqrt{5}$. If we scale the oval by a factor $\frac{2}{\sqrt{5}}$, it will reduce to the case 3. we have covered before with following parameters:

$$( k, a ) = (\frac{8}{\sqrt{5}}, 2 )\quad\implies\quad(\alpha,\beta,\sigma) = ( \frac{8}{3\sqrt{5}}, \frac13, 3 -\frac{8}{\sqrt{5}})$$

The area of the original oval is given by the integral in $(*1')$ scaled by a factor $\left(\frac{\sqrt{5}}{2}\right)^2$:

$$\text{Area} = \frac{5}{4} \int_{3 -\frac{8}{\sqrt{5}}}^1 \frac{\left( {\left( \frac{u}{3}+\frac{8}{3\,\sqrt{5}}\right) }^{2}-1\right) \,\left( \frac{u}{3}+\frac{8}{3\,\sqrt{5}}\right) -\frac{u\,\left( 1-{u}^{2}\right) }{3}}{\sqrt{\left( {\left( \frac{u}{3}+\frac{8}{3\,\sqrt{5}}\right) }^{2}-1\right) \,\left( 1-{u}^{2}\right) }} du$$

Throwing this expression to WA gives us $$3.2040508691540993294644483966405124836692674716427904\cdots$$

Minor progress: The integral in question is $$\int_{\gamma} y dx$$ where $\gamma$ runs around the Cartesian oval. To see this, break the integral up into $\gamma_1$ and $\gamma_2$ where $\gamma_1$ runs over the top of the oval and $\gamma_2$ over the bottom, both oriented from left to right. Then $$\int_{\gamma} y dx = \int_{\gamma_1} y dx - \int_{\gamma_2} y dx$$ and the second integral is the standard formula for the area between two curves.

The curve $\gamma$ is one of two connected components of the plane curve $$-711 - 184 x + 402 x^2 - 120 x^3 + 9 x^4 + 340 y + 80 x y - 12 x^2 y + 6 y^2 - 120 x y^2 + 18 x^2 y^2 - 12 y^3 + 9 y^4=0.$$

We can homogenize this to $$9 X^4 + 18 X^2 Y^2 + 9 Y^4 - 120 X^3 Z - 12 X^2 Y Z - 120 X Y^2 Z - 12 Y^3 Z + 402 X^2 Z^2 + 80 X Y Z^2 + 6 Y^2 Z^2 - 184 X Z^3 + 340 Y Z^3 - 711 Z^4=0$$

A degree $4$ curve usually has genus $3$. However, this curve has nodes at $(X:Y:Z) = (1:\pm i: 0)$, so it has genus $1$. If I cared enough about this problem, I would find a rational change of coordinates turning this into a cubic curve in Weierstrass form. The integral $y dx$ would then presumably be some sort of elliptic integral.

UPDATE:

Since distance is rotation symmetric, we can take any two points which are $\sqrt{5}$ apart and use them in place of $(2,3)$ and $(3,1)$. I choose the equation $$2 \sqrt{(x-r)^2 + y^2} + \sqrt{(x-r-\sqrt{5})^2 + y^2} = 4$$ where $r = (4-\sqrt{5})/3$. The reason for this choice is that it puts the point $(0,0)$ on the curve. Eliminating the radicals, this is $$128 x + (256 \sqrt{5} x)/3 + 16 x^2 - 96 \sqrt{5} x^2 - 48 x^3 + 24 \sqrt{5} x^3 + 9 x^4 - 128 y^2 - 32 \sqrt{5} y^2 - 48 x y^2 + 24 \sqrt{5} x y^2 + 18 x^2 y^2 + 9 y^4 \quad (\ast)$$ Again, we want to integrate $y dx$ around the connected component through $(0,0)$.

I make the substitution $u=x/(x^2+y^2)$, $v=y/(x^2+y^2)$. (Motivation: If we consider the homogenization of $(\ast)$, it is a degree $4$ equation passing through $(0:0:1)$ and with nodes at $(1:\pm i : 0 )$. Any conic through these points will intersect the curve at $3$ other points, since each node is a double intersection. A basis for the vector space of such conics is $XZ$, $YZ$, $X^2+Y^2$. So mapping the degree $4$ curve to $\mathbb{P}^2$ by $(X:Y:Z) \mapsto (XZ:YZ:X^2+Y^2)$ should have image a degree $3$ curve. Dehomogenizing, we want $(x,y) \mapsto (x/(x^2+y^2), y/(x^2+y^2))$.)

Conveniently, we can invert $u=x/(x^2+y^2)$, $v=y/(x^2+y^2)$ to give $x=u/(u^2+v^2)$, $y=v/(u^2+v^2)$. Plugging this into $(\ast)$ and clearing denominators, we now have the curve: $$(-384 - 96 \sqrt{5} + 384 u + 256 \sqrt{5} u) v^2 + (27 - 144 u + 72 \sqrt{5} u + 48 u^2 - 288 \sqrt{5} u^2 + 384 u^3 + 256 \sqrt{5} u^3)=0$$

Hooray! It's cubic, and almost in Weierstrass form. We now want to integrate, not $y dx$, but $$\frac{v}{u^2 + v^2} d \left( \frac{u}{u^2+v^2} \right) = \frac{v\left( (u^2+v^2 ) du - 2(u du + v dv)\right)}{(u^2+v^2)^3}.$$

Anyone want to turn this into a standard elliptic integral and put it out of its misery?

Simple minded approach.
Take an area of $3500 \times 3500$ pixels (not really) and count the pixels inside the curve.

DEPRECATED

VERY LATE UPDATE (2021).

Take an area of $3^{13} \times 3^{13}$ pixels (not really) and count the pixels inside the curve. The following method is supposed to work for any convex closed (Jordan) curve that divides the plane.
Start with a grid of $3 \times 3$ squares and refine the grid three times each iteration.
Leave out a square as soon as all four vertices of it are guaranteed to be inside or outside the curve.
Continue (recursively) with the remaining squares, until the maximum grid refinement is reached.
Pictures say more that a thousand words:

Free accompanying (Delphi Pascal) source code available at this place:

• MSE publications / references 2021

There is also a version available with binary instead of ternary refinement.
In addition, the software has been enhanced with Isoline Refinement. The meaning of the latter is best explained with a picture ($\color{red}{red}$ line):

With isoline refinement, a much better approximation of achille hui's result is found. Agreement is now within 12 digits. As a bonus, a value for the length of the perimeter is obtained as well:

Area = 3.20405086915288E+0000
Perimeter = 6.457419618 86746E+0000