The previous sections on approximating \( f \) by a finite element function \( u \) utilize the projection/Galerkin or least squares approaches to minimize the approximation error. We may, alternatively, use the collocation/interpolation method as described in the section Comparison with finite elements and interpolation/collocation. Here we shall compare these three approaches with what one does in the finite difference method when representing a given function on a mesh.
Approximating a given function \( f(x) \) on a mesh in a finite difference context will typically just sample \( f \) at the mesh points. If \( u_i \) is the value of the approximate \( u \) at the mesh point \( \xno{i} \), we have \( u_i = f(\xno{i}) \). The collocation/interpolation method using finite element basis functions gives exactly the same representation, as shown the section Comparison with finite elements and interpolation/collocation,
$$ u(\xno{i}) = c_i = f(\xno{i})\tp$$
How does a finite element Galerkin or least squares approximation differ from this straightforward interpolation of \( f \)? This is the question to be addressed next. We now limit the scope to P1 elements since this is the element type that gives formulas closest to those arising in the finite difference method.
The linear system arising from a Galerkin or least squares approximation reads in general
$$ \sum_{j\in\If} c_j (\baspsi_i,\baspsi_j) = (f,\baspsi_i),\quad i\in\If\tp $$ In the finite element approximation we choose \( \baspsi_i =\basphi_i \). With \( \basphi_i \) corresponding to P1 elements and a uniform mesh of element length \( h \) we have in the section Calculating the linear system calculated the matrix with entries \( (\basphi_i,\basphi_j) \). Equation number \( i \) reads
$$ \begin{equation} \frac{h}{6}(u_{i-1} + 4u_i + u_{i+1}) = (f,\basphi_i) \tp \tag{39} \end{equation} $$ The first and last equation, corresponding to \( i=0 \) and \( i=N \) are slightly different, see the section The structure of the coefficient matrix.
The finite difference counterpart to (39) is just \( u_i=f_i \) as explained in the section Finite difference approximation of given functions. To easier compare this result to the finite element approach to approximating functions, we can rewrite the left-hand side of (39) as
$$ \begin{equation} h(u_i + \frac{1}{6}(u_{i-1} - 2u_i + u_{i+1})) \tp \end{equation} $$ Thinking in terms of finite differences, we can write this expression using finite difference operator notation:
$$ [h(u + \frac{h^2}{6}D_x D_x u)]_i,$$ which is nothing but the standard discretization of
$$ h(u + \frac{h^2}{6}u'')\tp$$
Before interpreting the approximation procedure as solving a differential equation, we need to work out what the right-hand side is in the context of P1 elements. Since \( \basphi_i \) is the linear function that is 1 at \( \xno{i} \) and zero at all other nodes, only the interval \( [\xno{i-1},\xno{i+1}] \) contribute to the integral on the right-hand side. This integral is naturally split into two parts according to (25):
$$ (f,\basphi_i) = \int_{\xno{i-1}}^{\xno{i}} f(x)\frac{1}{h} (x - \xno{i-1}) dx + \int_{\xno{i}}^{\xno{i+1}} f(x)\frac{1}{h}(1 - (x - x_{i})) dx \tp $$ However, if \( f \) is not known we cannot do much else with this expression. It is clear that many values of \( f \) around \( \xno{i} \) contributes to the right-hand side, not just the single point value \( f(\xno{i}) \) as in the finite difference method.
To proceed with the right-hand side, we can turn to numerical integration schemes. The Trapezoidal method for \( (f,\basphi_i) \), based on sampling the integrand \( f\basphi_i \) at the node points \( \xno{i}=i h \) gives
$$ (f,\basphi_i) = \int_\Omega f\basphi_i dx\approx h\half( f(\xno{0})\basphi_i(\xno{0}) + f(\xno{N})\basphi_i(\xno{N})) + h\sum_{j=1}^{N-1} f(\xno{j})\basphi_i(\xno{j}) \tp $$ Since \( \basphi_i \) is zero at all these points, except at \( \xno{i} \), the Trapezoidal rule collapses to one term:
$$ \begin{equation} (f,\basphi_i) \approx hf(\xno{i}), \end{equation} $$ for \( i=1,\ldots,N-1 \), which is the same result as with collocation/interpolation, and of course the same result as in the finite difference method. For \( i=0 \) and \( i=N \) we get contribution from only one element so
$$ \begin{equation} (f,\basphi_i) \approx {\half}hf(\xno{i}),\quad i=0,\ i=N \tp \end{equation} $$
Simpson's rule with sample points also in the middle of the elements, at \( \xno{i+\half}=(\xno{i} + \xno{i+1})/2 \), can be written as
$$ \int_\Omega g(x)dx \approx \frac{\tilde h}{3}\left( g(\xno{0}) +
2\sum_{j=1}^{N-1} g(\xno{j})
+ 4\sum_{j=0}^{N-1} g(\xno{j+\half}) + f(\xno{2N})\right),
$$
where \( \tilde h= h/2 \) is the spacing between the sample points.
Our integrand is \( g=f\basphi_i \). For all the node points,
\( \basphi_i(\xno{j})=\delta_{ij} \), and therefore
\( \sum_{j=1}^{N-1} f(\xno{j})\basphi_i(\xno{j})=f(\xno{i}) \).
At the midpoints, \( \basphi_i(\xno{i\pm\half})=1/2 \) and
\( \basphi_i(\xno{j+\half})=0 \) for \( j>1 \) and \( j
$$
\begin{equation}
(f,\basphi_i) \approx
\frac{h}{3}(f_{i-\half} + f_i + f_{i+\half})
\tp
\end{equation}
$$
This result shows that, with Simpson's rule, the finite element method
operates with the average of \( f \) over three points, while the finite difference
method just applies \( f \) at one point. We may interpret this as
a "smearing" or smoothing of \( f \) by the finite element method.
We can now summarize our findings. With the approximation of
\( (f,\basphi_i) \) by the Trapezoidal rule, P1 elements give rise
to equations that can be expressed as a finite difference
discretization of
$$
\begin{equation}
u + \frac{h^2}{6} u'' = f,\quad u'(0)=u'(L)=0,
\end{equation}
$$
expressed with operator notation as
$$
\begin{equation}
[u + \frac{h^2}{6} D_x D_x u = f]_i\tp \end{equation}
$$
As \( h\rightarrow 0 \), the extra term proportional to \( u'' \) goes to zero,
and the two methods are then equal.
With the Simpson's rule, we may say that we solve
$$
\begin{equation}
[u + \frac{h^2}{6} D_x D_x u = \bar f]_i,
\end{equation}
$$
where \( \bar f_i \) means the average \( \frac{1}{3}(f_{i-1/2} + f_i + f_{i+1/2}) \).
The extra term \( \frac{h^2}{6} u'' \) represents a smoothing effect: with
just this term, we would find \( u \) by integrating \( f \) twice and thereby
smooth \( f \) considerably. In addition, the finite element
representation of \( f \) involves an average, or a smoothing, of \( f \) on
the right-hand side of the equation system. If \( f \) is a noisy
function, direct interpolation \( u_i=f_i \) may result in a noisy \( u \)
too, but with a Galerkin or least squares formulation and P1 elements,
we should expect that \( u \) is smoother than \( f \) unless \( h \) is very
small.
The interpretation that finite elements tend to smooth the solution
is valid in applications far beyond approximation of 1D functions.
With a simple trick, using numerical integration, we can easily produce
the result \( u_i=f_i \) with the Galerkin or least square formulation
with P1 elements. This is useful in many occasions when we deal
with more difficult differential equations and want the finite element
method to have properties like the finite difference method (solving
standard linear wave equations is one primary example).
We have already seen that applying the Trapezoidal rule to the
right-hand side \( (f,\basphi_i) \) simply gives \( f \) sampled at \( \xno{i} \).
Using the Trapezoidal rule on the matrix entries
\( A_{i,j}=(\basphi_i,\basphi_j) \) involves a sum
$$ \sum_k \basphi_i(\xno{k})\basphi_j(\xno{k}),$$
but \( \basphi_i(\xno{k})=\delta_{ik} \) and
\( \basphi_j(\xno{k})=\delta_{jk} \).
The product \( \basphi_i\basphi_j \) is then different from zero only
when sampled at \( \xno{i} \) and \( i=j \). The Trapezoidal
approximation to the integral
is then
$$ (\basphi_i,\basphi_j) \approx h,\quad i=j,$$
and zero if \( i\neq j \). This means that we have obtained a diagonal matrix!
The first and last diagonal elements, \( (\basphi_0,\basphi_0) \) and
\( (\basphi_N,\basphi_N) \) get contribution only from the first and last
element, respectively, resulting in the approximate integral value \( h/2 \).
The corresponding right-hand side also has a factor \( 1/2 \) for \( i=0 \) and \( i=N \).
Therefore, the least squares or Galerkin approach with P1 elements and
Trapezoidal integration results in
$$ c_i = f_i,\quad i\in\If\tp $$
Simpsons's rule can be used to achieve a similar result for P2 elements, i.e,
a diagonal coefficient matrix, but with the previously derived
average of \( f \) on the right-hand side.
Identical results to those above will arise if we perform elementwise
computations. The idea is to use the Trapezoidal rule on the reference
element for computing the element matrix and vector. When assembled,
the same equations \( c_i=f(\xno{i}) \) arise. Exercise 19: Use the Trapezoidal rule and P1 elements encourages you to carry out the
details.
The matrix with entries \( (\basphi_i,\basphi_j) \) typically arises from
terms proportional to \( u \) in a differential equation where \( u \) is the
unknown function. This matrix is often called the mass matrix,
because in the early days of the finite element method, the matrix
arose from the mass times acceleration term in Newton's second law of
motion. Making the mass matrix diagonal by, e.g., numerical
integration, as demonstrated above, is a widely used technique and is
called mass lumping. In time-dependent problems it can sometimes
enhance the numerical accuracy and computational efficiency of the
finite element method. However, there are also examples where mass
lumping destroys accuracy.
Making finite elements behave as finite differences
Computations in physical space
Elementwise computations
Terminology