Quadrature Rules on Triangles

Quadrature rules on the interval $[-1, 1]$ is pretty well-documented in many different textbooks (see for example Cheney and Kincaid). By transformations, we can easily integrate across the whole real line spectrum. Going into two dimensions, we can take tensor products to integrate rectangles easily. On the other hand, integration on the 2D simplex (i.e. triangles) is a bit harder. There really exists two types of quadrature rules for the triangle:

Tensor products + transformations
“Symmetric quadrature rules”

The latter offers better convergence with smaller points, and has better spread-out points than the tensor product version, but offers less technical support through software packages; we refer the reader to the following paper. We will discuss the first, and exposit a little more; I will be borrowing heavily from Karniadakis and Sherwin’s fluid dynamics textbook.

Suppose our triangle at (-1, -1), (1, -1), and (-1, 1) (so called reference triangle) uses $x, y$ as it’s coordinates. If we use the change of coordinates $u = 2\frac{1 + x}{1 - y} - 1, v = y$ , it’ll map the triangle to the square defined by (-1, -1), (1, -1), (-1, 1) and (1,1). Hence, we can express the integral using a change of variables (details omitted)

$\int_T f(x, y) \, dx dy = \int_{-1}^1 \int_{-1}^1 f(u, v) \frac{1 - v}{2} \, du dv$ .

We note that the weird factor is simply the change of variables Jacobian. Now it’s easy to use Gaussian quadrature on this! But there are better version, as we have the Jacobian term; looking up the quadrature rules, we see that Gauss-Jacobi fits our bill. We refer the reader to page 144 of the Sherwin textbook to see the full-details.

A simple implementation is below (this is far from the best). I guess I’ll talk about this at the end of the week.

from scipy.special import roots_jacobi as rj
def quadtriangle(f, w00, w10, x00, x10):
    total = 0
    for i in range(len(w00)):
        subsum = 0
        for j in range(len(w10)):
            subsum += w10[j] / 2 * f([(1 + x00[i]) * (1 - x10[j]) / 2 - 1, x10[j]])
        total += w00[i] * subsum
    return total
def f(x):
    return x[0]**2 + x[1]**2
p = 4
x00, w00 = rj(p + 1, 0, 0)
x10, w10 = rj(p + 1, 1, 0)
print(quadtriangle(f, w00, w10, x00, x10))

Quadrature rules

Leave a Reply Cancel reply