When going from the strong form of a PDE to the FEM form it seems one should always do this by first stating the variational form. To do this you multiply the strong form by an element in some (Sobolev) space and integrate over your region. This I can accept. What I don't understand is why one also has to use Green's formula (one or several times).

I've mostly been working with Poisson's equation, so if we take that (with homogenous Dirichlet boundary conditions) as an example, i e

$$ \begin{align} -\nabla^2u &= f,\quad u\in\Omega \\ u &= 0, \quad u\in\partial\Omega \end{align} $$

then it is claimed that the correct way to form the variational form is

$$ \begin{align} \int_\Omega fv\,\mathrm{d}\vec{x} &= -\int_\Omega\nabla^2 uv\,\mathrm{d}\vec{x} \\ &=\int_\Omega\nabla u\cdot\nabla v\,\mathrm{d}\vec{x} - \int_{\partial\Omega}\vec{n}\cdot\nabla u v\,\mathrm{d}\vec{s} \\ &=\int_\Omega\nabla u\cdot\nabla v\,\mathrm{d}\vec{x}. \end{align} $$

But what stops me from using the expression on the first line, isn't that also a variational form that can be used to get a FEM form? Isn't it corresponding to the bilinear and linear forms $b(u,v)=(\nabla^2 u, v)$ and $l(v)=(f, v)$? Is the problem here that if I use linear basis functions (shape functions) then I'll be in trouble because my stiffness matrix will be the null matrix (not invertible)? But what if I use non-linear shape functions? Do I still have to use Green's formula? If I don't have to: is it advisable? If I don't, do I then have a variational-but-not-weak formulation?

Now, let's say that I have a PDE with higher order derivatives, does that mean that there are many possible variational forms, depending on how I use Green's formula? And they all lead to (different) FEM approximations?


3 Answers 3


Short answer:

No, you don't have to do integration for certain FEMs. But in your case, you have to do that.

Long answer:

  • Let's say $u_h$ is the finite element solution. If you choose piecewise linear polynomial as your basis, then taking $\Delta$ on it will give you an order 1 distribution (think taking derivative on a Heaviside step function), and the integration of $-\Delta u_h\in H^{-1}$ multiplying with $v$ will only make sense when you take it as a duality pair rather than an $L^2$-inner product. You will neither get a null matrix, the Riesz representation theorem says that there is an element in $\varphi_{-\Delta u_h} \in H^1_0$ can characterize the duality pair by the inner product in $H^1$: $$ \langle-\Delta u_h ,v \rangle_{H^{-1},H^1_0} = \underbrace{\int_{\Omega}\nabla \varphi_{-\Delta u_h} \cdot \nabla v}_{\text{inner product in }H^1}. $$ Integrating by parts element by element for $u_h$ will shed a light on this duality pair: for $T$ an element in this triangulation $$ \int_{\Omega}\nabla u_h \cdot \nabla v = -\sum_{T}\left(\int_{T} \Delta u_h\,v + \int_{\partial T}\frac{\partial u_h}{\partial n}v\,dS\right), $$ this tells you that $-\Delta u_h$ should include inter-element flux jump in its duality pair representation, notice the integration on the boundary of each element is also a duality pair between $H^{1/2}$ and $H^{-1/2}$. Even if you use quadratic basis, which has a non-vanishing $\Delta$ on each element, you still can't write $(\Delta u, v)$ as an inner product, because of this inter-element flux jump's presence.

  • Integration by parts can be traced back to the Sobolev theory for elliptic pde using smooth function, where the $W^{k,p}$-spaces are all closure of smooth functions under the $W^{k,p}$ type of integral norm. Then people say what is the minimum regularity here that we can perform inner product. Also bearing in mind that an $H^1$-regular weak solution under certain condition is the $H^2$-strong solution (elliptic regularity). But piecewise continuous linear polynomial is not $H^2$, from this point of view, it doesn't make any sense to take inner product using $\Delta u_h$ either.

  • For certain FEMs, you don't have to do integration by parts. For example, Least-square finite element. Write the second order pde as a first order system: $$ \begin{cases} \boldsymbol{\sigma} = -\nabla u, \\ \nabla \cdot \boldsymbol{\sigma} = f. \end{cases} $$ Then you wanna minimize the least-square functional: $$ \mathcal{J}(v) = \|\boldsymbol{\sigma} + \nabla u\|_{L^2{\Omega}}^2 + \|\nabla \cdot \boldsymbol{\sigma} - f\|_{L^2{\Omega}}^2, $$ bearing the same spirit with Ritz-Galerkin functional, the finite element formulation of minimizing above functional in a finite element space does not require integration by parts.

  • $\begingroup$ when integration by parts is a must, what if we don't use it, but instead keep the highest order terms and do the FEA. Is that allowed and what impact? $\endgroup$
    – feynman
    Commented Aug 11, 2022 at 5:06
  • $\begingroup$ @feynman Then in conforming methods the finite element spaces should "conform" to higher derivative, which introduces much more degrees of freedom to ensure inter-element continuity, so that the scheme is consistent. If you wanna sacrifice consistency and use nonconforming methods, there will be various inter-element jump terms. Check the literature on FEM for non-divergence form elliptic equations. $\endgroup$
    – Shuhao Cao
    Commented Aug 12, 2022 at 4:15

Nothing stops you from doing that technically, but when you integrate by parts you get more flexibility with the solution space in that they need not have $H^2$ regularity (required for the non I.B.P formulation). The linear elements you suggest generally have enforced continuity between elements, and so could not be in $H^2$. The I.B.P formulation furthermore is symmetric, which has some of its own advantages as well.

  • 2
    $\begingroup$ Are you saying that linear shape functions gives a solution to the FEM formulation that is not in $H^2$ because differentiating this FEM solution twice (weakly) gives a sum of delta distributions, which is not in $L^2$? Does that mean that for pde:s of higher order than 2 I must use shape functions of higher order than 1 (at least if the test and trial spaces should be the same?)? $\endgroup$
    – Christian
    Commented Jun 30, 2013 at 1:24
  • 1
    $\begingroup$ What you're saying is essentially right. As for higher than second order PDE you don't necessarily have to use higher regularity spaces as writing down the mixed formulation (see Shuhao's answer) can help. You could also use other techniques like jump penalization to avoid this difficulty. For a classical FEM answer though, yes you would need higher regularity. $\endgroup$ Commented Jun 30, 2013 at 15:19
  • 2
    $\begingroup$ Let me stress the importance of symmetry. If a differential operator is self-adjoint, I expect to end end up with a symmetric matrix. Without integration by parts this will not be the case. $\endgroup$
    – Stefano M
    Commented Jul 2, 2013 at 8:02
  • 1
    $\begingroup$ Computational benefits was my primary thought in adding that, but are there also strong theoretical benefits of symmetry (aside from easier proofs of facts that likely still hold in the elliptic case, even if the discretization is nonsymmetric)? $\endgroup$ Commented Jul 2, 2013 at 15:08

Excellent answers already on this page, but there is still a (small) missing point.

The OP asked:

Now, let's say that I have a PDE with higher order derivatives, does that mean that there are many possible variational forms, depending on how I use Green's formula? And they all lead to (different) FEM approximations?

Integrating by parts (in the correct way) is important when you have Neumann type of boundary conditions. In fact it is by the i.b.p that you take into account the Neumann b.c in your variational formulation. The form of the Neumann b.c depends on how you integrate by parts, cf. this answer on integration by parts in linear elasticity. So even for second order elliptic PDE's, integration by parts has to be performed in a given way, in order to recover a variational formulation valid for Neumann or mixed boundary conditions. (And this of course irrespective of the fact that you discretize by FEM).

In mathematical physics, where the Neumann b.c. have a well defined meaning (heat flux, stress...), integration by parts is important in order to maintain the correct interpretation of the results. Even for homogeneous Dirichlet conditions and FEM this is true, since if we use a Lagrange multiplier method to impose the b.c.'s, the multipliers become physical quantities, like concentrated fluxes or forces.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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