# computing higher order derivatives with linear elements

Consider the following equation on $$(0,1)$$, with Dirichlet boundary conditions on both ends.

$$\frac{d}{dx}\left(k(x)\frac{du}{dx}\right) = 0$$

Let us solve this using simple linear finite elements. The computed displacement $$u$$ is linear between nodes and is globally piecewise linear. Between nodes, we can take one derivative to obtain $$u_x$$. Higher order derivatives $$u_{xx},u_{xxx}$$ will be identically zero. Since linear finite elements are being used, I wouldn't expect to be able to calculate beyond $$u_{x}$$. So far so good.

Now, if I calculate $$u_x$$ by a suitable minimization procedure (detailed below), will it be possible to calculate beyond $$u_x$$ (i.e. $$u_{xx},u_{xxx}$$ etc ). I think not. Ultimately one is tied to the finite element basis that one chooses. I'm asking this question because I'm looking for someone to confirm this.

Calculating $$u_x$$ by a minimization procedure:

Find a $$\phi(x)$$ such that minimizes

$$\pi(x) = \frac{1}{2}\int_0^1 (\phi(x) - u_x)^2 dx$$

Taking variations in $$\phi$$

$$\int_0^1 \delta\phi\,\phi = \int_0^1 \delta\phi\,u_x \tag{1}$$

We can put in a linear finite element expansion for $$\phi,\delta\phi$$ (i.e. $$\phi = \sum N_A(x)\phi_A$$) and thus can calculate nodal values of $$u_x$$. Now if we extend the same procedure further and try to find $$u_{xx}$$.

Find a $$\psi$$ that minimizes

$$\psi(x) = \frac{1}{2}\int_0^1 (\psi(x) - \phi_x)^2 dx$$

Taking variations wrt $$\psi$$ yields,

$$\int_0^1 \psi\delta\psi = \int_0^1 \phi_x \delta\psi \tag{2}$$

However, equation (1) implies $$\phi\approx u_x$$ and hence $$\phi_x \approx u_{xx} \approx 0$$ (since a linear finite element expansion was used for $$u$$). Which implies $$\psi\approx0$$. So, even if we use the minimization procedure, we cannot get higher order derivatives, one is tied to the finite element basis one chooses at the start.

Is this line of thinking correct?

What you are trying to do has a name: it's called "recovery" and is a particular case of "postprocessing" the solution.

In your approach of finding an approximation to derivatives by minimization, remember that intuitively you think that $$\phi\approx u_x$$, but in fact you would choose $$\phi$$ to be from a piecewise linear (and thus continuous) space whereas $$u_x$$ is a piecewise constant and discontinuous. That is, the two objects are from different spaces and so you shouldn't expect them to be approximately equal in a pointwise sense. $$\phi$$ is simply the closest continuous piecewise linear function to $$u_x$$.

In practice, however, we typically do this via projection, not minimization. That is, we compute $$\phi = \Pi_1 u_x$$ where $$\Pi_1$$ is a projection operation that defines $$\phi$$ via $$\int \psi \phi = \int \psi u_x \qquad \forall \psi \in W_h.$$ Here, $$W_h$$ is the finite element space of $$\phi$$ into which you want to project the derivative. This also works for second derivatives: $$p = \Pi_2 u_{xx}$$ if you define $$\Pi_2$$ via integration by parts: $$\int qp = -\int q_x u_x + \text{boundary terms}\qquad \forall q \in P_h.$$ One has to be a bit careful with the boundary terms because integration-by-parts can only be done on a cell-by-cell basis, and so the "boundary terms" in this formula are likely going to include the jump of the gradient of $$u_x$$ between cells times something that involves the test function $$q$$ at these cell interfaces.

Projecting second derivatives this way is as far as you can go because you can reduce second derivatives to first derivatives by integration by parts. You can no longer do this with third derivatives, but they are not often needed.

• Very neat application of integration by parts. I was told in grad school: when in doubt use integration by parts. I guess I didn't learn my lesson.
– NNN
Jan 3, 2023 at 4:43
• @NNN In the various fields of analysis, my half-serious statement to students is "if you don't know what to do, there are two options of which one will work 90% of the time: (i) you integrate by parts, (ii) you do a Taylor expansion." Jan 4, 2023 at 14:58
• Thanks for you explanation. I am wondering how software tools like ANSYS can make it happen that a recovery of second order derivatives (like stresses) are accurate for linear elements? Do you have an idea how that's done, without using several degrees of freedom for a node? Jan 31, 2023 at 7:35
• @user21 This procedure is called "gradient recovery" and is based on averaging the gradients of neighboring cells. The procedure was introduced in a paper by Zienkiewicz and Zhu in 1992, see here: onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620330702 Jan 31, 2023 at 17:48
• @WolfgangBangerth Interesting that this technique was introduced relatively recently in 1992. I thought it would have been introduced in the 60s or 70s when FEM started to develop.
– NNN
Feb 2, 2023 at 4:05

Your statement is true for the onedimentional case with linear elements.

However, if you consider higher dimensions with non-constant metrics (bilinear or trilinear) then your seconde derivatives won't necessarily vanish.

Moreover, if you drop (reinterpret) the element local perpective and assume smooth data across element boundaries, no one prevents you of using neighboring elements to calculate higher derivatives. This is quite common in finite volume methods.