Green's theorem

Let $${\mathbf E}({\mathbf x})=[E_x(x,y), E_y(x,y)]$$ be a smooth (differentiable) two-component vector field (or the pair of functions $$E_x(x,y), E_y(x,y)$$) on the two dimensional space then the line integral of the field projection onto the unite length vector anti-clockwise field $$\mathbf{m}(\mathbf{x})$$ always smoothly tangent to the closed curved over the arbitrary two dimensional closed curve $${\mathcal C}$$ equals the integral of the difference of the partial derivatives $$\frac{\partial E_y}{\partial x}$$, and $$\frac{\partial E_x}{\partial y}$$ over the plane region $$D$$ bounded inside the curve or otherwise the outside of the curve values of the field make virtually no contributions to the integral over the region providing that the field is sufficiently smooth that the second derivatives of the field components exists in the region i.e.

$$\oint\limits_{\mathcal C} \mathbf{E} \cdot d \mathbf{l}=\oint\limits_{\mathcal C} (E_x dx + E_y dy) = \iint_D \left(\frac{\partial E_y}{\partial x}-\frac{\partial E_x}{\partial y}\right)dx dy$$

where $$d \mathbf{l}= \mathbf{m} dl $$ and  $$D$$ is the region enclosed by the curve $$\mathcal C$$.

Proof
We can approximate the integral on the right side over the region by the finite sum by dividing densely the space around the region $$D$$ into small squares with the sides $$dx=dy$$ and the vertices $$[x_{i},y_{j}]$$ and approximating the bounding curve $$\mathcal C$$ of the region by the sides of squares which are the closet to the curve as well as the coordinate derivatives of the field $${\mathbf E}$$ by their difference quotients. We will keep the vertices coordinate names for the convenience even if they are equal and keep the square vertices coordinate indices $$i,j$$ even if they are limited by the region bounded by the curve.

We get

$$\iint_D \left(\frac{\partial E_y}{\partial x}-\frac{\partial E_x}{\partial y}\right)dx dy=\sum_{i,j} \left [ \frac{E_y(x_{i+1},y_{j}) - E_y(x_i,y_j)}{dx}-\frac{E_x({x_i,y_{j+1}}) - E_x(x_i,y_j)}{dy} \right ] dxdy+\Theta(dxdy),$$

Now the essential in proving the theorem is to focus on the contribution to the finite sum approximating the region integral from the one component of the $$\mathbf{E}$$ field itself and notice that because of the cancelation of the sign alternating term the sums reduce to only the end points. For example for $$E_y$$ and the fixed $$y$$-line $$j$$ and its length we have

$$\sum_{i} \left [ \frac{E_y(x_{i+1},y_{j}) - E_y(x_i,y_j)}{dx}\right ] dxdy =[E_y(x_n,y_j)-E_y(x_1,y_j)]dy $$,

Note that while $$dy$$ is an infinitesimal (small) linear element of the region boundary curve parallel to the $$y$$ axis and for the unite vector $$\mathbf{m}_y=[0,1]$$ parallel to it $$\mathbf{E} (x_{n},y_j) \cdot \mathbf{m}_y = E_y(x_{n},y_j)$$ and so for the second point with the minus sign the right side is an approximate to the growth $$\mathbf{E} \cdot d \mathbf{l}$$ of the counter-clockwise line integral $$\oint\limits_{\mathcal C} \mathbf{E} \cdot d \mathbf{l}$$ i.e.

$$[E_y(x_n,y_j)-E_y(x_1,y_j)]dy = \mathbf{E} \cdot d \mathbf{l} + \Theta(dy)$$.

Summing up all the all the $$dy$$ contributions over $$j$$ and repeating the considerations for the field component $$E_x$$ leading to the $$dx$$ contributions of the region boundary curve integral we get $$\iint_D \left(\frac{\partial E_y}{\partial x}-\frac{\partial E_x}{\partial y}\right)dx dy=\sum_{i,j} [E_y(x_{nj},y_{j})-E_y(x_{1j},y_{j})]dy-[E_x(x_{i},y_{ni})-E_x(x_{i},y_{1i})]dx+\Theta(dx)+\Theta(dy)$$

and so finally prove $$\oint\limits_{\mathcal C} \mathbf{E} \cdot d \mathbf{l}=\oint\limits_{\mathcal C} (E_x dx + E_y dy) = \iint_D \left(\frac{\partial E_y}{\partial x}-\frac{\partial E_x}{\partial y}\right)dx dy$$.