1
$\begingroup$

When setting up a finite element system you have to use quadrature to calculate the integrals. I'm having trouble understanding what order rule to use.

I know of some rules of thumb, for example with the 1D Hemholtz you end up with the matrix $M + S$, where $M_{ij} = (\phi_i, \phi_j)$ and $S_{ij} = (\phi_i', \phi_j')$. In this case if $\phi$ is a polynomial of degree $n$, then the integrals in $M$ will of degree $2n$ and the integrals in $S$ will be of degree $2(n-1)$. Using a Gauss rule of order $n$ is recommended here, which integrates $S$ exactly.

When using DG on the advection equation the "stiffness" matrix becomes $S_{ij} = (\phi_i', \phi_j)$. In this case I found that using an order $n$ Gauss quadrature rule (which should still integrate $S$ exactly) gives me a singular matrix, and that I have to use a rule of order $n+1$ to get an invertible matrix. Is there some way to know beforehand what order rule to use, or it just determined experimentally?

EDIT: As per comments below, I'm using DG with Langrangian elements on a uniform grid to solve the time dependent advection equation $u_t + au_x = 0$. I have periodic boundary conditions on the domain $[0, 2\pi]$ and the initial condition $u(x,0) = \sin(x)$.

This is my first attempt coding DG myself, so I'm using central flux and backwards Euler for the time step. Using basis functions of degree $n$ and $\Delta t = h^{n+1}$, the solutions converge to the exact solution in the $L^2$ norm as $O(h^{n+1})$.

$\endgroup$
  • $\begingroup$ I'm guessing it's the quadrature because if I add the extra quadrature point the matrix is invertible and the solutions converge at the theoretical rate (at least for my simple test case). In 1D the interface terms are not integrals and are independent of the quadrature rule. $\endgroup$ – Lukas Bystricky Aug 27 '15 at 12:38
  • $\begingroup$ Small aside: for the advection equation $u_x = f$, the stiffness matrix should have entries $S_{ij} = (\phi_i,\phi_j')$ (test functions correspond to rows, not columns -- this is assuming standard mathematical notation that the first index corresponds to the row). $\endgroup$ – Christian Clason Aug 27 '15 at 13:04
  • $\begingroup$ On the choice of quadrature rules for DG applied to conservation laws from the point of view optimal accuracy, see this paper dx.doi.org/10.1090/S0025-5718-1990-1010597-0 $\endgroup$ – cpraveen Jun 21 '16 at 4:02
1
$\begingroup$

Since you have periodic boundary conditions, the system matrix (arising from the volume and flux terms) cannot have full rank (the solution of the stationary equation is unique only up to constants). Hence the well-posedness of the backward Euler step relies on the invertibility of the mass matrix (or at least that the two kernels have only a trivial intersection).

Now you can fairly easily show that the mass matrix will be positive definite if the quadrature is exact. As the product of two polynomials of order $n$ is $2n$, Gauß quadrature of order $n$ is not sufficient and you need at least order $n+1$ (as you have already noted).

The converse is not as obvious, but there's a bit of "finite element folklore" (at least in the engineering community) that each Gauß quadrature point adds $d$ (the space dimension) to the rank of the elemental mass matrix, and full rank is equal to the number of degrees of freedom in each element. In 1D, this is $n+1$, so you need $n+1$ points.

Usually, it is sufficient that the highest-order term (i.e., the one in the bilinear form involving the hightest-order derivatives) is integrated exactly (ignoring variable coefficients) and surface terms are integrated using one order higher, but if the highest-order term is itself not invertible (e.g., because of pure Neumann or periodic boundary conditions), you need to make sure the mass matrix is integrated exactly as well.

$\endgroup$
  • $\begingroup$ What do you mean by highest order term? If we integrate the highest order term exactly won't all the lower order terms also be integrated exactly? $\endgroup$ – Lukas Bystricky Aug 27 '15 at 16:27
  • $\begingroup$ The highest-order term is the one involving the highest occurring derivatives (such as $(\nabla \phi_i,\nabla \phi_j)$ for the second-order equation $-\Delta y +\nabla y = f$). If this term has the highest derivatives, it will have the lowest order polynomials to integrate, so this is the least restriction on the quadrature order. $\endgroup$ – Christian Clason Aug 27 '15 at 16:51
  • $\begingroup$ You're welcome; I hope it's useful. $\endgroup$ – Christian Clason Aug 28 '15 at 8:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.