In a previous post, I have described how one can integrate polynomials over polygonal curves in two dimensions. This post extends the previous one by describing how one can integrate polynomials over polygonal domains (whose boundaries do not self-intersect).

Let $p_d(x,y)$ be a polynomial of degree $d$: $$ p_d(x,y) = \sum_{m,n = 0}^{d} c_{mn}x^m y ^n $$ For instance, if $d=2$: $$ p_2(x,y) = c_{00} + c_{10}x + c_{01}y + c_{20}x^2 + c_{11}xy + c_{02}y^2 \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_polyn_deg_2} $$ Let $\Omega$ be a two-dimensional polygonal domain and let $S$ be the positively oriented boundary of $\Omega$. We wish to compute the area integral of $p_d(x,y)$ over $\Omega$: $$ \int_{\Omega} p_d(x,y) dxdy \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_p_omega} $$ The curve $S$ is a chain of $n$ segments $S_i$, $i = 1,2,\ldots,n$. Let $L_i$ be the length of $S_i$ and let $(x_i,y_i)$ and $(x_i + \Delta{x_i},y_i + \Delta{y_i})$ be the coordinates of its end nodes (see figure 1).

To compute $\int_{\Omega} p_d(x,y) dxdy$, we can use the divergence theorem. First notice that: $$ p_d(x,y) = \frac{d}{dx}\left(\sum_{m,n = 0}^{d} c_{mn} \frac{x^{m+1} }{m+1}y^n\right) = \nabla\cdot\binom{q_{d+1}(x,y)}{0} \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_pd_div} $$ where $$ q_{d+1}(x,y) = \sum_{m,n = 0}^{d} c_{mn} \frac{x^{m+1} }{m+1}y^n \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_def_qdp1} $$ is a polynomial of degree $(d+1)$. Integrating both sides of equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_pd_div} yields: $$ \begin{eqnarray} \int_{\Omega} p_d(x,y) dxdy &=& \int_{\Omega} \nabla\cdot\binom{q_{d+1}(x,y)}{0} dxdy \nonumber\\[5pt] &=& \int_{S} \binom{q_{d+1}(x,y)}{0} \cdot {\bf n}(x,y) dl \nonumber\\[5pt] &=& \int_{S} q_{d+1}(x,y)n^1(x,y) dl \nonumber\\[5pt] &=& \sum_{i=1}^{n} \int_{S_i} q_{d+1}(x,y)n^1(x,y) dl \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_qdp1n1} \end{eqnarray} $$ where ${\bf n}(x,y)$ is the outward unit normal to the surface $S$ (remember $S$ is positively oriented) and $n^1(x,y)$ is its $x$ component. Over the extension of each segment, the normal ${\bf n}(x,y)$ is constant. Denoting the normal along the segment $S_i$ by ${\bf n}_i$, we can then write equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_qdp1n1} as: $$ \int_{\Omega} p_d(x,y) dxdy = \sum_{i=1}^{n} n^1_i \int_{S_i} q_{d+1}(x,y) dl \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_qdp1} $$ An expression for ${\bf n}_i$ can be obtained by rotating the vector $(\Delta{x_i},\Delta{y_i})$ by $\pi/2$ in the clockwise direction (remember that $(\Delta{x_i},\Delta{y_i}) = (x_{i+1},y_{i+1}) - (x_i,y_i)$ is the vector which connects the two end points of $S_i$) and normalizing the resulting vector: $$ \displaystyle{\bf n}_i = \frac{(\Delta{y_i},-\Delta{x_i})}{\sqrt{\Delta{y_i}^2 + \Delta{x_i}^2}} \Longrightarrow n^1_i = \frac{\Delta{y_i}}{\sqrt{\Delta{x_i}^2 + \Delta{y_i}^2}} = \frac{\Delta{y_i}}{L_i} \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_n_i} $$ Although the right-hand side of equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_qdp1} can be computed using the method described in the previous post mentioned above, we will use another approach which yields a simpler expression for the area integral we wish to compute.

To begin, notice that segments having $\Delta{y_i} = 0$ (horizontal segments) do not contribute to the sum in equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_qdp1} since for these segments $n^1_i = 0$. In other words, we need to consider only the segments $S_i$ for which $\Delta{y_i} \neq 0$.

We can parameterize each segment $S_i$ as shown below: $$ {\bf r}_i(s) = \left( \tilde{x}_i(s),\tilde{y}_i(s) \right) = (x_i + s\Delta{x_i}, y_i + s\Delta{y_i}),\quad s \in [0,1] \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_param_segm} $$ Let ${\bf r}_i'(s)$ be the derivative of ${\bf r}_i(s)$ with respect to $s$. The length of ${\bf r}_i'(s)$ is: $$ \big\|{\bf r}'_i(s)\big\| = \big\|(\Delta{x_i},\Delta{y_i})\big\| = \sqrt{\Delta{x_i}^2 + \Delta{y_i}^2} = L_i $$ We can now use the parameterization from equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_param_segm} as well as the definition of line integral to compute the right-hand side of equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_qdp1}: $$ \begin{eqnarray} \int_{\Omega} p_d(x,y) dxdy &=& \sum_{i=1}^{n} n^1_i \int_{S_i} q_{d+1}(x,y) dl \nonumber\\[5pt] &=& \sum_{i=1}^{n} n^1_i \int_{0}^{1} q_{d+1}(\tilde{x}_i(s),\tilde{y}_i(s)) \big\|{\bf r}'_i(s)\big\|ds \nonumber\\[5pt] &=& \sum_{i=1}^{n} n^1_i \int_{0}^{1} q_{d+1}(x_i + s\Delta{x_i}, y_i + s\Delta{y_i}) L_i ds \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_line_int_q} \end{eqnarray} $$ Consider now the following change of variable (from now on we will explicitly consider only segments with $\Delta{y_i} \neq 0$): $$ y = y_i + s\Delta{y_i} \Longrightarrow dy = \Delta{y_i}ds \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_y_func_s} $$ As for the $x$ part, \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_y_func_s} implies: $$ \displaystyle x_i + s\Delta{x_i} = x_i + (y - y_i)\frac{\Delta{x_i}}{\Delta{y_i}} = x_i + m_i(y - y_i), \quad m_i = \frac{\Delta{x_i}}{\Delta{y_i}} \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_x_func_s} $$ Using equations \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_y_func_s} and \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_x_func_s} on equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_line_int_q}, we obtain: $$ \begin{eqnarray} \int_{\Omega} p_d(x,y) dxdy &=& \sum_{\underset{\Delta{y_i} \neq 0}{i=1}}^{n} n^1_i L_i \int_{y_i}^{y_{i+1}} q_{d+1}(x_i + m_i(y - y_i), y) \frac{1}{\Delta{y_i}} dy \nonumber\\[5pt] &=& \sum_{\underset{\Delta{y_i} \neq 0}{i=1}}^{n} \frac{n^1_i L_i}{\Delta{y_i}} \int_{y_i}^{y_{i+1}} q_{d+1}(x_i + m_i(y - y_i), y) dy \nonumber\\[5pt] &=& \sum_{\underset{\Delta{y_i} \neq 0}{i=1}}^{n} \int_{y_i}^{y_{i+1}} q_{d+1}(x_i - m_iy_i + m_i y, y) dy \nonumber\\[5pt] &=& \sum_{\underset{\Delta{y_i} \neq 0}{i=1}}^{n} \sum_{m,n=0}^{d} c_{mn} \int_{y_i}^{y_{i+1}} \frac{ (x_i - m_iy_i + m_i y)^{m+1} }{m+1} y^n dy \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_pd_interm} \end{eqnarray} $$ where equations \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_n_i} and \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_def_qdp1} were used. Since: $$ (x_i - m_iy_i + m_iy)^{m+1} = \sum_{p=0}^{m+1} \binom{m+1}{p} (x_i - m_iy_i)^p m_i^{m+1-p} y^{m+1-p} $$ then equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_pd_interm} can be written as: $$ \int_{\Omega} p_d(x,y) dxdy = \sum_{\underset{\Delta{y_i} \neq 0}{i=1}}^{n} \sum_{m,n=0}^{d} \sum_{p=0}^{m+1} \beta_{mnip}\int_{y_i}^{y_{i+1}} y^{m+1-p+n} dy \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_pd_interm2} $$ where (the reason for defining $\beta_{mnip} = 0$ when $\Delta{y_i} = 0$ will become clear soon): $$ \beta_{mnip} = \begin{cases} c_{mn} \binom{m+1}{p} \frac{(x_i - m_iy_i)^p}{m+1} m_i^{m+1-p} & \textrm{ if } \Delta{y_i} \neq 0\\[5pt] 0 & \textrm{ otherwise } \end{cases} \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_beta} $$ where $m_i$, defined in equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_x_func_s}, is given by: $$ \displaystyle m_i = \frac{\Delta{x_i}}{\Delta{y_i}} $$ Finally, computing the integrals over $y$ on equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_int_pd_interm2} gives us our final result: $$ \boxed{ \displaystyle \int_{\Omega} p_d(x,y) dxdy = \sum_{i=1}^{n} \sum_{m,n=0}^{d} \sum_{p=0}^{m+1} \beta_{mnip}\frac{\left(y_{i+1}^{m+2-p+n} - y_i^{m+2-p+n}\right)}{m+2-p+n} } \label{post_f8cf6364586db30c14e1742f378bbc1c_eq_symb_integ_poly_over_domain} $$ On the equation above, the terms corresponding to segments with $\Delta{y_i} = 0$ vanish automatically since for those terms we have $\beta_{mnip} = 0$. This allows us to remove the restriction $\Delta{y_i} \neq 0$ from the sum over segments.

Unfortunately equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_symb_integ_poly_over_domain} is not very efficient for numerical computations (although it can be used as a last resource). Whenever possible, it is better to compute the integrals explicitly and hard-code them. For example, if for a given problem all polyomials that must be integrated have degree $d=2$ (as in equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_polyn_deg_2}), then one can use equation \eqref{post_f8cf6364586db30c14e1742f378bbc1c_eq_symb_integ_poly_over_domain} to obtain a formula which is easy to hard-code: $$ \begin{eqnarray} \int_{\Omega} p_2(x,y) dxdy = \sum_{i=1}^{n} \Delta{y_i} \bigg[ \; & & c_{00} \left( \frac{\Delta{x_i}}{2} + x_i \right) \nonumber\\[5pt] \; &+ & c_{10} \left( \frac{\Delta{x_i}^2}{6} + \frac{\Delta{x_i} x_i}{2} + \frac{x_i^2}{2} \right) \nonumber\\[5pt] \; &+ & c_{01} \left( \frac{\Delta{x_i}\Delta{y_i}}{3} + \frac{\Delta{x_i} y_i + \Delta{y_i} x_i}{2} + x_iy_i \right) \nonumber\\[5pt] \; &+ & c_{20}\left( \frac{\Delta{x_i}^3}{12} + \frac{x_i\Delta{x_i}^2}{3} + \frac{x_i^2\Delta{x_i}}{2} + \frac{x_i^3}{3} \right) \nonumber\\[5pt] \; &+ & c_{11}\bigg( \frac{\Delta{x_i}^2\Delta{y_i}}{8} + \frac{\Delta{x_i}\Delta{y_i} x_i}{3} + \frac{x_i^2 \Delta{y_i}}{4} \nonumber\\[5pt] \; & & \quad \quad + \frac{\Delta{x_i}^2 y_i}{6} + \frac{\Delta{x_i} x_i y_i}{2} + \frac{x_i^2 y_i}{2} \bigg) \nonumber\\[5pt] \; &+ & c_{02} \bigg( \frac{\Delta{x_i}\Delta{y_i}^2}{4} + \frac{2\Delta{x_i}\Delta{y_i} y_i}{3} + \frac{\Delta{x_i} y_i^2}{2} \nonumber\\[5pt] \; & & \quad \quad + \frac{x_i\Delta{y_i}^2}{3} + \Delta{y_i} y_i x_i + x_i y_i^2 \bigg) \bigg] \end{eqnarray} $$