Finding the Area of an Implicit Relation – Integration

areaintegrationmultivariable-calculus

Let's say we have the function: $$x^2+y^2+\sin(4x)+\sin(4y)=4$$

I haven't taken Calculus III, in fact I'm just taking Calculus I. Since I learned how to find the derivative of implicit relations I was wondering if you can find the area of it.

Eventually I went into finding the area of closed shapes, and I found Green's Theorem.

How could you use Green's theorem to determine a contour path, use the line integral, and finally solve for the area of the closed implicit shape, such as $x^2+y^2+\sin(4x)+\sin(4y)=4$.

Is it also effective? Or could their be easier methods.

Best Answer

I offer three solutions using (as requested in your own answer) the Python's SciPy package - (1) a simple double integral technique, (2) a polar integral, and (3) a path integral using Green's theorem. Of these, the second is certainly the best for reasons that we'll explore, the first is the easiest, and I only offer the third since it was specifically requested. The correct value is approximately 12.711.

A double integral

Let $\chi$ be the characteristic function of the set. That is, $$\chi(x,y) = \left\{ \begin{array}{ll} 1 & \text{if } x^2 + y^2 + \sin(4x)+\sin(4y)<4 \\ 0 & \text{otherwise}. \end{array}\right. $$ Then the area in question is $$\int_{-3}^{3} \int_{-3}^{3} \chi(x,y) \, dx \, dy.$$ The bounds for the integral can be any numbers so that the resulting rectangle contains the region of interest but (from the perspective of a numerical integrator, the smaller the better.

We can implement this in SciPy as follows:

from scipy import sin
from scipy.integrate import dblquad
def f(x,y): 
    if x**2 + y**2 + sin(4*x)+sin(4*y)<4:
        return 1
    else:
        return 0
dblquad(f, -3,3, lambda x: -3, lambda x: 3)

#Out: (12.707523995238732, 0.0006145836908579838)

Thus, the area should be (approximately) $12.7075$. The second value is an error bound, which I wouldn't trust too much, since such bounds make assumptions on the smoothness of the integrand. In fact, a warning message is issued. Thus, while conceptually simple and easy to program for any implicitly defined region, this technique yields a result that is not particularly reliable. Not too surprising, as we have a multi-dimensional integrand that is discontinuous on a somewhat complicated set.

Another way to check results like this is to examine the result in another system. This can be done with a one-liner in Mathematica:

NIntegrate[Boole[x^2 + y^2 + Sin[4x] + Sin[4y] < 4], {x,-3,3}, {y,-3,3}]
(* Out: 12.6825 *)

The result differs significantly but, again, a warning message is issued and the result should be considered suspect.

A polar integral

Another approach is to use the polar expression for area: $$\int_{\alpha}^{\beta} \frac{1}{2} f(\theta)^2 \, d\theta,$$ where $\alpha\leq\theta\leq\beta$ and (for a given $\theta$) $0\leq r\leq f(\theta)$. In our case, $\theta$ runs from $0$ to $2\pi$ and, given a $\theta$, we must determine the distance of the corresponding point on the boundary of our region.

enter image description here

The radius $f(\theta)$ can be computed fairly easily for this region using a numerical equation solver; we are fortunate that the region displays sufficient regularity to allow this. This scheme can be implemented in SciPy as follows:

from scipy import sin, cos, pi
from scipy.optimize import brentq
from scipy.integrate import quad
def f(r,theta): return (r*cos(theta))**2 + (r*sin(theta))**2 + 
  sin(4*r*cos(theta))+sin(4*r*sin(theta)) - 4
def rr(theta): return brentq(lambda r: f(r,theta),1,3)
def integrand(theta): return 0.5*rr(theta)**2
quad(integrand, 0, 2*pi)

#Out: (12.710967450506423, 7.130808983545281e-08)

Note that there is no message issued and the error estimate is much better. Not surprising, as we have an easily computed, continuous, one-dimensional integrand. Also, it's easily implemented in Mathematica:

f[theta_?NumericQ] := r /. FindRoot[
  (r*Cos[theta])^2 + (r*Sin[theta])^2 + 
    Sin[4 r*Cos[theta]] + Sin[4 r*Sin[theta]] == 4, 
  {r, 1, 3}];
NIntegrate[f[theta]^2/2, {theta, 0, 2 Pi}] 
(* Out: 12.710967451423258 *)

The difference between the two systems is on the order of $10^{-10}$. The one disadvantage of this approach is that it is not very broadly applicable as the region needs to be amenable to a polar description.

Green's theorem

Finally, you can use Green's theorem as you propose. While the integral is one dimensional, the only way to parametrize the boundary is to use interpolation. Then, you're stuck computing a numerical derivative of the interpolation in the integrand, which leads to it's own problems.

At any rate, the idea behind this approach arises from Green's theorem

$$ \iint_{\Omega} \left(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}\right) \, dA = \oint_{\partial\Omega} P \, dx + Q \, dy $$ with $Q(x,y)=x$ and $P(x,y)=0$ we get $$\iint_{\Omega} 1 \, dA = \oint_{\partial\Omega} x \, dy.$$ Now, the integral on the left clearly represents the area of $\Omega$ and the integral on the right can be easily computed as $$\oint_{\partial\Omega} x \, dy = \int_a^b x(t)y'(t)\,dt,$$ if we have a parametrization $(x(t),y(t))$ of the boundary of $\Omega$. We can get that parametrization using interpolation.

To implement this in SciPy, first we'll generate an image of the region.

import matplotlib.pyplot as plt
import numpy as np

def f(x,y): return x**2 + y**2 + np.sin(4*x) + np.sin(4*y)
delta = 0.01
x = np.arange(-3,3,delta)
y = np.arange(-3,3,delta)
x,y=np.meshgrid(x,y)
z = f(x,y)
contour_pic = plt.contour(x,y,z,[4])

Calling plt.show() should generate the following image.

enter image description here

Now, we need the points in the image for interpolation.

xy = contour_pic.collections[0].get_paths()[0].vertices

The interpolation goes as follows

from scipy.interpolate import interp1d
t = np.linspace(0,1,xy.shape[0])
fx = interp1d(t,xy[:,0], fill_value=xy[0,0], bounds_error=False, kind='cubic')
fy = interp1d(t,xy[:,1], fill_value=xy[0,1], bounds_error=False, kind='cubic')

With these points, we can now compute the integral.

from scipy.misc import derivative
from scipy.integrate import quad

def integrand(t): return fy(t)*derivative(fx,t,0.002)
quad(integrand, 0,1)

# Out: (12.675772468194559, 0.05846541344794076)

The result does not look particularly close to the previous results, we have a poor error bound, and a warning message is issued. None of this is surprising, as the integrand is much more complicated than the polar version.

Related Question