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.
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.
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.
Best Answer
@Valerie:
Does this figure help?
Just take the average of the $y$ (vertical) values at the left value of $x$ and right value of $x$--the green points. If this average is above the squares midpoint, then the larger area is below.