Me ofrecen tres soluciones (como se pide en su propia respuesta) el Python paquete SciPy - (1) un simple doble integral técnica, (2) un polar integral, y (3) un camino integral usando el Verde del teorema. De estos, el segundo es sin duda la mejor por las razones que vamos a explorar, la primera es la más sencilla, y solo me ofrecen el tercer puesto fue solicitado específicamente. El valor correcto es de aproximadamente 12.711.
Una integral doble
Deje $\chi$ ser la función característica del conjunto. Es decir,
$$\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.
$$
A continuación, el área en cuestión es
$$\int_{-3}^{3} \int_{-3}^{3} \chi(x,y) \, dx \, dy.$$
Los límites de la integral puede ser de cualquiera de los números de modo que el rectángulo resultante contiene la región de interés, pero (desde la perspectiva de un integrador numérico, el más pequeño, mejor.
Podemos implementar esto en SciPy de la siguiente manera:
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)
Por lo tanto, el área debe ser (aproximadamente) $12.7075$. El segundo valor es un error de enlazado, que yo no confiar demasiado, ya que tales límites hacer suposiciones sobre la suavidad de el integrando. De hecho, genera un mensaje de advertencia. Por lo tanto, conceptualmente simple y fácil de programar para cualquier definidas implícitamente la región, esta técnica se obtiene un resultado que no es del todo fiable. No es demasiado sorprendente, ya que contamos con un multi-dimensional integrando que es discontinua en un poco complicada.
Otra forma de comprobar los resultados como este es examinar el resultado en otro sistema. Esto se puede hacer con un one-liner en Mathematica:
NIntegrate[Boole[x^2 + y^2 + Sin[4x] + Sin[4y] < 4], {x,-3,3}, {y,-3,3}]
(* Out: 12.6825 *)
El resultado difiere significativamente, pero, de nuevo, un mensaje de advertencia se emite y el resultado debe ser considerado sospechoso.
Una integral polar
Otro enfoque es el uso de la expresión polar para el área de:
$$\int_{\alpha}^{\beta} \frac{1}{2} f(\theta)^2 \, d\theta,$$
donde $\alpha\leq\theta\leq\beta$ (para un determinado $\theta$) $0\leq r\leq f(\theta)$. En nuestro caso, $\theta$ pistas de $0$ $2\pi$y, dado un $\theta$, se debe determinar la distancia del punto correspondiente en el límite de nuestra región.
![enter image description here]()
El radio de $f(\theta)$ puede ser calculada fácilmente para esta región el uso de un numérica de la ecuación de problemas; tenemos la suerte de que la región muestra suficiente regularidad como para permitir esto. Este esquema puede ser aplicado en SciPy de la siguiente manera:
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)
Tenga en cuenta que no hay ningún mensaje emitido y el error de estimación es mucho mejor. No es de extrañar, como hemos fácil de calcular, continuo, integrando dimensiones. Además, es fácil de implementar en 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 *)
La diferencia entre los dos sistemas está en el orden de $10^{-10}$. La única desventaja de este enfoque es que no es muy ampliamente aplicable ya que la región necesita para ser susceptibles de un polar descripción.
Verde del teorema de
Finalmente, usted puede usar el Verde del teorema como usted propone. Mientras que la integral es unidimensional, la única manera de parametrizar el límite es el uso de la interpolación. Entonces, usted está atascado, el cómputo numérico derivado de la interpolación en el integrando, lo que conduce a su propios problemas.
En cualquier caso, la idea detrás de este enfoque surge a partir del teorema de Green
$$
\iint_{\Omega} \left(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}\right) \, dA
=
\oint_{\partial\Omega} P \, dx + Q \, dy
$$
con $Q(x,y)=x$ $P(x,y)=0$ tenemos
$$\iint_{\Omega} 1 \, dA = \oint_{\partial\Omega} x \, dy.$$
Ahora, la integral de la izquierda representa claramente el área de $\Omega$ y la integral de la derecha puede ser fácilmente calculada como
$$\oint_{\partial\Omega} x \, dy = \int_a^b x(t)y'(t)\,dt,$$
si tenemos una parametrización de la $(x(t),y(t))$ del límite de $\Omega$. Podemos conseguir que la parametrización mediante la interpolación.
Para implementar esto en SciPy, primero vamos a generar una imagen de la región.
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])
Llamando plt.show() debe generar la siguiente imagen.
![enter image description here]()
Ahora, necesitamos los puntos de la imagen para la interpolación.
xy = contour_pic.collections[0].get_paths()[0].vertices
La interpolación es el siguiente
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')
Con estos puntos, podemos ahora calcular la 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)
El resultado no parece especialmente cerca de los resultados anteriores, tenemos un mal error de enlazado, y genera un mensaje de advertencia. Nada de esto es sorprendente, ya que el integrando es mucho más complicado que la versión polar.
