Study Guide: Introduction to Finite Element Methods

Why finite elements?

• Can with ease solve PDEs in domains with complex geometry
• Can with ease provide higher-order approximations
• Has (in simpler stationary problems) a rigorus mathematical analysis framework (not much considered here)

Domain for flow around a dolphin

The flow

Basic ingredients of the finite element method

• Transform the PDE problem to a variational form
• Define function approximation over finite elements
• Use a machinery to derive linear systems
• Solve linear systems

Our learning strategy

• Start with approximation of functions, not PDEs
• Introduce finite element approximations
• See later how this is applied to PDEs

Reason: the finite element method has many concepts and a jungle of details. This strategy minimizes the mixing of ideas, concepts, and technical details.

Approximation in vector spaces

Approximation set-up

General idea of finding an approximation \( u(x) \) to some given \( f(x) \):

 
$$ \begin{equation} u(x) = \sum_{i=0}^N c_i\baspsi_i(x) \tag{1} \end{equation} $$

where

• \( \baspsi_i(x) \) are prescribed functions
• \( c_i \), \( i=0,\ldots,N \) are unknown coefficients to be determined

How to determine the coefficients?

We shall address three approaches:

• The least squares method
• The projection (or Galerkin) method
• The interpolation (or collocation) method

Approximation of planar vectors; problem

Given a vector \( \f = (3,5) \), find an approximation to \( \f \) directed along a given line.

Approximation of planar vectors; vector space terminology

 
$$ \begin{equation} V = \mbox{span}\,\{ \psib_0\} \end{equation} $$

• \( \psib_0 \) is a basis vector in the space \( V \)
• Seek \( \u = c_0\psib_0\in V \)
• Determine \( c_0 \) such that \( \u \) is the "best" approximation to \( \f \)
• Visually, "best" is obvious

Define

• the error \( \e = \f - \u \)
• the (Eucledian) scalar product of two vectors: \( (\u,\v) \)
• the norm of \( \e \): \( ||\e|| = \sqrt{(\e, \e)} \)

The least squares method; principle

• Idea: find \( c_0 \) such that \( ||\e|| \) is minimized
• Actually, we always minimize \( E=||\e||^2 \)

 
$$ \begin{equation*} \frac{\partial E}{\partial c_0} = 0 \end{equation*} $$

The least squares method; calculations

 
$$ \begin{equation} E(c_0) = (\e,\e) = (\f,\f) - 2c_0(\f,\psib_0) + c_0^2(\psib_0,\psib_0) \end{equation} $$

 
$$ \begin{equation} \frac{\partial E}{\partial c_0} = -2(\f,\psib_0) + 2c_0 (\psib_0,\psib_0) = 0 \tag{2} \end{equation} $$

 
$$ \begin{equation} c_0 = \frac{(\f,\psib_0)}{(\psib_0,\psib_0)} \tag{3} \end{equation} $$

 
$$ \begin{equation} c_0 = \frac{3a + 5b}{a^2 + b^2} \end{equation} $$

Observation for later: the vanishing derivative (2) can be alternatively written as

 
$$ \begin{equation} (\e, \psib_0) = 0 \tag{4} \end{equation} $$

The projection (or Galerkin) method

• Backgrund: minimizing \( ||\e||^2 \) implies that \( \e \) is orthogonal to any vector \( \v \) in the space \( V \) (visually clear, but can easily be computed too)
• Alternative idea: demand \( (\e, \v) = 0,\quad\forall\v\in V \)
• Equivalent statement: \( (\e, \psib_0)=0 \) (see notes for why)
• Insert \( \e = \f - c_0\psib_0 \) and solve for \( c_0 \)
• Same equation for \( c_0 \) and hence same solution as in the least squares method

Approximation of general vectors

Given a vector \( \f \), find an approximation \( \u\in V \):

 
$$ \begin{equation*} V = \hbox{span}\,\{\psib_0,\ldots,\psib_N\} \end{equation*} $$

• We have a set of linearly independent basis vectors \( \psib_0,\ldots,\psib_N \)
• Any \( \u\in V \) can then be written as \( \u = \sum_{j=0}^Nc_j\psib_j \)

The least squares method

Idea: find \( c_0,\ldots,c_N \) such that \( E= ||\e||^2 \) is minimized, \( \e=\f-\u \).

 
$$ \begin{align*} E(c_0,\ldots,c_N) &= (\e,\e) = (\f -\sum_jc_j\psib_j,\f -\sum_jc_j\psib_j) \nonumber\\ &= (\f,\f) - 2\sum_{j=0}^Nc_j(\f,\psib_j) + \sum_{p=0}^N\sum_{q=0}^N c_pc_q(\psib_p,\psib_q) \end{align*} $$

 
$$ \begin{equation*} \frac{\partial E}{\partial c_i} = 0,\quad i=0,\ldots,N \end{equation*} $$

After some work we end up with a linear system

 
$$ \begin{align} \sum_{j=0}^N A_{i,j}c_j &= b_i,\quad i=0,\ldots,N\\ A_{i,j} &= (\psib_i,\psib_j)\\ b_i &= (\psib_i, \f) \end{align} $$

The projection (or Galerkin) method

Can be shown that minimizing \( ||\e|| \) implies that \( \e \) is orthogonal to all \( \v\in V \):

 
$$ (\e,\v)=0,\quad \forall\v\in V $$

 
which implies that \( \e \) most be orthogonal to each basis vector:

 
$$ \begin{equation} (\e,\psib_i)=0,\quad i=0,\ldots,N \tag{5} \end{equation} $$

This orthogonality condition is the principle of the projection (or Galerkin) method. Leads to the same linear system as in the least squares method.

Approximation of functions

Let \( V \) be a function space spanned by a set of basis functions \( \baspsi_0,\ldots,\baspsi_N \),

 
$$ \begin{equation*} V = \hbox{span}\,\{\baspsi_0,\ldots,\baspsi_N\} \end{equation*} $$

Find \( u\in V \) as a linear combination of the basis functions:

 
$$ \begin{equation} u = \sum_{j\in\If} c_j\baspsi_j,\quad\If = \{0,1,\ldots,N\} \tag{6} \end{equation} $$

The least squares method

• Extend the ideas from the vector case: minimize the (square) norm of the error.
• What norm? \( (f,g) = \int_\Omega f(x)g(x)\, dx \)

 
$$ \begin{equation} E = (e,e) = (f-u,f-u) = (f(x)-\sum_{j\in\If} c_j\baspsi_j(x), f(x)-\sum_{j\in\If} c_j\baspsi_j(x)) \tag{7} \end{equation} $$

 
$$ \begin{equation} E(c_0,\ldots,c_N) = (f,f) -2\sum_{j\in\If} c_j(f,\baspsi_i) + \sum_{p\in\If}\sum_{q\in\If} c_pc_q(\baspsi_p,\baspsi_q) \end{equation} $$

 
$$ \begin{equation*} \frac{\partial E}{\partial c_i} = 0,\quad i=\in\If \end{equation*} $$

After computations identical to the vector case, we get a linear system

 
$$ \begin{align} \sum_{j\in\If}^N A_{i,j}c_j &= b_i,\quad i\in\If \tag{8}\\ A_{i,j} &= (\baspsi_i,\baspsi_j) \tag{9}\\ b_i &= (f,\baspsi_i) \tag{10} \end{align} $$

The projection (or Galerkin) method

As before, minimizing \( (e,e) \) is equivalent to the projection (or Galerkin) method

 
$$ \begin{equation} (e,v)=0,\quad\forall v\in V \tag{11} \end{equation} $$

 
which means, as before,

 
$$ \begin{equation} (e,\baspsi_i)=0,\quad i\in\If \tag{12} \end{equation} $$

With the same algebra as in the multi-dimensional vector case, we get the same linear system as arose from the least squares method.

Example: linear approximation; problem

 
$$ \begin{equation*} V = \hbox{span}\,\{1, x\} \end{equation*} $$

 
That is, \( \baspsi_0(x)=1 \), \( \baspsi_1(x)=x \), and \( N=1 \). We seek

 
$$ \begin{equation*} u=c_0\baspsi_0(x) + c_1\baspsi_1(x) = c_0 + c_1x \end{equation*} $$

Example: linear approximation; solution

 
$$ \begin{align} A_{0,0} &= (\baspsi_0,\baspsi_0) = \int_1^21\cdot 1\, dx = 1\\ A_{0,1} &= (\baspsi_0,\baspsi_1) = \int_1^2 1\cdot x\, dx = 3/2\\ A_{1,0} &= A_{0,1} = 3/2\\ A_{1,1} &= (\baspsi_1,\baspsi_1) = \int_1^2 x\cdot x\,dx = 7/3 \end{align} $$

 
$$ \begin{align} b_1 &= (f,\baspsi_0) = \int_1^2 (10(x-1)^2 - 1)\cdot 1 \, dx = 7/3\\ b_2 &= (f,\baspsi_1) = \int_1^2 (10(x-1)^2 - 1)\cdot x\, dx = 13/3 \end{align} $$

 
Solution of 2x2 linear system:

 
$$ \begin{equation} c_0 = -38/3,\quad c_1 = 10,\quad u(x) = 10x - \frac{38}{3} \end{equation} $$

Example: linear approximation; plot

Implementation of the least squares method; ideas

Consider symbolic computation of the linear system, where

• \( f(x) \) is given as a sympy expression f (involving the symbol x),
• psi is a list of \( \sequencei{\baspsi} \),
• Omega is a 2-tuple/list holding the domain \( \Omega \)

Carry out the integrations, solve the linear system, and return \( u(x)=\sum_jc_j\baspsi_j(x) \)

Implementation of the least squares method; symbolic code

import sympy as sp

def least_squares(f, psi, Omega):
    N = len(psi) - 1
    A = sp.zeros((N+1, N+1))
    b = sp.zeros((N+1, 1))
    x = sp.Symbol('x')
    for i in range(N+1):
        for j in range(i, N+1):
            A[i,j] = sp.integrate(psi[i]*psi[j],
                                  (x, Omega[0], Omega[1]))
            A[j,i] = A[i,j]
        b[i,0] = sp.integrate(psi[i]*f, (x, Omega[0], Omega[1]))
    c = A.LUsolve(b)
    u = 0
    for i in range(len(psi)):
        u += c[i,0]*psi[i]
    return u, c

Observe: symmetric coefficient matrix so we can halve the integrations.

Implementation of the least squares method; numerical code

• Symbolic integration may be impossible and/or very slow
• Turn to pure numerical computations in those cases
• Supply Python functions f(x), psi(x,i), and a mesh x

def least_squares_numerical(f, psi, N, x,
                            integration_method='scipy',
                            orthogonal_basis=False):
    import scipy.integrate
    A = np.zeros((N+1, N+1))
    b = np.zeros(N+1)
    Omega = [x[0], x[-1]]
    dx = x[1] - x[0]

    for i in range(N+1):
        j_limit = i+1 if orthogonal_basis else N+1
        for j in range(i, j_limit):
            print '(%d,%d)' % (i, j)
            if integration_method == 'scipy':
                A_ij = scipy.integrate.quad(
                    lambda x: psi(x,i)*psi(x,j),
                    Omega[0], Omega[1], epsabs=1E-9, epsrel=1E-9)[0]
            elif ...
            A[i,j] = A[j,i] = A_ij

        if integration_method == 'scipy':
            b_i = scipy.integrate.quad(
                lambda x: f(x)*psi(x,i), Omega[0], Omega[1],
                epsabs=1E-9, epsrel=1E-9)[0]
        elif ...
        b[i] = b_i

    c = b/np.diag(A) if orthogonal_basis else np.linalg.solve(A, b)
    u = sum(c[i]*psi(x, i) for i in range(N+1))
    return u, c

Implementation of the least squares method; plotting

Compare \( f \) and \( u \) visually:

def comparison_plot(f, u, Omega, filename='tmp.pdf'):
    x = sp.Symbol('x')
    # Turn f and u to ordinary Python functions
    f = sp.lambdify([x], f, modules="numpy")
    u = sp.lambdify([x], u, modules="numpy")
    resolution = 401  # no of points in plot
    xcoor  = linspace(Omega[0], Omega[1], resolution)
    exact  = f(xcoor)
    approx = u(xcoor)
    plot(xcoor, approx)
    hold('on')
    plot(xcoor, exact)
    legend(['approximation', 'exact'])
    savefig(filename)

All code in module approx1D.py

Implementation of the least squares method; application

>>> from approx1D import *
>>> x = sp.Symbol('x')
>>> f = 10*(x-1)**2-1
>>> u, c = least_squares(f=f, psi=[1, x], Omega=[1, 2])
>>> comparison_plot(f, u, Omega=[1, 2])

Perfect approximation; parabola approximating parabola

• What if we add \( \baspsi_2=x^2 \) to the space \( V \)?
• That is, approximating a parabola by any parabola?
• (Hopefully we get the exact parabola!)

>>> from approx1D import *
>>> x = sp.Symbol('x')
>>> f = 10*(x-1)**2-1
>>> u, c = least_squares(f=f, psi=[1, x, x**2], Omega=[1, 2])
>>> print u
10*x**2 - 20*x + 9
>>> print sp.expand(f)
10*x**2 - 20*x + 9

Perfect approximation; the general result

• What if we use \( \psi_i(x)=x^i \) for \( i=0,\ldots,N=40 \)?
• The output from least_squares is \( c_i=0 \) for \( i>2 \)

Perfect approximation; proof of the general result

If \( f\in V \), \( f=\sum_{j\in\If}d_j\baspsi_j \), for some \( \sequencei{d} \). Then

 
$$ \begin{equation*} b_i = (f,\baspsi_i) = \sum_{j\in\If}d_j(\baspsi_j, \baspsi_i) = \sum_{j\in\If} d_jA_{i,j} \end{equation*} $$

 
The linear system \( \sum_j A_{i,j}c_j = b_i \), \( i\in\If \), is then

 
$$ \begin{equation*} \sum_{j\in\If}c_jA_{i,j} = \sum_{j\in\If}d_jA_{i,j},\quad i\in\If \end{equation*} $$

 
which implies that \( c_i=d_i \) for \( i\in\If \) and \( u \) is identical to \( f \).

Finite-precision/numerical computations

The previous computations were symbolic. What if we solve the linear system numerically with standard arrays?

• Column 2: sympy.mpmath.fp.matrix and sympy.mpmath.fp.lu_solve
• Column 3: numpy arrays with numpy.float32 entries
• Column 4: numpy arrays with numpy.float64 entries

Ill-conditioning (1)

Observations:

• Significant round-off errors in the numerical computations (!)
• But if we plot the approximations they look good (!)

Problem: The basis functions \( x^i \) become almost linearly dependent for large \( N \).

Ill-conditioning (2)

• Almost linearly dependent basis functions give almost singular matrices
• Such matrices are said to be ill conditioned, and Gaussian elimination is severely affected by round-off errors
• The basis \( 1, x, x^2, x^3, x^4, \ldots \) is a bad basis
• Polynomials are fine as basis, but the more orthogonal they are, \( (\baspsi_i,\baspsi_j)\approx 0 \), the better

Fourier series approximation; problem and code

Consider

 
$$ \begin{equation*} V = \hbox{span}\,\{ \sin \pi x, \sin 2\pi x,\ldots,\sin (N+1)\pi x\} \end{equation*} $$

N = 3
from sympy import sin, pi
psi = [sin(pi*(i+1)*x) for i in range(N+1)]
f = 10*(x-1)**2 - 1
Omega = [0, 1]
u, c = least_squares(f, psi, Omega)
comparison_plot(f, u, Omega)

Fourier series approximation; plot

\( N=3 \) vs \( N=11 \):

Fourier series approximation; improvements

• Considerably improvement by \( N=11 \)
• But always discrepancy of \( f(0)-u(0)=9 \) at \( x=0 \), because all the \( \baspsi_i(0)=0 \) and hence \( u(0)=0 \)
• Possible remedy: add a term that leads to correct boundary values

 
$$ \begin{equation} u(x) = f(0)(1-x) + xf(1) + \sum_{j\in\If} c_j\baspsi_j(x) \end{equation} $$

 
The extra term ensures \( u(0)=f(0) \) and \( u(1)=f(1) \) and is a strikingly good help to get a good approximation!

Fourier series approximation; final results

\( N=3 \) vs \( N=11 \):

Orthogonal basis functions

This choice of sine functions as basis functions is popular because

• the basis functions are orthogonal: \( (\baspsi_i,\baspsi_j)=0 \)
• implying that \( A_{i,j} \) is a diagonal matrix
• implying that we can solve for \( c_i = 2\int_0^1 f(x)\sin ((i+1)\pi x) dx \)

In general for an orthogonal basis, \( A_{i,j} \) is diagonal and we can easily solve for \( c_i \):

 
$$ c_i = \frac{b_i}{A_{i,i}} = \frac{(f,\baspsi_i)}{(\baspsi_i,\baspsi_i)} $$

The collocation or interpolation method; ideas and math

Here is another idea for approximating \( f(x) \) by \( u(x)=\sum_jc_j\baspsi_j \):

• Force \( u(\xno{i}) = f(\xno{i}) \) at some selected collocation points \( \sequencei{x} \)
• Then \( u \) interpolates \( f \)
• The method is known as interpolation or collocation

 
$$ \begin{equation} u(\xno{i}) = \sum_{j\in\If} c_j \baspsi_j(\xno{i}) = f(\xno{i}) \quad i\in\If,N \end{equation} $$

This is a linear system with no need for integration:

 
$$ \begin{align} \sum_{j\in\If} A_{i,j}c_j &= b_i,\quad i\in\If\\ A_{i,j} &= \baspsi_j(\xno{i})\\ b_i &= f(\xno{i}) \end{align} $$

No symmetric matrix: \( \baspsi_j(\xno{i})\neq \baspsi_i(\xno{j}) \) in general

The collocation or interpolation method; implementation

points holds the interpolation/collocation points

def interpolation(f, psi, points):
    N = len(psi) - 1
    A = sp.zeros((N+1, N+1))
    b = sp.zeros((N+1, 1))
    x = sp.Symbol('x')
    # Turn psi and f into Python functions
    psi = [sp.lambdify([x], psi[i]) for i in range(N+1)]
    f = sp.lambdify([x], f)
    for i in range(N+1):
        for j in range(N+1):
            A[i,j] = psi[j](points[i])
        b[i,0] = f(points[i])
    c = A.LUsolve(b)
    u = 0
    for i in range(len(psi)):
        u += c[i,0]*psi[i](x)
    return u

The collocation or interpolation method; approximating a parabola by linear functions

• Potential difficulty: how to choose \( \xno{i} \)?
• The results are sensitive to the points!

\( (4/3,5/3) \) vs \( (1,2) \):

Lagrange polynomials; motivation and ideas

Motivation:

• The interpolation/collocation method avoids integration
• With a diagonal matrix \( A_{i,j} = \baspsi_j(\xno{i}) \) we can solve the linear system by hand

The Lagrange interpolating polynomials \( \baspsi_j \) have the property that

 
$$ \baspsi_i(\xno{j}) =\delta_{ij},\quad \delta_{ij} = \left\lbrace\begin{array}{ll} 1, & i=j\\ 0, & i\neq j \end{array}\right. $$

Hence, \( c_i = f(x_i) \) and

 
$$ \begin{equation} u(x) = \sum_{j\in\If} f(\xno{i})\baspsi_i(x) \end{equation} $$

• Lagrange polynomials and interpolation/collocation look convenient
• Lagrange polynomials are very much used in the finite element method

Lagrange polynomials; formula and code

 
$$ \begin{equation} \baspsi_i(x) = \prod_{j=0,j\neq i}^N \frac{x-\xno{j}}{\xno{i}-\xno{j}} = \frac{x-x_0}{\xno{i}-x_0}\cdots\frac{x-\xno{i-1}}{\xno{i}-\xno{i-1}}\frac{x-\xno{i+1}}{\xno{i}-\xno{i+1}} \cdots\frac{x-x_N}{\xno{i}-x_N} \tag{13} \end{equation} $$

def Lagrange_polynomial(x, i, points):
    p = 1
    for k in range(len(points)):
        if k != i:
            p *= (x - points[k])/(points[i] - points[k])
    return p

Lagrange polynomials; successful example

Lagrange polynomials; a less successful example

Lagrange polynomials; oscillatory behavior

12 points, degree 11, plot of two of the Lagrange polynomials - note that they are zero at all points except one.

Problem: strong oscillations near the boundaries for larger \( N \) values.

Lagrange polynomials; remedy for strong oscillations

The oscillations can be reduced by a more clever choice of interpolation points, called the Chebyshev nodes:

 
$$ \begin{equation} \xno{i} = \half (a+b) + \half(b-a)\cos\left( \frac{2i+1}{2(N+1)}pi\right),\quad i=0\ldots,N \end{equation} $$

 
on an interval \( [a,b] \).

Lagrange polynomials; recalculation with Chebyshev nodes

Lagrange polynomials; less oscillations with Chebyshev nodes

12 points, degree 11, plot of two of the Lagrange polynomials - note that they are zero at all points except one.

Finite element basis functions

The basis functions have so far been global: \( \baspsi_i(x) \neq 0 \) almost everywhere

In the finite element method we use basis functions with local support

• Local support: \( \baspsi_i(x) \neq 0 \) for \( x \) in a small subdomain of \( \Omega \)
• Typically hat-shaped
• \( u(x) \) based on these \( \baspsi_i \) is a piecewise polynomial defined over many (small) subdomains
• We introduce \( \basphi_i \) as the name of these finite element hat functions (and for now choose \( \baspsi_i=\basphi_i \))

The linear combination of hat functions is a piecewise linear function

Elements and nodes

Split \( \Omega \) into non-overlapping subdomains called elements:

 
$$ \begin{equation} \Omega = \Omega^{(0)}\cup \cdots \cup \Omega^{(N_e)} \end{equation} $$

On each element, introduce points called nodes: \( \xno{0},\ldots,\xno{N_n} \)

• The finite element basis functions are named \( \basphi_i(x) \)
• \( \basphi_i=1 \) at node \( i \) and 0 at all other nodes
• \( \basphi_i \) is a Lagrange polynomial on each element
• For nodes at the boundary between two elements, \( \basphi_i \) is made up of a Lagrange polynomial over each element

Example on elements with two nodes (P1 elements)

Data structure: nodes holds coordinates or nodes, elements holds the node numbers in each element

nodes = [0, 1.2, 2.4, 3.6, 4.8, 5]
elements = [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5]]

Illustration of two basis functions on the mesh

Example on elements with three nodes (P2 elements)

nodes = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]
elements = [[0, 1, 2], [2, 3, 4], [4, 5, 6], [6, 7, 8]]

Some corresponding basis functions (P2 elements)

Examples on elements with four nodes per element (P3 elements)

d = 3  # d+1 nodes per element
num_elements = 4
num_nodes = num_elements*d + 1
nodes = [i*0.5 for i in range(num_nodes)]
elements = [[i*d+j for j in range(d+1)] for i in range(num_elements)]

Some corresponding basis functions (P3 elements)

The numbering does not need to be regular from left to right

nodes = [1.5, 5.5, 4.2, 0.3, 2.2, 3.1]
elements = [[2, 1], [4, 5], [0, 4], [3, 0], [5, 2]]

Interpretation of the coefficients \( c_i \)

Important property: \( c_i \) is the value of \( u \) at node \( i \), \( \xno{i} \):

 
$$ \begin{equation} u(\xno{i}) = \sum_{j\in\If} c_j\basphi_j(\xno{i}) = c_i\basphi_i(\xno{i}) = c_i \tag{14} \end{equation} $$

because \( \basphi_j(\xno{i}) =0 \) if \( i\neq j \)

Properties of the basis functions

• \( \basphi_i(x) \neq 0 \) only on those elements that contain global node \( i \)
• \( \basphi_i(x)\basphi_j(x) \neq 0 \) if and only if \( i \) and \( j \) are global node numbers in the same element

Since \( A_{i,j}=\int\basphi_i\basphi_j\dx \), most of the elements in the coefficient matrix will be zero

How to construct quadratic \( \basphi_i \) (P2 elements)

1. Associate Lagrange polynomials with the nodes in an element
2. When the polynomial is 1 on the element boundary, combine it with the polynomial in the neighboring element

Example on linear \( \basphi_i \) (P1 elements)

 
$$ \begin{equation} \basphi_i(x) = \left\lbrace\begin{array}{ll} 0, & x < \xno{i-1}\\ (x - \xno{i-1})/h & \xno{i-1} \leq x < \xno{i}\\ 1 - (x - x_{i})/h, & \xno{i} \leq x < \xno{i+1}\\ 0, & x\geq \xno{i+1} \end{array} \right. \tag{15} \end{equation} $$

Example on cubic \( \basphi_i \) (P3 elements)

Calculating the linear system for \( c_i \)

Computing a specific matrix entry (1)

\( A_{2,3}=\int_\Omega\basphi_2\basphi_3 dx \): \( \basphi_2\basphi_3\neq 0 \) only over element 2. There,

 
$$ \basphi_3(x) = (x-x_2)/h,\quad \basphi_2(x) = 1- (x-x_2)/h$$

 
$$ A_{2,3} = \int_\Omega \basphi_2\basphi_{3}\dx = \int_{\xno{2}}^{\xno{3}} \left(1 - \frac{x - \xno{2}}{h}\right) \frac{x - x_{2}}{h} \dx = \frac{h}{6} $$

Computing a specific matrix entry (2)

 
$$ A_{2,2} = \int_{\xno{1}}^{\xno{2}} \left(\frac{x - \xno{1}}{h}\right)^2\dx + \int_{\xno{2}}^{\xno{3}} \left(1 - \frac{x - \xno{2}}{h}\right)^2\dx = \frac{h}{3} $$

Calculating a general row in the matrix; figure

 
$$ A_{i,i-1} = \int_\Omega \basphi_i\basphi_{i-1}\dx = \hbox{?}$$

Calculating a general row in the matrix; details

 
$$ \begin{align*} A_{i,i-1} &= \int_\Omega \basphi_i\basphi_{i-1}\dx\\ &= \underbrace{\int_{\xno{i-2}}^{\xno{i-1}} \basphi_i\basphi_{i-1}\dx}_{\basphi_i=0} + \int_{\xno{i-1}}^{\xno{i}} \basphi_i\basphi_{i-1}\dx + \underbrace{\int_{\xno{i}}^{\xno{i+1}} \basphi_i\basphi_{i-1}\dx}_{\basphi_{i-1}=0}\\ &= \int_{\xno{i-1}}^{\xno{i}} \underbrace{\left(\frac{x - x_{i}}{h}\right)}_{\basphi_i(x)} \underbrace{\left(1 - \frac{x - \xno{i-1}}{h}\right)}_{\basphi_{i-1}(x)} \dx = \frac{h}{6} \end{align*} $$

• \( A_{i,i+1}=A_{i,i-1} \) due to symmetry
• \( A_{i,i}=h/3 \) (same calculation as for \( A_{2,2} \))
• \( A_{0,0}=A_{N,N}=h/3 \) (only one element)

Calculation of the right-hand side

 
$$ \begin{equation} b_i = \int_\Omega\basphi_i(x)f(x)\dx = \int_{\xno{i-1}}^{\xno{i}} \frac{x - \xno{i-1}}{h} f(x)\dx + \int_{x_{i}}^{\xno{i+1}} \left(1 - \frac{x - x_{i}}{h}\right) f(x) \dx \tag{16} \end{equation} $$

Need a specific \( f(x) \) to do more...

Specific example with two elements; linear system and solution

• \( f(x)=x(1-x) \) on \( \Omega=[0,1] \)
• Two equal-sized elements \( [0,0.5] \) and \( [0.5,1] \)

 
$$ \begin{equation*} A = \frac{h}{6}\left(\begin{array}{ccc} 2 & 1 & 0\\ 1 & 4 & 1\\ 0 & 1 & 2 \end{array}\right),\quad b = \frac{h^2}{12}\left(\begin{array}{c} 2 - 3h\\ 12 - 14h\\ 10 -17h \end{array}\right) \end{equation*} $$

 
$$ \begin{equation*} c_0 = \frac{h^2}{6},\quad c_1 = h - \frac{5}{6}h^2,\quad c_2 = 2h - \frac{23}{6}h^2 \end{equation*} $$

Specific example with two elements; plot

 
$$ \begin{equation*} u(x)=c_0\basphi_0(x) + c_1\basphi_1(x) + c_2\basphi_2(x)\end{equation*} $$

Specific example: what about four elements?

Assembly of elementwise computations

Split the integrals into elementwise integrals

 
$$ \begin{equation} A_{i,j} = \int_\Omega\basphi_i\basphi_jdx = \sum_{e} \int_{\Omega^{(e)}} \basphi_i\basphi_jdx,\quad A^{(e)}_{i,j}=\int_{\Omega^{(e)}} \basphi_i\basphi_jdx \tag{17} \end{equation} $$

Important:

• \( A^{(e)}_{i,j}\neq 0 \) if and only if \( i \) and \( j \) are nodes in element \( e \) (otherwise no overlap between the basis functions)
• all the nonzero elements in \( A^{(e)}_{i,j} \) are collected in an element matrix

The element matrix

 
$$ \tilde A^{(e)} = \{ \tilde A^{(e)}_{r,s}\},\quad \tilde A^{(e)}_{r,s} = \int_{\Omega^{(e)}}\basphi_{q(e,r)}\basphi_{q(e,s)}dx, \quad r,s\in\Ifd=\{0,\ldots,d\} $$

• \( r,s \) run over local node numbers in an element; \( i,j \) run over global node numbers
• \( i=q(e,r) \): mapping of local node number \( r \) in element \( e \) to the global node number \( i \) (math equivalent to i=elements[e][r])
• Add \( \tilde A^{(e)}_{r,s} \) into the global \( A_{i,j} \) (assembly)

 
$$ \begin{equation} A_{q(e,r),q(e,s)} := A_{q(e,r),q(e,s)} + \tilde A^{(e)}_{r,s},\quad r,s\in\Ifd \end{equation} $$

Illustration of the matrix assembly: regularly numbered P1 elements

Animation

Illustration of the matrix assembly: regularly numbered P3 elements

Animation

Illustration of the matrix assembly: irregularly numbered P1 elements

Animation

Assembly of the right-hand side

 
$$ \begin{equation} b_i = \int_\Omega f(x)\basphi_i(x)dx = \sum_{e} \int_{\Omega^{(e)}} f(x)\basphi_i(x)dx,\quad b^{(e)}_{i}=\int_{\Omega^{(e)}} f(x)\basphi_i(x)dx \end{equation} $$

Important:

• \( b_i^{(e)}\neq 0 \) if and only if global node \( i \) is a node in element \( e \) (otherwise \( \basphi_i=0 \))
• The \( d+1 \) nonzero \( b_i^{(e)} \) can be collected in an element vector \( \tilde b_r^{(e)}=\{ \tilde b_r^{(e)}\} \), \( r\in\Ifd \)

Assembly:

 
$$ \begin{equation} b_{q(e,r)} := b_{q(e,r)} + \tilde b^{(e)}_{r},\quad r,s\in\Ifd \end{equation} $$

Mapping to a reference element

Instead of computing

 
$$ \begin{equation*} \tilde A^{(e)}_{r,s} = \int_{\Omega^{(e)}}\basphi_{q(e,r)}(x)\basphi_{q(e,s)}(x)dx = \int_{x_L}^{x_R}\basphi_{q(e,r)}(x)\basphi_{q(e,s)}(x)dx \end{equation*} $$

 
we now map \( [x_L, x_R] \) to a standardized reference element domain \( [-1,1] \) with local coordinate \( X \)

Affine mapping

 
$$ \begin{equation} x = \half (x_L + x_R) + \half (x_R - x_L)X \tag{18} \end{equation} $$

 
or rewritten as

 
$$ \begin{equation} x = x_m + {\half}hX, \qquad x_m=(x_L+x_R)/2 \tag{19} \end{equation} $$

Integral transformation

Reference element integration: just change integration variable from \( x \) to \( X \). Introduce local basis function

 
$$ \begin{equation} \refphi_r(X) = \basphi_{q(e,r)}(x(X)) \end{equation} $$

 
$$ \begin{equation} \tilde A^{(e)}_{r,s} = \int_{\Omega^{(e)}}\basphi_{q(e,r)}(x)\basphi_{q(e,s)}(x)dx = \int\limits_{-1}^1 \refphi_r(X)\refphi_s(X)\underbrace{\frac{dx}{dX}}_{\det J = h/2}dX = \int\limits_{-1}^1 \refphi_r(X)\refphi_s(X)\det J\,dX \end{equation} $$

 
$$ \begin{equation} \tilde b^{(e)}_{r} = \int_{\Omega^{(e)}}f(x)\basphi_{q(e,r)}(x)dx = \int\limits_{-1}^1 f(x(X))\refphi_r(X)\det J\,dX \tag{20} \end{equation} $$

Advantages of the reference element

• Always the same domain for integration: \( [-1,1] \)
• We only need formulas for \( \refphi_r(X) \) over one element (no piecewise polynomial definition)
• \( \refphi_r(X) \) is the same for all elements: no dependence on element length and location, which is "factored out" in the mapping and \( \det J \)

Standardized basis functions for P1 elements

 
$$ \begin{align} \refphi_0(X) &= \half (1 - X) \tag{21}\\ \refphi_1(X) &= \half (1 + X) \tag{22} \end{align} $$

Standardized basis functions for P2 elements

P2 elements:

 
$$ \begin{align} \refphi_0(X) &= \half (X-1)X\\ \refphi_1(X) &= 1 - X^2\\ \refphi_2(X) &= \half (X+1)X \end{align} $$

Easy to generalize to arbitrary order!

Integration over a reference element; element matrix

P1 elements and \( f(x)=x(1-x) \).

 
$$ \begin{align} \tilde A^{(e)}_{0,0} &= \int_{-1}^1 \refphi_0(X)\refphi_0(X)\frac{h}{2} dX\nonumber\\ &=\int_{-1}^1 \half(1-X)\half(1-X) \frac{h}{2} dX = \frac{h}{8}\int_{-1}^1 (1-X)^2 dX = \frac{h}{3} \tag{23}\\ \tilde A^{(e)}_{1,0} &= \int_{-1}^1 \refphi_1(X)\refphi_0(X)\frac{h}{2} dX\nonumber\\ &=\int_{-1}^1 \half(1+X)\half(1-X) \frac{h}{2} dX = \frac{h}{8}\int_{-1}^1 (1-X^2) dX = \frac{h}{6}\\ \tilde A^{(e)}_{0,1} &= \tilde A^{(e)}_{1,0} \tag{24}\\ \tilde A^{(e)}_{1,1} &= \int_{-1}^1 \refphi_1(X)\refphi_1(X)\frac{h}{2} dX\nonumber\\ &=\int_{-1}^1 \half(1+X)\half(1+X) \frac{h}{2} dX = \frac{h}{8}\int_{-1}^1 (1+X)^2 dX = \frac{h}{3} \tag{25} \end{align} $$

Integration over a reference element; element vector

 
$$ \begin{align} \tilde b^{(e)}_{0} &= \int_{-1}^1 f(x(X))\refphi_0(X)\frac{h}{2} dX\nonumber\\ &= \int_{-1}^1 (x_m + \half hX)(1-(x_m + \half hX)) \half(1-X)\frac{h}{2} dX \nonumber\\ &= - \frac{1}{24} h^{3} + \frac{1}{6} h^{2} x_{m} - \frac{1}{12} h^{2} - \half h x_{m}^{2} + \half h x_{m} \tag{26}\\ \tilde b^{(e)}_{1} &= \int_{-1}^1 f(x(X))\refphi_1(X)\frac{h}{2} dX\nonumber\\ &= \int_{-1}^1 (x_m + \half hX)(1-(x_m + \half hX)) \half(1+X)\frac{h}{2} dX \nonumber\\ &= - \frac{1}{24} h^{3} - \frac{1}{6} h^{2} x_{m} + \frac{1}{12} h^{2} - \half h x_{m}^{2} + \half h x_{m} \end{align} $$

\( x_m \): element midpoint.

Tedious calculations! Let's use symbolic software

>>> import sympy as sp
>>> x, x_m, h, X = sp.symbols('x x_m h X')
>>> sp.integrate(h/8*(1-X)**2, (X, -1, 1))
h/3
>>> sp.integrate(h/8*(1+X)*(1-X), (X, -1, 1))
h/6
>>> x = x_m + h/2*X
>>> b_0 = sp.integrate(h/4*x*(1-x)*(1-X), (X, -1, 1))
>>> print b_0
-h**3/24 + h**2*x_m/6 - h**2/12 - h*x_m**2/2 + h*x_m/2

Can printe out in LaTeX too (convenient for copying into reports):

>>> print sp.latex(b_0, mode='plain')
- \frac{1}{24} h^{3} + \frac{1}{6} h^{2} x_{m}
- \frac{1}{12} h^{2} - \half h x_{m}^{2}
+ \half h x_{m}

Implementation

• Coming functions appear in fe_approx1D.py
• Functions can operate in symbolic or numeric mode
• The code documents all steps in finite element calculations!

Compute finite element basis functions in the reference element

Let \( \refphi_r(X) \) be a Lagrange polynomial of degree d:

import sympy as sp
import numpy as np

def phi_r(r, X, d):
    if isinstance(X, sp.Symbol):
        h = sp.Rational(1, d)  # node spacing
        nodes = [2*i*h - 1 for i in range(d+1)]
    else:
        # assume X is numeric: use floats for nodes
        nodes = np.linspace(-1, 1, d+1)
    return Lagrange_polynomial(X, r, nodes)

def Lagrange_polynomial(x, i, points):
    p = 1
    for k in range(len(points)):
        if k != i:
            p *= (x - points[k])/(points[i] - points[k])
    return p

def basis(d=1):
    """Return the complete basis."""
    X = sp.Symbol('X')
    phi = [phi_r(r, X, d) for r in range(d+1)]
    return phi

Compute the element matrix

def element_matrix(phi, Omega_e, symbolic=True):
    n = len(phi)
    A_e = sp.zeros((n, n))
    X = sp.Symbol('X')
    if symbolic:
        h = sp.Symbol('h')
    else:
        h = Omega_e[1] - Omega_e[0]
    detJ = h/2  # dx/dX
    for r in range(n):
        for s in range(r, n):
            A_e[r,s] = sp.integrate(phi[r]*phi[s]*detJ, (X, -1, 1))
            A_e[s,r] = A_e[r,s]
    return A_e

Example on symbolic vs numeric element matrix

>>> from fe_approx1D import *
>>> phi = basis(d=1)
>>> phi
[1/2 - X/2, 1/2 + X/2]
>>> element_matrix(phi, Omega_e=[0.1, 0.2], symbolic=True)
[h/3, h/6]
[h/6, h/3]
>>> element_matrix(phi, Omega_e=[0.1, 0.2], symbolic=False)
[0.0333333333333333, 0.0166666666666667]
[0.0166666666666667, 0.0333333333333333]

Compute the element vector

def element_vector(f, phi, Omega_e, symbolic=True):
    n = len(phi)
    b_e = sp.zeros((n, 1))
    # Make f a function of X
    X = sp.Symbol('X')
    if symbolic:
        h = sp.Symbol('h')
    else:
        h = Omega_e[1] - Omega_e[0]
    x = (Omega_e[0] + Omega_e[1])/2 + h/2*X  # mapping
    f = f.subs('x', x)  # substitute mapping formula for x
    detJ = h/2  # dx/dX
    for r in range(n):
        b_e[r] = sp.integrate(f*phi[r]*detJ, (X, -1, 1))
    return b_e

Note f.subs('x', x): replace x by \( x(X) \) such that f contains X

Fallback on numerical integration if symbolic integration fails

• Element matrix: only polynomials and sympy always succeeds
• Element vector: \( \int f\refphi \dx \) can fail (sympy then returns an Integral object instead of a number)

def element_vector(f, phi, Omega_e, symbolic=True):
        ...
        I = sp.integrate(f*phi[r]*detJ, (X, -1, 1))  # try...
        if isinstance(I, sp.Integral):
            h = Omega_e[1] - Omega_e[0]  # Ensure h is numerical
            detJ = h/2
            integrand = sp.lambdify([X], f*phi[r]*detJ)
            I = sp.mpmath.quad(integrand, [-1, 1])
        b_e[r] = I
        ...

Linear system assembly and solution

def assemble(nodes, elements, phi, f, symbolic=True):
    N_n, N_e = len(nodes), len(elements)
    zeros = sp.zeros if symbolic else np.zeros
    A = zeros((N_n, N_n))
    b = zeros((N_n, 1))
    for e in range(N_e):
        Omega_e = [nodes[elements[e][0]], nodes[elements[e][-1]]]

        A_e = element_matrix(phi, Omega_e, symbolic)
        b_e = element_vector(f, phi, Omega_e, symbolic)

        for r in range(len(elements[e])):
            for s in range(len(elements[e])):
                A[elements[e][r],elements[e][s]] += A_e[r,s]
            b[elements[e][r]] += b_e[r]
    return A, b

Linear system solution

if symbolic:
    c = A.LUsolve(b)           # sympy arrays, symbolic Gaussian elim.
else:
    c = np.linalg.solve(A, b)  # numpy arrays, numerical solve

Note: the symbolic computation of A and b and the symbolic solution can be very tedious.

Example on computing symbolic approximations

>>> h, x = sp.symbols('h x')
>>> nodes = [0, h, 2*h]
>>> elements = [[0, 1], [1, 2]]
>>> phi = basis(d=1)
>>> f = x*(1-x)
>>> A, b = assemble(nodes, elements, phi, f, symbolic=True)
>>> A
[h/3,   h/6,   0]
[h/6, 2*h/3, h/6]
[  0,   h/6, h/3]
>>> b
[     h**2/6 - h**3/12]
[      h**2 - 7*h**3/6]
[5*h**2/6 - 17*h**3/12]
>>> c = A.LUsolve(b)
>>> c
[                           h**2/6]
[12*(7*h**2/12 - 35*h**3/72)/(7*h)]
[  7*(4*h**2/7 - 23*h**3/21)/(2*h)]

Example on computing numerical approximations

>>> nodes = [0, 0.5, 1]
>>> elements = [[0, 1], [1, 2]]
>>> phi = basis(d=1)
>>> x = sp.Symbol('x')
>>> f = x*(1-x)
>>> A, b = assemble(nodes, elements, phi, f, symbolic=False)
>>> A
[ 0.166666666666667, 0.0833333333333333,                  0]
[0.0833333333333333,  0.333333333333333, 0.0833333333333333]
[                 0, 0.0833333333333333,  0.166666666666667]
>>> b
[          0.03125]
[0.104166666666667]
[          0.03125]
>>> c = A.LUsolve(b)
>>> c
[0.0416666666666666]
[ 0.291666666666667]
[0.0416666666666666]

The structure of the coefficient matrix

>>> d=1; N_e=8; Omega=[0,1]  # 8 linear elements on [0,1]
>>> phi = basis(d)
>>> f = x*(1-x)
>>> nodes, elements = mesh_symbolic(N_e, d, Omega)
>>> A, b = assemble(nodes, elements, phi, f, symbolic=True)
>>> A
[h/3,   h/6,     0,     0,     0,     0,     0,     0,   0]
[h/6, 2*h/3,   h/6,     0,     0,     0,     0,     0,   0]
[  0,   h/6, 2*h/3,   h/6,     0,     0,     0,     0,   0]
[  0,     0,   h/6, 2*h/3,   h/6,     0,     0,     0,   0]
[  0,     0,     0,   h/6, 2*h/3,   h/6,     0,     0,   0]
[  0,     0,     0,     0,   h/6, 2*h/3,   h/6,     0,   0]
[  0,     0,     0,     0,     0,   h/6, 2*h/3,   h/6,   0]
[  0,     0,     0,     0,     0,     0,   h/6, 2*h/3, h/6]
[  0,     0,     0,     0,     0,     0,     0,   h/6, h/3]

Note: do this by hand to understand what is going on!

General result: the coefficient matrix is sparse

• Sparse = most of the entries are zeros
• Below: P1 elements

 
$$ \begin{equation} A = \frac{h}{6} \left( \begin{array}{cccccccccc} 2 & 1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ 1 & 4 & 1 & \ddots & & & & & \vdots \\ 0 & 1 & 4 & 1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & 1 & 4 & 1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & 1 & 4 & 1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 1 & 2 \end{array} \right) \end{equation} $$

Exemplifying the sparsity for P2 elements

 
$$ \begin{equation} A = \frac{h}{30} \left( \begin{array}{ccccccccc} 4 & 2 & - 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 2 & 16 & 2 & 0 & 0 & 0 & 0 & 0 & 0\\- 1 & 2 & 8 & 2 & - 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 16 & 2 & 0 & 0 & 0 & 0\\ 0 & 0 & - 1 & 2 & 8 & 2 & - 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 2 & 16 & 2 & 0 & 0\\ 0 & 0 & 0 & 0 & - 1 & 2 & 8 & 2 & - 1 \\0 & 0 & 0 & 0 & 0 & 0 & 2 & 16 & 2\\0 & 0 & 0 & 0 & 0 & 0 & - 1 & 2 & 4 \end{array} \right) \end{equation} $$

Matrix sparsity pattern for regular/random numbering of P1 elements

• Left: number nodes and elements from left to right
• Right: number nodes and elements arbitrarily

Matrix sparsity pattern for regular/random numbering of P3 elements

• Left: number nodes and elements from left to right
• Right: number nodes and elements arbitrarily

Sparse matrix storage and solution

The minimum storage requirements for the coefficient matrix \( A_{i,j} \):

• P1 elements: only 3 nonzero entires per row
• P2 elements: only 5 nonzero entires per row
• P3 elements: only 7 nonzero entires per row
• It is important to utilize sparse storage and sparse solvers
• In Python: scipy.sparse package

Approximate \( f\sim x^9 \) by various elements; code

Compute a mesh with \( N_e \) elements, basis functions of degree \( d \), and approximate a given symbolic expression \( f(x) \) by a finite element expansion \( u(x) = \sum_jc_j\basphi_j(x) \):

import sympy as sp
from fe_approx1D import approximate
x = sp.Symbol('x')

approximate(f=x*(1-x)**8, symbolic=False, d=1, N_e=4)
approximate(f=x*(1-x)**8, symbolic=False, d=2, N_e=2)
approximate(f=x*(1-x)**8, symbolic=False, d=1, N_e=8)
approximate(f=x*(1-x)**8, symbolic=False, d=2, N_e=4)

Approximate \( f\sim x^9 \) by various elements; plot

Comparison of finite element and finite difference approximation

• Finite difference approximation of a function \( f(x) \): simply choose \( u_i = f(x_i) \) (interpolation)
• Galerkin/projection and least squares method: must derive and solve a linear system
• What is really the difference in \( u \)?

Interpolation/collocation with finite elements

Let \( \{\xno{i}\}_{i\in\If} \) be the nodes in the mesh. Collocation means

 
$$ \begin{equation} u(\xno{i})=f(\xno{i}),\quad i\in\If, \end{equation} $$

 
which translates to

 
$$ \sum_{j\in\If} c_j \basphi_j(\xno{i}) = f(\xno{i}),$$

 
but \( \basphi_j(\xno{i})=0 \) if \( i\neq j \) so the sum collapses to one term \( c_i\basphi_i(\xno{i}) = c_i \), and we have the result

 
$$ \begin{equation} c_i = f(\xno{i}) \end{equation} $$

Same result as the standard finite difference approach, but finite elements define \( u \) also between the \( \xno{i} \) points

Galerkin/project and least squares vs collocation/interpolation or finite differences

• Scope: work with P1 elements
• Use projection/Galerkin or least squares (equivalent)
• Interpret the resulting linear system as finite difference equations

The P1 finite element machinery results in a linear system where equation no \( i \) is

 
$$ \begin{equation} \frac{h}{6}(u_{i-1} + 4u_i + u_{i+1}) = (f,\basphi_i) \tag{27} \end{equation} $$

Note:

• We have used \( u_i \) for \( c_i \) to make notation similar to finite differences
• The finite difference counterpart is just \( u_i=f_i \)

Expressing the left-hand side in finite difference operator notation

Rewrite the left-hand side of finite element equation no \( i \):

 
$$ \begin{equation} h(u_i + \frac{1}{6}(u_{i-1} - 2u_i + u_{i+1})) = [h(u + \frac{h^2}{6}D_x D_x u)]_i \end{equation} $$

 
This is the standard finite difference approximation of

 
$$ h(u + \frac{h^2}{6}u'')$$

Treating the right-hand side; Trapezoidal rule

 
$$ (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 $$

 
Cannot do much unless we specialize \( f \) or use numerical integration.

Trapezoidal rule using the nodes:

 
$$ (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}) $$

 
\( \basphi_i(\xno{j})=\delta_{ij} \), so this formula collapses to one term:

 
$$ \begin{equation} (f,\basphi_i) \approx hf(\xno{i}),\quad i=1,\ldots,N-1\thinspace. \end{equation} $$

Same result as in collocation (interpolation) and the finite difference method!

Treating the right-hand side; Simpson's rule

 
$$ \int_\Omega g(x)dx \approx \frac{h}{6}\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), $$

 
Our case: \( g=f\basphi_i \). The sums collapse because \( \basphi_i=0 \) at most of the points.

 
$$ \begin{equation} (f,\basphi_i) \approx \frac{h}{3}(f_{i-\half} + f_i + f_{i+\half}) \end{equation} $$

Conclusions:

• While the finite difference method just samples \( f \) at \( x_i \), the finite element method applies an average (smoothing) of \( f \) around \( x_i \)
• On the left-hand side we have a term \( \sim hu'' \), and \( u'' \) also contribute to smoothing
• There is some inherent smoothing in the finite element method

Finite element approximation vs finite differences

With Trapezoidal integration of \( (f,\basphi_i) \), the finite element metod essentially solve

 
$$ \begin{equation} u + \frac{h^2}{6} u'' = f,\quad u'(0)=u'(L)=0, \end{equation} $$

 
by the finite difference method

 
$$ \begin{equation} [u + \frac{h^2}{6} D_x D_x u = f]_i \end{equation} $$

With Simpson integration of \( (f,\basphi_i) \) we essentially solve

 
$$ \begin{equation} [u + \frac{h^2}{6} D_x D_x u = \bar f]_i, \end{equation} $$

 
where

 
$$ \bar f_i = \frac{1}{3}(f_{i-1/2} + f_i + f_{i+1/2}) $$

Note: as \( h\rightarrow 0 \), \( hu''\rightarrow 0 \) and \( \bar f_i\rightarrow f_i \).

Making finite elements behave as finite differences

• Can we adjust the finite element method so that we do not get the extra \( hu'' \) smoothing term and averaging of \( f \)?
• This is sometimes important in time-dependent problems to incorporate good properties of finite differences into finite elements

Result:

• Compute all integrals by the Trapezoidal method and P1 elements
• Specifically, the coefficient matrix becomes diagonal ("lumped") - no linear system (!)
• Loss of accuracy? The Trapezoidal rule has error \( \Oof{h^2} \), the same as the approximation error in P1 elements

Limitations of the nodes and element concepts

So far,

• Nodes: points for defining \( \basphi_i \) and computing \( u \) values
• Elements: subdomain (containing a few nodes)
• This is a common notion of nodes and elements

One problem:

• Our algorithms need nodes at the element boundaries
• This is often not desirable, so we need to throw the nodes and elements arrays away and find a more generalized element concept

A generalized element concept

• We introduce cell for the subdomain that we up to now called element
• A cell has vertices (interval end points)
• Nodes are, almost as before, points where we want to compute unknown functions
• Degrees of freedom is what the \( c_j \) represent (usually function values at nodes)

The concept of a finite element

1. a reference cell in a local reference coordinate system
2. a set of basis functions \( \refphi_r \) defined on the cell
3. a set of degrees of freedom (e.g., function values) that uniquely determine the basis functions such that \( \refphi_r=1 \) for degree of freedom number \( r \) and \( \refphi_r=0 \) for all other degrees of freedom
4. a mapping between local and global degree of freedom numbers (dof map)
5. a geometric mapping of the reference cell onto to cell in the physical domain: \( [-1,1]\ \Rightarrow\ [x_L,x_R] \)

Implementation; basic data structures

• Cell vertex coordinates: vertices (equals nodes for P1 elements)
• Element vertices: cell[e][r] holds global vertex number of local vertex no r in element e (same as elements for P1 elements)
• dof_map[e,r] maps local dof r in element e to global dof number (same as elements for Pd elements)

The assembly process now applies dof_map:

A[dof_map[e][r], dof_map[e][s]] += A_e[r,s]
b[dof_map[e][r]] += b_e[r]

Implementation; example with P2 elements

vertices = [0, 0.4, 1]
cells = [[0, 1], [1, 2]]
dof_map = [[0, 1, 2], [2, 3, 4]]

Implementation; example with P0 elements

Example: Same mesh, but \( u \) is piecewise constant in each cell (P0 element). Same vertices and cells, but

dof_map = [[0], [1]]

May think of one node in the middle of each element.

Example on doing the algorithmic steps

# Use modified fe_approx1D module
from fe_approx1D_numint import *

x = sp.Symbol('x')
f = x*(1 - x)

N_e = 10
# Create mesh with P3 (cubic) elements
vertices, cells, dof_map = mesh_uniform(N_e, d=3, Omega=[0,1])

# Create basis functions on the mesh
phi = [basis(len(dof_map[e])-1) for e in range(N_e)]

# Create linear system and solve it
A, b = assemble(vertices, cells, dof_map, phi, f)
c = np.linalg.solve(A, b)

# Make very fine mesh and sample u(x) on this mesh for plotting
x_u, u = u_glob(c, vertices, cells, dof_map,
                resolution_per_element=51)
plot(x_u, u)

Approximating a parabola by P0 elements

The approximate function automates the steps in the previous slide:

from fe_approx1D_numint import *
x=sp.Symbol("x")
for N_e in 4, 8:
    approximate(x*(1-x), d=0, N_e=N_e, Omega=[0,1])

Computing the error of the approximation; principles

 
$$ L^2 \hbox{ error: }\quad ||e||_{L^2} = \left(\int_{\Omega} e^2 dx\right)^{1/2}$$

Accurate approximation of the integral:

• Sample \( u(x) \) at many points in each element (call u_glob, returns x and u)
• Use the Trapezoidal rule based on the samples
• It is important to integrate \( u \) accurately over the elements
• (In a finite difference method we would just sample the mesh point values)

Computing the error of the approximation; details

# Given c, compute x and u values on a very fine mesh
x, u = u_glob(c, vertices, cells, dof_map,
              resolution_per_element=101)
# Compute the error on the very fine mesh
e = f(x) - u
e2 = e**2
# Vectorized Trapezoidal rule
E = np.sqrt(0.5*np.sum((e2[:-1] + e2[1:])*(x[1:] - x[:-1]))

How does the error depend on \( h \) and \( d \)?

Theory and experiments show that the least squares or projection/Galerkin method in combination with Pd elements of equal length \( h \) has an error

 
$$ \begin{equation} ||e||_{L^2} = Ch^{d+1} \tag{28} \end{equation} $$

 
where \( C \) depends on \( f \), but not on \( h \) or \( d \).

Cubic Hermite polynomials; definition

• Can we construct \( \basphi_i(x) \) with continuous derivatives? Yes!

Consider a reference cell \( [-1,1] \). We introduce two nodes, \( X=-1 \) and \( X=1 \). The degrees of freedom are

• 0: value of function at \( X=-1 \)
• 1: value of first derivative at \( X=-1 \)
• 2: value of function at \( X=1 \)
• 3: value of first derivative at \( X=1 \)

Cubic Hermite polynomials; derivation

4 constraints on \( \refphi_r \) (1 for dof \( r \), 0 for all others):

• \( \refphi_0(\Xno{0}) = 1 \), \( \refphi_0(\Xno{1}) = 0 \), \( \refphi_0'(\Xno{0}) = 0 \), \( \refphi_0'(\Xno{1}) = 0 \)
• \( \refphi_1'(\Xno{0}) = 1 \), \( \refphi_1'(\Xno{1}) = 0 \), \( \refphi_1(\Xno{0}) = 0 \), \( \refphi_1(\Xno{1}) = 0 \)
• \( \refphi_2(\Xno{1}) = 1 \), \( \refphi_2(\Xno{0}) = 0 \), \( \refphi_2'(\Xno{0}) = 0 \), \( \refphi_2'(\Xno{1}) = 0 \)
• \( \refphi_3'(\Xno{1}) = 1 \), \( \refphi_3'(\Xno{0}) = 0 \), \( \refphi_3(\Xno{0}) = 0 \), \( \refphi_3(\Xno{1}) = 0 \)

This gives 4 linear, coupled equations for each \( \refphi_r \) to determine the 4 coefficients in the cubic polynomial

Cubic Hermite polynomials; result

 
$$ \begin{align} \refphi_0(X) &= 1 - \frac{3}{4}(X+1)^2 + \frac{1}{4}(X+1)^3\\ \refphi_1(X) &= -(X+1)(1 - \half(X+1))^2\\ \refphi_2(X) &= \frac{3}{4}(X+1)^2 - \half(X+1)^3\\ \refphi_3(X) &= -\half(X+1)(\half(X+1)^2 - (X+1))\\ \end{align} $$

Numerical integration

• \( \int_\Omega f\basphi_idx \) must in general be computed by numerical integration
• Numerical integration is often used for the matrix too

Common form of a numerical integration rule:

 
$$ \begin{equation} \int_{-1}^{1} g(X)dX \approx \sum_{j=0}^M w_jg(\bar X_j), \end{equation} $$

 
where

• \( \bar X_j \) are integration points
• \( w_j \) are integration weights

Different rules correspond to different choices of points and weights

The Midpoint rule

Simplest possibility: the Midpoint rule,

 
$$ \begin{equation} \int_{-1}^{1} g(X)dX \approx 2g(0),\quad \bar X_0=0,\ w_0=2, \end{equation} $$

Exact for linear integrands

Newton-Cotes rules

• Idea: use a fixed, uniformly distributed set of points in \( [-1,1] \)
• The points often coincides with nodes
• Very useful for making \( \basphi_i\basphi_j=0 \) and get diagonal ("mass") matrices ("lumping")

The Trapezoidal rule:

 
$$ \begin{equation} \int_{-1}^{1} g(X)dX \approx g(-1) + g(1),\quad \bar X_0=-1,\ \bar X_1=1,\ w_0=w_1=1, \tag{29} \end{equation} $$

Simpson's rule:

 
$$ \begin{equation} \int_{-1}^{1} g(X)dX \approx \frac{1}{3}\left(g(-1) + 4g(0) + g(1)\right), \end{equation} $$

 
where

 
$$ \begin{equation} \bar X_0=-1,\ \bar X_1=0,\ \bar X_2=1,\ w_0=w_2=\frac{1}{3},\ w_1=\frac{4}{3} \end{equation} $$

Gauss-Legendre rules with optimized points

• Optimize the location of points to get higher accuracy
• Gauss-Legendre rules (quadrature) adjust points and weights to integrate polynomials exactly

 
$$ \begin{align} M=1&:\quad \bar X_0=-\frac{1}{\sqrt{3}},\ \bar X_1=\frac{1}{\sqrt{3}},\ w_0=w_1=1\\ M=2&:\quad \bar X_0=-\sqrt{\frac{3}{{5}}},\ \bar X_0=0,\ \bar X_2= \sqrt{\frac{3}{{5}}},\ w_0=w_2=\frac{5}{9},\ w_1=\frac{8}{9} \end{align} $$

• \( M=1 \): integrates 3rd degree polynomials exactly
• \( M=2 \): integrates 5th degree polynomials exactly
• In general, \( M \)-point rule integrates a polynomial of degree \( 2M+1 \) exactly.

See numint.py for a large collection of Gauss-Legendre rules.

Approximation of functions in 2D

Inner product in 2D:

 
$$ \begin{equation} (f,g) = \int_\Omega f(x,y)g(x,y) dx dy \end{equation} $$

Least squares and project/Galerkin lead to a linear system

 
$$ \begin{align*} \sum_{j\in\If} A_{i,j}c_j &= b_i,\quad i\in\If\\ A_{i,j} &= (\baspsi_i,\baspsi_j)\\ b_i &= (f,\baspsi_i) \end{align*} $$

Challenge: How to construct 2D basis functions \( \baspsi_i(x,y) \)?

2D basis functions as tensor products of 1D functions

Use a 1D basis for \( x \) variation and a similar for \( y \) variation:

 
$$ \begin{align} V_x &= \mbox{span}\{ \hat\baspsi_0(x),\ldots,\hat\baspsi_{N_x}(x)\} \tag{30}\\ V_y &= \mbox{span}\{ \hat\baspsi_0(y),\ldots,\hat\baspsi_{N_y}(y)\} \tag{31} \end{align} $$

The 2D vector space can be defined as a tensor product \( V = V_x\otimes V_y \) with basis functions

 
$$ \baspsi_{p,q}(x,y) = \hat\baspsi_p(x)\hat\baspsi_q(y) \quad p\in\Ix,q\in\Iy\tp $$

Tensor products

Given two vectors \( a=(a_0,\ldots,a_M) \) and \( b=(b_0,\ldots,b_N) \) their outer tensor product, also called the dyadic product, is \( p=a\otimes b \), defined through

 
$$ p_{i,j}=a_ib_j,\quad i=0,\ldots,M,\ j=0,\ldots,N\tp$$

 
Note: \( p \) has two indices (as a matrix or two-dimensional array)

Example: 2D basis as tensor product of 1D spaces,

 
$$ \baspsi_{p,q}(x,y) = \hat\baspsi_p(x)\hat\baspsi_q(y), \quad p\in\Ix,q\in\Iy$$

Double or single index?

The 2D basis can employ a double index and double sum:

 
$$ u = \sum_{p\in\Ix}\sum_{q\in\Iy} c_{p,q}\baspsi_{p,q}(x,y) $$

Or just a single index:

 
$$ u = \sum_{j\in\If} c_j\baspsi_j(x,y)$$

 
with

 
$$ \baspsi_i(x,y) = \hat\baspsi_p(x)\hat\baspsi_q(y), \quad i=p N_y + q\hbox{ or } i=q N_x + p $$

Example on 2D (bilinear) basis functions; formulas

In 1D we use the basis

 
$$ \{ 1, x \} $$

2D tensor product (all combinations):

 
$$ \baspsi_{0,0}=1,\quad \baspsi_{1,0}=x, \quad \baspsi_{0,1}=y, \quad \baspsi_{1,1}=xy $$

 
or with a single index:

 
$$ \baspsi_0=1,\quad \baspsi_1=x, \quad \baspsi_2=y,\quad\baspsi_3 =xy $$

See notes for details of a hand-calculation.

Example on 2D (bilinear) basis functions; plot

Quadratic \( f(x,y) = (1+x^2)(1+2y^2) \) (left), bilinear \( u \) (right):

Implementation; principal changes to the 1D code

Very small modification of approx1D.py:

• Omega = [[0, L_x], [0, L_y]]
• Symbolic integration in 2D
• Construction of 2D (tensor product) basis functions

Implementation; 2D integration

import sympy as sp

integrand = psi[i]*psi[j]
I = sp.integrate(integrand,
                 (x, Omega[0][0], Omega[0][1]),
                 (y, Omega[1][0], Omega[1][1]))

# Fall back on numerical integration if symbolic integration
# was unsuccessful
if isinstance(I, sp.Integral):
    integrand = sp.lambdify([x,y], integrand)
    I = sp.mpmath.quad(integrand,
                       [Omega[0][0], Omega[0][1]],
                       [Omega[1][0], Omega[1][1]])

Implementation; 2D basis functions

Tensor product of 1D "Taylor-style" polynomials \( x^i \):

def taylor(x, y, Nx, Ny):
    return [x**i*y**j for i in range(Nx+1) for j in range(Ny+1)]

Tensor product of 1D sine functions \( \sin((i+1)\pi x) \):

def sines(x, y, Nx, Ny):
    return [sp.sin(sp.pi*(i+1)*x)*sp.sin(sp.pi*(j+1)*y)
            for i in range(Nx+1) for j in range(Ny+1)]

Complete code in approx2D.py

Implementation; application

\( f(x,y) = (1+x^2)(1+2y^2) \)

>>> from approx2D import *
>>> f = (1+x**2)*(1+2*y**2)
>>> psi = taylor(x, y, 1, 1)
>>> Omega = [[0, 2], [0, 2]]
>>> u, c = least_squares(f, psi, Omega)
>>> print u
8*x*y - 2*x/3 + 4*y/3 - 1/9
>>> print sp.expand(f)
2*x**2*y**2 + x**2 + 2*y**2 + 1

Implementation; trying a perfect expansion

Add higher powers to the basis such that \( f\in V \):

>>> psi = taylor(x, y, 2, 2)
>>> u, c = least_squares(f, psi, Omega)
>>> print u
2*x**2*y**2 + x**2 + 2*y**2 + 1
>>> print u-f
0

Expected: \( u=f \) when \( f\in V \)

Generalization to 3D

Key idea:

 
$$ V = V_x\otimes V_y\otimes V_z$$

 
$$ \begin{align*} \baspsi_{p,q,r}(x,y,z) &= \hat\baspsi_p(x)\hat\baspsi_q(y)\hat\baspsi_r(z)\\ u(x,y,z) &= \sum_{p\in\Ix}\sum_{q\in\Iy}\sum_{r\in\Iz} c_{p,q,r} \baspsi_{p,q,r}(x,y,z) \end{align*} $$

Finite elements in 2D and 3D

The two great advantages of the finite element method:

• Can handle complex-shaped domains in 2D and 3D
• Can easily provide higher-order polynomials in the approximation

Finite elements in 1D: mostly for learning, insight, debugging

Examples on cell types

2D:

• triangles
• quadrilaterals

3D:

• tetrahedra
• hexahedra

Rectangular domain with 2D P1 elements

Deformed geometry with 2D P1 elements

Rectangular domain with 2D Q1 elements

Basis functions over triangles in the physical domain

The P1 triangular 2D element: \( u \) is linear \( ax + by + c \) over each triangular cell

Basic features of 2D P1 elements

• \( \basphi_r(X,Y) \) is a linear function over each element
• Cells = triangles
• Vertices = corners of the cells
• Nodes = vertices
• Degrees of freedom = function values at the nodes

Linear mapping of reference element onto general triangular cell

\( \basphi_i \): pyramid shape, composed of planes

• \( \basphi_i(X,Y) \) varies linearly over an element
• \( \basphi_i=1 \) at vertex (node) \( i \), 0 at all other vertices (nodes)

Element matrices and vectors

• As in 1D, the contribution from one cell to the matrix involves just a few numbers, collected in the element matrix and vector
• \( \basphi_i\basphi_j\neq 0 \) only if \( i \) and \( j \) are degrees of freedom (vertices/nodes) in the same element
• The 2D P1 has a \( 3\times 3 \) element matrix

Basis functions over triangles in the reference cell

 
$$ \begin{align} \refphi_0(X,Y) &= 1 - X - Y\\ \refphi_1(X,Y) &= X\\ \refphi_2(X,Y) &= Y \end{align} $$

Higher-degree \( \refphi_r \) introduce more nodes (dof = node values)

2D P1, P2, P3, P4, P5, and P6 elements

P1 elements in 1D, 2D, and 3D

P2 elements in 1D, 2D, and 3D

• Interval, triangle, tetrahedron: simplex element (plural quick-form: simplices)
• Side of the cell is called face
• Thetrahedron has also edges

Affine mapping of the reference cell; formula

Mapping of local \( \X = (X,Y) \) coordinates in the reference cell to global, physical \( \x = (x,y) \) coordinates:

 
$$ \begin{equation} \x = \sum_{r} \refphi_r^{(1)}(\X)\xdno{q(e,r)} \tag{32} \end{equation} $$

where

• \( r \) runs over the local vertex numbers in the cell
• \( \xdno{i} \) are the \( (x,y) \) coordinates of vertex \( i \)
• \( \refphi_r^{(1)} \) are P1 basis functions

This mapping preserves the straight/planar faces and edges.

Affine mapping of the reference cell; figure

Isoparametric mapping of the reference cell

Idea: Use the basis functions of the element (not only the P1 functions) to map the element

 
$$ \begin{equation} \x = \sum_{r} \refphi_r(\X)\xdno{q(e,r)} \tag{33} \end{equation} $$

Advantage: higher-order polynomial basis functions now map the reference cell to a curved triangle or tetrahedron.

Computing integrals

Integrals must be transformed from \( \Omega^{(e)} \) (physical cell) to \( \tilde\Omega^r \) (reference cell):

 
$$ \begin{align} \int_{\Omega^{(e)}}\basphi_i (\x) \basphi_j (\x) \dx &= \int_{\tilde\Omega^r} \refphi_i (\X) \refphi_j (\X) \det J\, \dX\\ \int_{\Omega^{(e)}}\basphi_i (\x) f(\x) \dx &= \int_{\tilde\Omega^r} \refphi_i (\X) f(\x(\X)) \det J\, \dX \end{align} $$

 
where \( \dx = dx dy \) or \( \dx = dxdydz \) and \( \det J \) is the determinant of the Jacobian of the mapping \( \x(\X) \).

 
$$ \begin{equation} J = \left[\begin{array}{cc} \frac{\partial x}{\partial X} & \frac{\partial x}{\partial Y}\\ \frac{\partial y}{\partial X} & \frac{\partial y}{\partial Y} \end{array}\right], \quad \det J = \frac{\partial x}{\partial X}\frac{\partial y}{\partial Y} - \frac{\partial x}{\partial Y}\frac{\partial y}{\partial X} \tag{34} \end{equation} $$

Affine mapping (32): \( \det J=2\Delta \), \( \Delta = \hbox{cell volume} \)

!slide

Remark on going from 1D to 2D/3D

Differential equation models

Our aim is to extend the ideas for approximating \( f \) by \( u \), or solving

 
$$ u = f $$

to real differential equations like[[[

 
$$ -u'' + bu = f,\quad u(0)=1,\ u'(L)=D $$

Three methods are addressed:

1. least squares
2. Galerkin/projection
3. collocation (interpolation)

Method 2 will be totally dominating!

Abstract differential equation

 
$$ \begin{equation} \mathcal{L}(u) = 0,\quad x\in\Omega \end{equation} $$

Examples (1D problems):

 
$$ \begin{align} \mathcal{L}(u) &= \frac{d^2u}{dx^2} - f(x), \tag{35}\\ \mathcal{L}(u) &= \frac{d}{dx}\left(\dfc(x)\frac{du}{dx}\right) + f(x), \tag{36}\\ \mathcal{L}(u) &= \frac{d}{dx}\left(\dfc(u)\frac{du}{dx}\right) - au + f(x), \tag{37}\\ \mathcal{L}(u) &= \frac{d}{dx}\left(\dfc(u)\frac{du}{dx}\right) + f(u,x) \tag{38} \end{align} $$

Abstract boundary conditions

 
$$ \begin{equation} {\cal B}_0(u)=0,\ x=0,\quad {\cal B}_1(u)=0,\ x=L \end{equation} $$

Examples:

 
$$ \begin{align} \mathcal{B}_i(u) &= u - g,\quad &\hbox{Dirichlet condition}\\ \mathcal{B}_i(u) &= -\dfc \frac{du}{dx} - g,\quad &\hbox{Neumann condition}\\ \mathcal{B}_i(u) &= -\dfc \frac{du}{dx} - h(u-g),\quad &\hbox{Robin condition} \end{align} $$

Reminder about notation

• \( \uex(x) \) is the symbol for the exact solution of \( \mathcal{L}(\uex)=0 \)
• \( u(x) \) denotes an approximate solution
• We seek \( u\in V \)
• \( V = \hbox{span}\{ \baspsi_0(x),\ldots,\baspsi_N(x)\} \), \( V \) has basis \( \sequencei{\baspsi} \)
• \( \If =\{0,\ldots,N\} \) is an index set
• \( u(x) = \sum_{j\in\If} c_j\baspsi_j(x) \)
• Inner product: \( (u,v) = \int_\Omega uv\dx \)
• Norm: \( ||u||=\sqrt{(u,u)} \)

New topics

Much is similar to approximating a function (solving \( u=f \)), but two new topics are needed:

• Variational formulation of the differential equation problem (including integration by parts)
• Handling of boundary conditions

Residual-minimizing principles

• When solving \( u=f \) we knew the error \( e=f-u \) and could use principles for minimizing the error
• When solving \( \mathcal{L}(\uex)=0 \) we do not know \( \uex \) and cannot work with the error \( e=\uex - u \)
• We only have the error in the equation: the residual \( R \)

Inserting \( u=\sum_jc_j\baspsi_j \) in \( \mathcal{L}=0 \) gives a residual

 
$$ \begin{equation} R = \mathcal{L}(u) = \mathcal{L}(\sum_j c_j \baspsi_j) \neq 0 \end{equation} $$

Goal: minimize \( R \) wrt \( \sequencei{c} \) (and hope it makes a small \( e \) too)

 
$$ R=R(c_0,\ldots,c_N; x)$$

The least squares method

Idea: minimize

 
$$ \begin{equation} E = ||R||^2 = (R,R) = \int_{\Omega} R^2 dx \end{equation} $$

Minimization wrt \( \sequencei{c} \) implies

 
$$ \begin{equation} \frac{\partial E}{\partial c_i} = \int_{\Omega} 2R\frac{\partial R}{\partial c_i} dx = 0\quad \Leftrightarrow\quad (R,\frac{\partial R}{\partial c_i})=0,\quad i\in\If \tag{39} \end{equation} $$

\( N+1 \) equations for \( N+1 \) unknowns \( \sequencei{c} \)

The Galerkin method

Idea: make \( R \) orthogonal to \( V \),

 
$$ \begin{equation} (R,v)=0,\quad \forall v\in V \tag{40} \end{equation} $$

This implies

 
$$ \begin{equation} (R,\baspsi_i)=0,\quad i\in\If \tag{41} \end{equation} $$

\( N+1 \) equations for \( N+1 \) unknowns \( \sequencei{c} \)

The Method of Weighted Residuals

Generalization of the Galerkin method: demand \( R \) orthogonal to some space \( W \), possibly \( W\neq V \):

 
$$ \begin{equation} (R,v)=0,\quad \forall v\in W \tag{42} \end{equation} $$

If \( \{w_0,\ldots,w_N\} \) is a basis for \( W \):

 
$$ \begin{equation} (R,w_i)=0,\quad i\in\If \tag{43} \end{equation} $$

• \( N+1 \) equations for \( N+1 \) unknowns \( \sequencei{c} \)
• Weighted residual with \( w_i = \partial R/\partial c_i \) gives least squares

Terminology: test and trial Functions

• \( \baspsi_j \) used in \( \sum_jc_j\baspsi_j \) is called trial function
• \( \baspsi_i \) or \( w_i \) used as weight in Galerkin's method is called test function

The collocation method

Idea: demand \( R=0 \) at \( N+1 \) points

 
$$ \begin{equation} R(\xno{i}; c_0,\ldots,c_N)=0,\quad i\in\If \tag{44} \end{equation} $$

Note: The collocation method is a weighted residual method with delta functions as weights

 
$$ 0 = \int_\Omega R(x;c_0,\ldots,c_N) \delta(x-\xno{i})\dx = R(\xno{i}; c_0,\ldots,c_N)$$

 
$$ \begin{equation} \hbox{property of } \delta(x):\quad \int_{\Omega} f(x)\delta (x-\xno{i}) dx = f(\xno{i}),\quad \xno{i}\in\Omega \tag{45} \end{equation} $$

Examples on using the principles

The first model problem

 
$$ \begin{equation} -u''(x) = f(x),\quad x\in\Omega=[0,L],\quad u(0)=0,\ u(L)=0 \tag{46} \end{equation} $$

Basis functions:

 
$$ \begin{equation} \baspsi_i(x) = \sinL{i},\quad i\in\If \tag{47} \end{equation} $$

The residual:

 
$$ \begin{align} R(x;c_0,\ldots,c_N) &= u''(x) + f(x),\nonumber\\ &= \frac{d^2}{dx^2}\left(\sum_{j\in\If} c_j\baspsi_j(x)\right) + f(x),\nonumber\\ &= -\sum_{j\in\If} c_j\baspsi_j''(x) + f(x) \tag{48} \end{align} $$

Boundary conditions

Since \( u(0)=u(L)=0 \) we must ensure that all \( \baspsi_i(0)=\baspsi_i(L)=0 \). Then

 
$$ u(0) = \sum_jc_j\baspsi_j(0) = 0,\quad u(L) = \sum_jc_j\baspsi_j(L) $$

• \( u \) known: Dirichlet boundary condition
• \( u' \) known: Neumann boundary condition
• Must have \( \baspsi_i=0 \) where Dirichlet conditions apply

The least squares method; principle

 
$$ (R,\frac{\partial R}{\partial c_i}) = 0,\quad i\in\If $$

 
$$ \begin{equation} \frac{\partial R}{\partial c_i} = \frac{\partial}{\partial c_i} \left(\sum_{j\in\If} c_j\baspsi_j''(x) + f(x)\right) = \baspsi_i''(x) \end{equation} $$

Because:

 
$$ \frac{\partial}{\partial c_i}\left(c_0\baspsi_0'' + c_1\baspsi_1'' + \cdots + c_{i-1}\baspsi_{i-1}'' + \color{blue}{c_i\baspsi_{i}''} + c_{i+1}\baspsi_{i+1}'' + \cdots + c_N\baspsi_N'' \right) = \baspsi_{i}'' $$

The least squares method; equation system

 
$$ \begin{equation} (\sum_j c_j \baspsi_j'' + f,\baspsi_i'')=0,\quad i\in\If \end{equation} $$

Rearrangement:

 
$$ \begin{equation} \sum_{j\in\If}(\baspsi_i'',\baspsi_j'')c_j = -(f,\baspsi_i''),\quad i\in\If \end{equation} $$

This is a linear system

 
$$ \begin{equation*} \sum_{j\in\If}A_{i,j}c_j = b_i,\quad i\in\If \end{equation*} $$

 
with

 
$$ \begin{align} A_{i,j} &= (\baspsi_i'',\baspsi_j'')\nonumber\\ & = \pi^4(i+1)^2(j+1)^2L^{-4}\int_0^L \sinL{i}\sinL{j}\, dx\nonumber\\ &= \left\lbrace \begin{array}{ll} {1\over2}L^{-3}\pi^4(i+1)^4 & i=j \\ 0, & i\neq j \end{array}\right. \\ b_i &= -(f,\baspsi_i'') = (i+1)^2\pi^2L^{-2}\int_0^Lf(x)\sinL{i}\, dx \end{align} $$

Orthogonality of the basis functions gives diagonal matrix

Useful property:

 
$$ \begin{equation} \int\limits_0^L \sinL{i}\sinL{j}\, dx = \delta_{ij},\quad \quad\delta_{ij} = \left\lbrace \begin{array}{ll} \half L & i=j \\ 0, & i\neq j \end{array}\right. \end{equation} $$

\( \Rightarrow\ (\baspsi_i'',\baspsi_j'') = \delta_{ij} \), i.e., diagonal \( A_{i,j} \), and we can easily solve for \( c_i \):

 
$$ \begin{equation} c_i = \frac{2L}{\pi^2(i+1)^2}\int_0^Lf(x)\sinL{i}\, dx \tag{49} \end{equation} $$

Least squares method; solution

Let's sympy do the work (\( f(x)=2 \)):

from sympy import *
import sys

i, j = symbols('i j', integer=True)
x, L = symbols('x L')
f = 2
a = 2*L/(pi**2*(i+1)**2)
c_i = a*integrate(f*sin((i+1)*pi*x/L), (x, 0, L))
c_i = simplify(c_i)
print c_i

 
$$ \begin{equation} c_i = 4 \frac{L^{2} \left(\left(-1\right)^{i} + 1\right)}{\pi^{3} \left(i^{3} + 3 i^{2} + 3 i + 1\right)},\quad u(x) = \sum_{k=0}^{N/2} \frac{8L^2}{\pi^3(2k+1)^3}\sinL{2k}\tp \end{equation} $$

Fast decay: \( c_2 = c_0/27 \), \( c_4=c_0/125 \) - only one term might be good enough:

 
$$ \begin{equation*} u(x) \approx \frac{8L^2}{\pi^3}\sin\left(\pi\frac{x}{L}\right)\tp \end{equation*} $$

The Galerkin method; principle

\( R=u''+f \):

 
$$ \begin{equation*} (u''+f,v)=0,\quad \forall v\in V, \end{equation*} $$

 
or

 
$$ \begin{equation} (u'',v) = -(f,v),\quad\forall v\in V \end{equation} $$

This is a variational formulation of the differential equation problem.

\( \forall v\in V \) means for all basis functions:

 
$$ \begin{equation} (\sum_{j\in\If} c_j\baspsi_j'', \baspsi_i)=-(f,\baspsi_i),\quad i\in\If \end{equation} $$

The Galerkin method; solution

Since \( \baspsi_i''\propto \baspsi_i \), Galerkin's method gives the same linear system and the same solution as the least squares method (in this particular example).

The collocation method

\( R=0 \) (i.e.,the differential equation) must be satisfied at \( N+1 \) points:

 
$$ \begin{equation} -\sum_{j\in\If} c_j\baspsi_j''(\xno{i}) = f(\xno{i}),\quad i\in\If \end{equation} $$

This is a linear system \( \sum_j A_{i,j}=b_i \) with entries

 
$$ \begin{equation*} A_{i,j}=-\baspsi_j''(\xno{i})= (j+1)^2\pi^2L^{-2}\sin\left((j+1)\pi \frac{x_i}{L}\right), \quad b_i=2 \end{equation*} $$

Choose: \( N=0 \), \( x_0=L/2 \)

 
$$ c_0=2L^2/\pi^2 $$

Comparison of the methods

• Exact solution: \( u(x)=x(L-x) \)
• Galerkin or least squares (\( N=0 \)): \( u(x)=8L^2\pi^{-3}\sin (\pi x/L) \)
• Collocation method (\( N=0 \)): \( u(x)=2L^2\pi^{-2}\sin (\pi x/L) \).
• Max error in Galerkin/least sq.: \( -0.008L^2 \)
• Max error in collocation: \( 0.047L^2 \)

Useful techniques

Integration by parts

Second-order derivatives will hereafter be integrated by parts

 
$$ \begin{align} \int_0^L u''(x)v(x) dx &= - \int_0^Lu'(x)v'(x)dx + [vu']_0^L\nonumber\\ &= - \int_0^Lu'(x)v'(x) dx + u'(L)v(L) - u'(0)v(0) \tag{50} \end{align} $$

Motivation:

• Lowers the order of derivatives
• Gives more symmetric forms (incl. matrices)
• Enables easy handling of Neumann boundary conditions
• Finite element basis functions \( \basphi_i \) have discontinuous derivatives (at cell boundaries) and are not suited for terms with \( \basphi_i'' \)

Boundary function; principles

• What about nonzero Dirichlet conditions? Say \( u(L)=D \)
• We always require \( \baspsi_i(L)=0 \) (i.e., \( \baspsi_i=0 \) where Dirichlet conditions applies)
• Problem: \( u(L) = \sum_j c_j\baspsi_j(L)=\sum_j c_j\cdot 0=0\neq D \) - always
• Solution: \( u(x) = B(x) + \sum_j c_j\baspsi_j(x) \)
• \( B(x) \): user-constructed boundary function that fulfills the Dirichlet conditions
• If \( u(L)=D \), \( B(L)=D \)
• No restrictions of how \( B(x) \) varies in the interior of \( \Omega \)

Boundary function; example (1)

Dirichlet conditions: \( u(0)=C \) and \( u(L)=D \). Choose for example

 
$$ B(x) = \frac{1}{L}(C(L-x) + Dx):\qquad B(0)=C,\ B(L)=D $$

 
$$ \begin{equation} u(x) = B(x) + \sum_{j\in\If} c_j\baspsi_j(x), \tag{51} \end{equation} $$

 
$$ u(0) = B(0)= C,\quad u(L) = B(L) = D $$

Boundary function; example (2)

Dirichlet condition: \( u(L)=D \). Choose for example

 
$$ B(x) = D:\qquad B(L)=D $$

 
$$ \begin{equation} u(x) = B(x) + \sum_{j\in\If} c_j\baspsi_j(x), \tag{51} \end{equation} $$

 
$$ u(L) = B(L) = D $$

Impact of the boundary function on the space where we seek the solution

• \( \sequencei{\baspsi} \) is a basis for \( V \)
• \( \sum_{j\in\If}c_j\baspsi_j(x)\in V \)
• But \( u\not\in V \)!
• Reason: say \( u(0)=C \) and \( u\in V \) (any \( v\in V \) has \( v(0)=C \), then \( 2u\not\in V \) because \( 2u(0)=2C \)
• When \( u(x) = B(x) + \sum_{j\in\If}c_j\baspsi_j(x) \), \( B \neq 0 \), \( B\not\in V \) (in general) and \( u\not\in V \), but \( (u-B)\in V \) since \( \sum_{j}c_j\baspsi_j\in V \)

Abstract notation for variational formulations

The finite element literature (and much FEniCS documentation) applies an abstract notation for the variational formulation:

Example on abstract notation

 
$$ -u''=f, \quad u'(0)=C,\ u(L)=D,\quad u=D + \sum_jc_j\baspsi_j$$

Variational formulation:

 
$$ \int_{\Omega} u' v'dx = \int_{\Omega} fvdx\quad - v(0)C \hbox{or}\quad (u',v') = (f,v) - v(0)C \quad\forall v\in V $$

Abstract formulation: finn \( (u-B)\in V \) such that

 
$$ a(u,v) = L(v)\quad \forall v\in V$$

We identify

 
$$ a(u,v) = (u',v'),\quad L(v) = (f,v) -v(0)C $$

Bilinear and linear forms

• \( a(u,v) \) is a bilinear form
• \( L(v) \) is a linear form

Linear form means

 
$$ L(\alpha_1 v_1 + \alpha_2 v_2) =\alpha_1 L(v_1) + \alpha_2 L(v_2), $$

Bilinear form means

 
$$ \begin{align*} a(\alpha_1 u_1 + \alpha_2 u_2, v) &= \alpha_1 a(u_1,v) + \alpha_2 a(u_2, v), \\ a(u, \alpha_1 v_1 + \alpha_2 v_2) &= \alpha_1 a(u,v_1) + \alpha_2 a(u, v_2) \end{align*} $$

In nonlinear problems: Find \( (u-B)\in V \) such that \( F(u;v)=0\ \forall v\in V \)

The linear system associated with abstract form

 
$$ a(u,v) = L(v)\quad \forall v\in V\quad\Leftrightarrow\quad a(u,\baspsi_i) = L(\baspsi_i)\quad i\in\If$$

We can now derive the corresponding linear system once and for all:

 
$$ a(\sum_{j\in\If} c_j \baspsi_j,\baspsi_i)c_j = L(\baspsi_i)\quad i\in\If$$

Because of linearity,

 
$$ \sum_{j\in\If} \underbrace{a(\baspsi_j,\baspsi_i)}_{A_{i,j}}c_j = \underbrace{L(\baspsi_i)}_{b_i}\quad i\in\If$$

Equivalence with minimization problem

If \( a(u,v)=a(v,u) \),

 
$$ a(u,v)=L(v)\quad\forall v\in V,$$

is equivalent to minimizing the functional

 
$$ F(v) = {\half}a(v,v) - L(v) $$

 
over all functions \( v\in V \). That is,

 
$$ F(u)\leq F(v)\quad \forall v\in V\tp $$

• Much used in the early days of finite elements
• Still much used in structural analysis and elasticity
• Not as general as Galerkin's method (since \( a(u,v)=a(v,u) \))

Examples on variational formulations

Variable coefficient; problem

 
$$ \begin{equation} -\frac{d}{dx}\left( \dfc(x)\frac{du}{dx}\right) = f(x),\quad x\in\Omega =[0,L],\ u(0)=C,\ u(L)=D \end{equation} $$

• Variable coefficient \( \dfc(x) \)
• Nonzero Dirichlet conditions at \( x=0 \) and \( x=L \)
• Must have \( \baspsi_i(0)=\baspsi_i(L)=0 \)
• \( V = \hbox{span}\{\baspsi_0,\ldots,\baspsi_N\} \)
• \( v\in V \): \( v(0)=v(L)=0 \)

 
$$ u(x) = B(x) + \sum_{j\in\If} c_j\baspsi_i(x) $$

 
$$ B(x) = C + \frac{1}{L}(D-C)x$$

Variable coefficient; variational formulation (1)

 
$$ R = -\frac{d}{dx}\left( a\frac{du}{dx}\right) -f $$

Galerkin's method:

 
$$ (R, v) = 0,\quad \forall v\in V, $$

or with integrals:

 
$$ \int_{\Omega} \left(\frac{d}{dx}\left( \dfc\frac{du}{dx}\right) -f\right)v \dx = 0,\quad \forall v\in V \tp $$

Variable coefficient; variational formulation (2)

Integration by parts:

 
$$ -\int_{\Omega} \frac{d}{dx}\left( \dfc(x)\frac{du}{dx}\right) v \dx = \int_{\Omega} \dfc(x)\frac{du}{dx}\frac{dv}{dx}\dx - \left[\dfc\frac{du}{dx}v\right]_0^L \tp $$

Boundary terms vanish since \( v(0)=v(L)=0 \)

Variable coefficient; linear system (the easy way)

With

 
$$ a(u,v) = (\dfc u', v),\quad L(v) = (f,v) $$

we can just use the formula for the linear system:

 
$$ \begin{align*} A_{i,j} &= a(\baspsi_j,\baspsi_i) = (\dfc \baspsi_j', \baspsi_i') = \int_\Omega \dfc \baspsi_j' \baspsi_i'\dx = \int_\Omega \baspsi_i' \dfc \baspsi_j'\dx = a(\baspsi_i,\baspsi_j) = A_{j,i}\\ b_i &= (f,\baspsi_i) = \int_\Omega f\baspsi_i\dx \end{align*} $$

Variable coefficient; linear system (full derivation)

\( v=\baspsi_i \) and \( u=B + \sum_jc_j\baspsi_j \):

 
$$ (\dfc B' + \dfc \sum_{j\in\If} c_j \baspsi_j', \baspsi_i') = (f,\baspsi_i), \quad i\in\If \tp $$

Reorder to form linear system:

 
$$ \sum_{j\in\If} (\dfc\baspsi_j', \baspsi_i')c_j = (f,\baspsi_i) + (a(D-C)L^{-1}, \baspsi_i'), \quad i\in\If \tp $$

This is \( \sum_j A_{i,j}c_j=b_i \) with

 
$$ \begin{align*} A_{i,j} &= (a\baspsi_j', \baspsi_i') = \int_{\Omega} \dfc(x)\baspsi_j'(x) \baspsi_i'(x)\dx\\ b_i &= (f,\baspsi_i) + (a(D-C)L^{-1},\baspsi_i')= \int_{\Omega} \left(f(x)\baspsi_i(x) + \dfc(x)\frac{D-C}{L}\baspsi_i'(x)\right) \dx \end{align*} $$

First-order derivative in the equation and boundary condition; problem

 
$$ \begin{equation} -u''(x) + bu'(x) = f(x),\quad x\in\Omega =[0,L],\ u(0)=C,\ u'(L)=E \end{equation} $$

New features:

• first-order derivative \( u' \) in the equation
• boundary condition with \( u' \): \( u'(L)=E \)

Initial steps:

• Must force \( \baspsi_i(0)=0 \) because of Dirichlet condition at \( x=0 \)
• Boundary function: \( B(x)=C(L-x) \) or just \( B(x)=C \)
• No requirements on \( \baspsi_i(L) \) (no Dirichlet condition at \( x=L \))

First-order derivative in the equation and boundary condition; details

 
$$ u = C + \sum_{j\in\If} c_j \baspsi_i(x)$$

Galerkin's method: multiply by \( v \), integrate over \( \Omega \), integrate by parts.

 
$$ (-u'' + bu' - f, v) = 0,\quad\forall v\in V$$

 
$$ (u',v') + (bu',v) = (f,v) + [u' v]_0^L, \quad\forall v\in V$$

Now, \( [u' v]_0^L = u'(L)v(L) = E v(L) \) because \( v(0)=0 \) and \( u'(L)=E \):

 
$$ (u'v') + (bu',v) = (f,v) + Ev(L), \quad\forall v\in V$$

First-order derivative in the equation and boundary condition; observations

 
$$ (u'v') + (bu',v) = (f,v) + Ev(L), \quad\forall v\in V,$$

Important:

• The boundary term can be used to implement Neumann conditions
• Forgetting the boundary term implies the condition \( u'=0 \) (!)
• Such conditions are called natural boundary conditions

First-order derivative in the equation and boundary condition; abstract notation

Abstract notation:

 
$$ a(u,v)=L(v)\quad\forall v\in V$$

Here:

 
$$ \begin{align*} a(u,v)&=(u',v') + (bu',v)\\ L(v)&= (f,v) + E v(L) \end{align*} $$

First-order derivative in the equation and boundary condition; linear system

Insert \( u=C+\sum_jc_j\baspsi_j \) and \( v=\baspsi_i \):

 
$$ \sum_{j\in\If} \underbrace{((\baspsi_j',\baspsi_i') + (b\baspsi_j',\baspsi_i))}_{A_{i,j}} c_j = \underbrace{(f,\baspsi_i) + E \baspsi_i(L)}_{b_i} $$

Observation: \( A_{i,j} \) is not symmetric because of the term

 
$$ (b\baspsi_j',\baspsi_i)=\int_{\Omega} b\baspsi_j'\baspsi_i dx \neq \int_{\Omega} b \baspsi_i' \baspsi_jdx = (\baspsi_i',b\baspsi_j) $$

Terminology: natural and essential boundary conditions

 
$$ (u',v') + (bu',v) = (f,v) + u'(L)v(L) - u'(0)v(0)$$

• Note: forgetting the boundary terms implies \( u'(L)=u'(0)=0 \) (unless prescribe a Dirichlet condition)
• Conditions on \( u' \) are simply inserted in the variational form and called natural conditions
• Conditions on \( u \) at \( x=0 \) requires modifying \( V \) (through \( \baspsi_i(0)=0 \)) and are known as essential conditions

Nonlinear coefficient; problem

Problem:

 
$$ \begin{equation} -(\dfc(u)u')' = f(u),\quad x\in [0,L],\ u(0)=0,\ u'(L)=E \end{equation} $$

• \( V \): basis \( \sequencei{\baspsi} \) with \( \baspsi_i(0)=0 \) because of \( u(0)=0 \)
• How does the nonlinear coefficients \( \dfc(u) \) and \( f(u) \) impact the variational formulation?
• (Not much!)

Nonlinear coefficient; variational formulation

Galerkin: multiply by \( v \), integrate, integrate by parts

 
$$ \int_0^L \dfc(u)\frac{du}{dx}\frac{dv}{dx}\dx = \int_0^L f(u)v\dx + [\dfc(u)vu']_0^L\quad\forall v\in V $$

• \( \dfc(u(0))v(0)u'(0)=0 \) since \( v(0) \)
• \( \dfc(u(L))v(L)u'(L) = \dfc(u(L))v(L)E \) since \( u'(L)=E \)

 
$$ \int_0^L \dfc(u)\frac{du}{dx}\frac{dv}{dx}v\dx = \int_0^L f(u)v\dx + \dfc(u(L))v(L)E\quad\forall v\in V $$

or

 
$$ (\dfc(u)u', v') = (f(u),v) + \dfc(u(L))v(L)E\quad\forall v\in V $$

Nonlinear coefficient; where does the nonlinearity cause challenges?

• Abstract notation: no \( a(u,v) \) and \( L(v) \) because \( a \) and \( L \) are nonlinear
• Instead: \( F(u;v)=0\ \forall v\in V \)
• What about forming a linear system? We get a nonlinear system of algebraic equations
• Must use methods like Picard iteration or Newton's method to solve nonlinear algebraic equations
• But: the variational formulation was not much affected by nonlinearities

Computing with Dirichlet and Neumann conditions; problem

 
$$ \begin{equation*} -u''(x)=f(x),\quad x\in \Omega=[0,1],\quad u'(0)=C,\ u(1)=D \end{equation*} $$

• Use a global polynomial basis \( \baspsi_i\sim x^i \) on \( [0,1] \)
• Because of \( u(1)=D \): \( \baspsi_i(1)=0 \)
• Basis: \( \baspsi_i(x)=(1-x)^{i+1} \), \( i\in\If \)
• \( B(x)=Dx \)

Computing with Dirichlet and Neumann conditions; details

 
$$ A_{i,j} = (\baspsi_j',\baspsi_i') = \int_{0}^1 \baspsi_i'(x)\baspsi_j'(x)dx = \int_0^1 (i+1)(j+1)(1-x)^{i+j} dx, $$

Choose \( f(x)=2 \):

 
$$ \begin{align*} b_i &= (2,\baspsi_i) - (D,\baspsi_i') -C\baspsi_i(0)\\ &= \int_0^1 \left( 2(1-x)^{i+1} - D(i+1)(1-x)^i\right)dx -C\baspsi_i(0) \end{align*} $$

Can easily do the integrals with sympy. \( N=1 \):

 
$$ \begin{equation*} \left(\begin{array}{cc} 1 & 1\\ 1 & 4/3 \end{array}\right) \left(\begin{array}{c} c_0\\ c_1 \end{array}\right) = \left(\begin{array}{c} -C+D+1\\ 2/3 -C + D \end{array}\right) \end{equation*} $$

 
$$ c_0=-C+D+2, \quad c_1=-1,$$

 
$$ u(x) = 1 -x^2 + D + C(x-1)\quad\hbox{(exact solution)} $$

When the numerical method is exact

Assume that apart from boundary conditions, \( \uex \) lies in the same space \( V \) as where we seek \( u \):

 
$$ \begin{align*} u &= B + F,\quad F\in V a(B+F, v) &= L(v)\quad\forall v\in V \uex & = B + E,\quad E\in V a(B+E, v) &= L(v)\quad\forall v\in V \end{align*} $$

Subtract: \( a(F-E,v)=0\ \Rightarrow\ E=F \) and \( u = \uex \)

Computing with finite elements

Tasks:

• Address the model problem \( -u''(x)=2 \), \( u(0)=u(L)=0 \)
• Uniform finite element mesh with P1 elements
• Show all finite element computations in detail

Variational formulation, finite element mesh, and basis

 
$$ -u''(x) = 2,\quad x\in (0,L),\ u(0)=u(L)=0,$$

Variational formulation:

 
$$ (u',v') = (2,v)\quad\forall v\in V $$

Since \( u(0)=0 \) and \( u(L)=0 \), we must force

 
$$ v(0)=v(L)=0,\quad \baspsi_i(0)=\baspsi_i(L)=0$$

Use finite element basis, but exclude \( \basphi_0 \) and \( \basphi_{N_n} \) since these are not 0 on the boundary:

 
$$ \baspsi_i=\basphi_{i+1},\quad i=0,\ldots,N=N_n-2$$

Introduce index mapping \( \nu(j) \): \( \baspsi_i = \basphi_{\nu(i)} \)

 
$$ u = \sum_{j\in\If}c_j\basphi_{\nu(i)},\quad i=0,\ldots,N,\quad \nu(j) = j+1$$

Irregular numbering: more complicated \( \nu(j) \) table

Computation in the global physical domain; formulas

 
$$ \begin{equation*} A_{i,j}=\int_0^L\basphi_{i+1}'(x)\basphi_{j+1}'(x) dx,\quad b_i=\int_0^L2\basphi_{i+1}(x) dx \end{equation*} $$

Many will prefer to change indices to obtain a \( \basphi_i'\basphi_j' \) product: \( i+1\rightarrow i \), \( j+1\rightarrow j \)

 
$$ \begin{equation*} A_{i-1,j-1}=\int_0^L\basphi_{i}'(x)\basphi_{j}'(x) \dx,\quad b_{i-1}=\int_0^L2\basphi_{i}(x) \dx \end{equation*} $$

Computation in the global physical domain; details

 
$$ \basphi_i = \pm h^{-1} $$

 
$$ A_{i-1,i-1} = h^{-2}2h = 2h^{-1},\quad A_{i-1,i-2} = h^{-1}(-h^{-1})h = -h^{-1},\quad A_{i-1,i}=A_{i-1,i-2}$$

 
$$ b_{i-1} = 2({\half}h + {\half}h) = 2h$$

Computation in the global physical domain; linear system

 
$$ \begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} 2 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & -1 & 2 \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} 2h \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ 2h \end{array} \right) \tag{53} \end{equation} $$

Comparison with a finite difference discretization

• Recall: \( c_i = u(\xno{i+1})\equiv u_{i+1} \)
• Write out a general equation at node \( i-1 \), expressed by \( u_i \)

 
$$ \begin{equation} -\frac{1}{h}u_{i-1} + \frac{2}{h}u_{i} - \frac{1}{h}u_{i+1} = 2h \tag{54} \end{equation} $$

The standard finite difference method for \( -u''=2 \) is

 
$$ -\frac{1}{h^2}u_{i-1} + \frac{2}{h^2}u_{i} - \frac{1}{h^2}u_{i+1} = 2 $$

(Remains to study the equations involving boundary values)

Cellwise computations; formulas

• Repeat the previous example, but apply the cellwise algorithm
• Work with one cell at a time
• Transform physical cell to reference cell \( X\in [-1,1] \)

 
$$ \begin{equation*} A_{i-1,j-1}^{(e)}=\int_{\Omega^{(e)}} \basphi_i'(x)\basphi_j'(x) \dx = \int_{-1}^1 \frac{d}{dx}\refphi_r(X)\frac{d}{dx}\refphi_s(X) \frac{h}{2} \dX, \end{equation*} $$

 
$$ \refphi_0(X)=\half(1-X),\quad\refphi_1(X)=\half(1+X)$$

 
$$ \frac{d\refphi_0}{dX} = -\half,\quad \frac{d\refphi_1}{dX} = \half $$

From the chain rule

 
$$ \frac{d\refphi_r}{dx} = \frac{d\refphi_r}{dX}\frac{dX}{dx} = \frac{2}{h}\frac{d\refphi_r}{dX}$$

Cellwise computations; details

 
$$ \begin{equation*} A_{i-1,j-1}^{(e)}=\int_{\Omega^{(e)}} \basphi_i'(x)\basphi_j'(x) \dx = \int_{-1}^1 \frac{2}{h}\frac{d\refphi_r}{dX}\frac{2}{h}\frac{d\refphi_s}{dX} \frac{h}{2} \dX = \tilde A_{r,s}^{(e)} \end{equation*} $$

 
$$ \begin{equation*} b_{i-1}^{(e)} = \int_{\Omega^{(e)}} 2\basphi_i(x) \dx = \int_{-1}^12\refphi_r(X)\frac{h}{2} \dX = \tilde b_{r}^{(e)}, \quad i=q(e,r),\ r=0,1 \end{equation*} $$

Must run through all \( r,s=0,1 \) and \( r=0,1 \) and compute each entry in the element matrix and vector:

 
$$ \begin{equation} \tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{rr} 1 & -1\\ -1 & 1 \end{array}\right),\quad \tilde b^{(e)} = h\left(\begin{array}{c} 1\\ 1 \end{array}\right)\tp \tag{55} \end{equation} $$

Example:

 
$$ \tilde A^{(e)}_{0,1} = \int_{-1}^1 \frac{2}{h}\frac{d\refphi_0}{dX}\frac{2}{h}\frac{d\refphi_1}{dX} \frac{h}{2} \dX = \frac{2}{h}(-\half)\frac{2}{h}\half\frac{h}{2} \int_{-1}^1\dX = -\frac{1}{h} $$

Cellwise computations; details of boundary cells

• The boundary cells involve only one unknown
• \( \Omega^{(0)} \): left node value known, only a contribution from right node
• \( \Omega^{(N_e)} \): right node value known, only a contribution from left node

For \( e=0 \) and \( =N_e \):

 
$$ \tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{r} 1 \end{array}\right),\quad \tilde b^{(e)} = h\left(\begin{array}{c} 1 \end{array}\right) $$

Only one degree of freedom ("node") in these cells (\( r=0 \) counts the only dof)

Cellwise computations; assembly

4 P1 elements:

vertices = [0, 0.5, 1, 1.5, 2]
cells = [[0, 1], [1, 2], [2, 3], [3, 4]]
dof_map = [[0], [0, 1], [1, 2], [2]]       # only 1 dof in elm 0, 3

Python code for the assembly algorithm:

# Ae[e][r,s]: element matrix, be[e][r]: element vector
# A[i,j]: coefficient matrix, b[i]: right-hand side

for e in range(len(Ae)):
    for r in range(Ae[e].shape[0]):
        for s in range(Ae[e].shape[1]):
            A[dof_map[e,r],dof_map[e,s]] += Ae[e][i,j]
        b[dof_map[e,r]] += be[e][i,j]

Result: same linear system as arose from computations in the physical domain

General construction of a boundary function

• Now we address nonzero Dirichlet conditions
• \( B(x) \) is not always easy to construct (extend to the interior of \( \Omega \)), especially not in 2D and 3D
• With finite element \( \basphi_i \), \( B(x) \) can be constructed in a completely general way
• \( \Ifb \): set of indices with nodes where \( u \) is known
• \( U_i \): Dirichlet value of \( u \) at node \( i \), \( i\in\Ifb \)

 
$$ \begin{equation} B(x) = \sum_{j\in\Ifb} U_j\basphi_j(x) \end{equation} $$

Suppose we have a Dirichlet condition \( u(\xno{k})=U_k \), \( k\in\Ifb \):

 
$$ u(\xno{k}) = \sum_{j\in\Ifb} U_j\underbrace{\basphi_j(x)}_{\neq 0 \hbox{ only for }j=k} + \sum_{j\in\If} c_j\underbrace{\basphi_{\nu(j)}(\xno{k})}_{=0,\ k\not\in\If} = U_k $$

Example with two Dirichlet values; variational formulation

 
$$ -u''=2, \quad u(0)=C,\ u(L)=D $$

 
$$ \int_0^L u'v'\dx = \int_0^L2v\dx\quad\forall v\in V$$

 
$$ (u',v') = (2,v)\quad\forall v\in V$$

Example with two Dirichlet values; boundary function

 
$$ \begin{equation} B(x) = \sum_{j\in\Ifb} U_j\basphi_j(x) \end{equation} $$

Here \( \Ifb = \{0,N_n\} \), \( U_0=C \), \( U_{N_n}=D \),

 
$$ \baspsi_i = \basphi_{\nu(i)}, \quad \nu(i)=i+1,\quad i\in\If = \{0,\ldots,N=N_n-2\} $$

 
$$ \begin{equation} u(x) = C\basphi_0(x) + D\basphi_{N_n}(x) + \sum_{j\in\If}c_j\basphi_{\nu(j)} \end{equation} $$

Example with two Dirichlet values; details

Insert \( u = B + \sum_j c_j\baspsi_j \) in variational formulation:

 
$$ (u',v') = (2,v)\quad\Rightarrow\quad (\sum_jc_j\baspsi_j',\baspsi_i') = (2-B',\baspsi_i)\quad \forall v\in V$$

 
$$ \begin{align*} u(x) &= \underbrace{C\cdot\basphi_0 + D\basphi_{N_n}}_{B(x)} + \sum_{j\in\If} c_j\basphi_{j+1}\\ &= C\cdot\basphi_0 + D\basphi_{N_n} + c_0\basphi_1 + c_1\basphi_2 +\cdots + c_N\basphi_{N_n-1} \end{align*} $$

 
$$ A_{i-1,j-1} = \int_0^L \basphi_i'(x)\basphi_j'(x) \dx,\quad b_{i-1} = \int_0^L (f(x) - C\basphi_{0}'(x) - D\basphi_{N_n}'(x)) \basphi_i(x) \dx $$

 
for \( i,j = 1,\ldots,N+1=N_n-1 \).

New boundary terms from \( -\int B'\basphi_i\dx \): \( C/2 \) for \( i=1 \) and \( -D/2 \) for \( i=N_n-1 \)

Example with two Dirichlet values; cellwise computations

• Element matrices as in the previous example (with \( u=0 \) on the boundary)
• New element vector in the first and last cell

From the last cell:

 
$$ \tilde b_0^{(N_e)} = \int_{-1}^1 \left(f - D\frac{2}{h} \frac{d\refphi_1}{dX}\right) \refphi_0\frac{h}{2} \dX = (\frac{h}{2}(2 - D\frac{2}{h}\half) \int_{-1}^1 \refphi_0 \dX = h - D/2 $$

From the first cell:

 
$$ \tilde b_0^{(0)} = \int_{-1}^1 \left(f - C\frac{2}{h} \frac{d\refphi_0}{dX}\right) \refphi_1\frac{h}{2} \dX = (\frac{h}{2}(2 + C\frac{2}{h}\half) \int_{-1}^1 \refphi_1 \dX = h + C/2\tp $$

Modification of the linear system; ideas

• Method 1: incorporate Dirichlet values through a \( B(x) \) function and demand \( \baspsi_i=0 \) where Dirichlet values apply
• Method 2: drop \( B(x) \), drop demands to \( \baspsi_i \), just assemble as if there were no Dirichlet conditions, and modify the linear system instead

Method 2: always \( \baspsi_i = \basphi_i \) and

 
$$ \begin{equation} u(x) = \sum_{j\in\If}c_j\basphi_j(x),\quad \If=\{0,\ldots,N=N_n\} \tag{56} \end{equation} $$

Modification of the linear system; original system

 
$$ -u''=2,\quad u(0)=0,\ u(L)=D$$

Assemble as if there were no Dirichlet conditions:

 
$$ \begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} 1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & -1 & 1 \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} h \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h\\ h \end{array} \right) \tag{57} \end{equation} $$

Modification of the linear system; row replacement

• Dirichlet condition \( u(\xno{k})= U_k \) means \( c_k=U_k \) (since \( c_k=u(\xno{k}) \))
• Replace first row by \( c_0=0 \)
• Replace last row by \( c_N=D \)

 
$$ \begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} h & 0 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & h \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} 0 \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h\\ D \end{array} \right) \tag{58} \end{equation} $$

Modification of the linear system; element matrix/vector

In cell 0 we know \( u \) for local node (degree of freedom) \( r=0 \). Replace the first cell equation by \( \tilde c_0 = 0 \):

 
$$ \begin{equation} \tilde A^{(0)} = A = \frac{1}{h}\left(\begin{array}{rr} h & 0\\ -1 & 1 \end{array}\right),\quad \tilde b^{(0)} = \left(\begin{array}{c} 0\\ h \end{array}\right) \tag{59} \end{equation} $$

In cell \( N_e \) we know \( u \) for local node \( r=1 \). Replace the last equation in the cell system by \( \tilde c_1=D \):

 
$$ \begin{equation} \tilde A^{(N_e)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & -1\\ 0 & h \end{array}\right),\quad \tilde b^{(N_e)} = \left(\begin{array}{c} h\\ D \end{array}\right) \tag{60} \end{equation} $$

Symmetric modification of the linear system; algorithm

• The modification above destroys symmetry of the matrix: e.g., \( A_{0,1}\neq A_{1,0} \)
• Symmetry is often important in 2D and 3D (faster computations)
• A more complex modification can preserve symmetry!

Algorithm for incorporating \( c_i=U_i \) in a symmetric way:

1. Subtract column \( i \) times \( U_i \) from the right-hand side
2. Zero out column and row no \( i \)
3. Place 1 on the diagonal
4. Set \( b_i=U_i \)

Symmetric modification of the linear system; example

 
$$ \begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} 1 & 0 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & 0 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & 1 \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} 0 \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h +D/h\\ D \end{array} \right) \tag{61} \end{equation} $$

Symmetric modification of the linear system; element level

Symmetric modification applied to \( \tilde A^{(N_e)} \):

 
$$ \begin{equation} \tilde A^{(N_e)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & 0\\ 0 & 1 \end{array}\right),\quad \tilde b^{(N-1)} = \left(\begin{array}{c} h + D/h\\ D \end{array}\right) \tag{62} \end{equation} $$

Boundary conditions: specified derivative

 
$$ -u''=f,\quad u'(0)=C,\ u(L)=D$$

• \( \baspsi_i(L)=0 \) because of Dirichlet condition \( u(L)=D \)
• No demand to \( \baspsi_i(0) \)

The variational formulation

Galerkin's method:

 
$$ \begin{equation*} \int_0^L(u''(x)+f(x))\baspsi_i(x) dx = 0,\quad i\in\If \end{equation*} $$

Integration of \( u''\baspsi_i \) by parts:

 
$$ \begin{equation*} \int_0^Lu'(x)\baspsi_i'(x) \dx -(u'(L)\baspsi_i(L) - u'(0)\baspsi_i(0)) - \int_0^L f(x)\baspsi_i(x) \dx =0, \quad i\in\If \end{equation*} $$

• \( u'(L){\baspsi_i(L)}=0 \) since \( \baspsi_i(L)=0 \)
• \( u'(0)\baspsi_i(0) = C\baspsi_i(0) \) since \( u'(0)=C \)

Method 1: Boundary function and exclusion of Dirichlet degrees of freedom

• \( \baspsi_i = \basphi_i \), \( i\in\If =\{0,\ldots,N=N_n-1\} \)
• \( B(x)=D\basphi_{N_n}(x) \), \( u= B + \sum_{j=0}^N c_j\basphi_j \)

 
$$ \begin{equation*} \int_0^Lu'(x)\basphi_i'(x) dx = \int_0^L f(x)\basphi_i(x) dx - C\basphi_i(0),\quad i\in\If \end{equation*} $$

 
$$ \begin{equation} \sum_{j=0}^{N=N_n-1}\left( \int_0^L \basphi_i'(x)\basphi_j'(x) dx \right)c_j = \int_0^L\left(f(x)\basphi_i(x) -D\basphi_N'(x)\basphi_i(x)\right) dx - C\basphi_i(0) \tag{63} \end{equation} $$

 
for \( i=0,\ldots,N=N_n-1 \).

Method 2: Use all \( \basphi_i \) and insert the Dirichlet condition in the linear system

• Now \( \baspsi_i=\basphi_i \), \( i=0,\ldots,N=N_n \)
• \( \basphi_N(L)\neq 0 \), so \( u'(L)\basphi_N(L)\neq 0 \)
• However, the term \( u'(L)\basphi_N(L) \) in \( b_N \) will be erased when we insert the Dirichlet value in \( b_N=D \)

We can forget about the term \( u'(L)\basphi_i(L) \)!

 
$$ \begin{equation*} u(x) = \sum_{j=0}^{N=N_n} c_j\basphi_j(x) \end{equation*} $$

 
$$ \begin{equation} \sum_{j=0}^{N=N_n}\left( \int_0^L \basphi_i'(x)\basphi_j'(x) dx \right)c_j = \int_0^L f(x)\basphi_i(x)\basphi_i(x) dx - C\basphi_i(0) \tag{63} \end{equation} $$

Assemble entries for \( i=0,\ldots,N=N_n \) and then modify the last equation to \( c_N=D \)

How the Neumann condition impacts the element matrix and vector

The extra term \( C\basphi_0(0) \) affects only the element vector from the first cells since \( \basphi_0=0 \) on all other cells.

 
$$ \begin{equation} \tilde A^{(0)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & 1\\ -1 & 1 \end{array}\right),\quad \tilde b^{(0)} = \left(\begin{array}{c} h - C\\ h \end{array}\right) \tag{65} \end{equation} $$

The finite element algorithm

The differential equation problem defines the integrals in the variational formulation.

Request these functions from the user:

integrand_lhs(phi, r, s, x)
boundary_lhs(phi, r, s, x)
integrand_rhs(phi, r, x)
boundary_rhs(phi, r, x)

Must also have a mesh with vertices, cells, and dof_map

Python pseudo code; the element matrix and vector

<Declare global matrix, global rhs: A, b>

# Loop over all cells
for e in range(len(cells)):

    # Compute element matrix and vector
    n = len(dof_map[e])  # no of dofs in this element
    h = vertices[cells[e][1]] - vertices[cells[e][0]]
    <Declare element matrix, element vector: A_e, b_e>

    # Integrate over the reference cell
    points, weights = <numerical integration rule>
    for X, w in zip(points, weights):
        phi = <basis functions + derivatives at X>
        detJ = h/2
        x = <affine mapping from X>
        for r in range(n):
            for s in range(n):
                A_e[r,s] += integrand_lhs(phi, r, s, x)*detJ*w
            b_e[r] += integrand_rhs(phi, r, x)*detJ*w

    # Add boundary terms
    for r in range(n):
        for s in range(n):
            A_e[r,s] += boundary_lhs(phi, r, s, x)*detJ*w
        b_e[r] += boundary_rhs(phi, r, x)*detJ*w

Python pseudo code; boundary conditions and assembly

for e in range(len(cells)):
    ...

    # Incorporate essential boundary conditions
    for r in range(n):
        global_dof = dof_map[e][r]
        if global_dof in essbc_dofs:
            # dof r is subject to an essential condition
            value = essbc_docs[global_dof]
            # Symmetric modification
            b_e -= value*A_e[:,r]
            A_e[r,:] = 0
            A_e[:,r] = 0
            A_e[r,r] = 1
            b_e[r] = value

    # Assemble
    for r in range(n):
        for s in range(n):
            A[dof_map[e][r], dof_map[e][r]] += A_e[r,s]
        b[dof_map[e][r] += b_e[r]

<solve linear system>

Variational formulations in 2D and 3D

Integration by parts

• \( \partial\Omega_N \): Neumann conditions \( -a\frac{\partial u}{\partial n} = g \)
• \( \partial\Omega_D \): Dirichlet conditions \( u = u_0 \)
• \( v\in V \) must vanish on \( \partial\Omega_D \) (in method 1)

Example on integration by parts; problem

 
$$ \begin{align} \v\cdot\nabla u + \alpha u &= \nabla\cdot\left( a\nabla u\right) + f, \quad & \x\in\Omega\\ u &= u_0,\quad &\x\in\partial\Omega_D\\ -a\frac{\partial u}{\partial n} &= g,\quad &\x\in\partial\Omega_N \end{align} $$

• Known: \( a \), \( \alpha \), \( f \), \( u_0 \), and \( g \).
• Second-order PDE: must have exactly one boundary condition at each point of the boundary

Method 1 with boundary function and \( \baspsi_i=0 \) on \( \partial\Omega_D \):

 
$$ u(\x) = B(\x) + \sum_{j\in\If} c_j\baspsi_j(\x),\quad B(\x)=u_0(\x) $$

Example on integration by parts; details (1)

Galerkin's method: multiply by \( v\in V \) and integrate over \( \Omega \),

 
$$ \int_{\Omega} (\v\cdot\nabla u + \alpha u)v\dx = \int_{\Omega} \nabla\cdot\left( a\nabla u\right)\dx + \int_{\Omega}fv \dx $$

Integrate second-order term by parts:

 
$$ \int_{\Omega} \nabla\cdot\left( a\nabla u\right) v \dx = -\int_{\Omega} a\nabla u\cdot\nabla v\dx + \int_{\partial\Omega} a\frac{\partial u}{\partial n} v\ds, $$

Resulting variational form:

 
$$ \int_{\Omega} (\v\cdot\nabla u + \alpha u)v\dx = -\int_{\Omega} a\nabla u\cdot\nabla v\dx + \int_{\partial\Omega} a\frac{\partial u}{\partial n} v\ds + \int_{\Omega} fv \dx $$

Example on integration by parts; details (2)

Note: \( v\neq 0 \) only on \( \partial\Omega_N \):

 
$$ \int_{\partial\Omega} a\frac{\partial u}{\partial n} v\ds = \int_{\partial\Omega_N} \underbrace{a\frac{\partial u}{\partial n}}_{-g} v\ds = -\int_{\partial\Omega_N} gv\ds $$

The final variational form:

 
$$ \int_{\Omega} (\v\cdot\nabla u + \alpha u)v\dx = -\int_{\Omega} a\nabla u\cdot\nabla v \dx - \int_{\partial\Omega_N} g v\ds + \int_{\Omega} fv \dx $$

Or with inner product notation:

 
$$ (\v\cdot\nabla u, v) + (\alpha u,v) = - (a\nabla u,\nabla v) - (g,v)_{N} + (f,v) $$

\( (g,v)_{N} \): line or surface integral over \( \partial\Omega_N \).

Example on integration by parts; linear system

 
$$ u = B + \sum_{j\in\If} c_j\baspsi_j,\quad B = u_0 $$

 
$$ A_{i,j} = (\v\cdot\nabla \baspsi_j, \baspsi_i) + (\alpha \baspsi_j ,\baspsi_i) + (a\nabla \baspsi_j,\nabla \baspsi_i) $$

 
$$ b_i = (g,\baspsi_i)_{N} + (f,\baspsi_i) - (\v\cdot\nabla u_0, \baspsi_i) + (\alpha u_0 ,\baspsi_i) + (a\nabla u_0,\nabla \baspsi_i) $$

Transformation to a reference cell in 2D/3D (1)

 
$$ \begin{equation} \int_{{\Omega}^{(e)}} a(\x)\nabla\basphi_i\cdot\nabla\basphi_j\dx \end{equation} $$

Mapping from reference to physical coordinates:

 
$$ \x(\X) $$

with Jacobian \( J \),

 
$$ J_{i,j}=\frac{\partial x_j}{\partial X_i} $$

• \( \dx \rightarrow \det J\dX \).
• Must express \( \nabla\basphi_i \) by an expression with \( \refphi_r \), \( i=q(e,r) \): \( \nabla\refphi_r(\X) \)
• We want \( \nabla_{\x}\refphi_r(\X) \) (derivatives wrt \( \x \))
• What we readily have is \( \nabla_{\X}\refphi_r(\X) \) (derivative wrt \( \X \))
• Need to transform \( \nabla_{\X}\refphi_r(\X) \) to \( \nabla_{\x}\refphi_r(\X) \)

Transformation to a reference cell in 2D/3D (2)

Can derive

 
$$ \begin{align*} \nabla_{\X}\refphi_r &= J\cdot\nabla_{\x}\basphi_i\\ \nabla_{\x}\basphi_i &= \nabla_{\x}\refphi_r(\X) = J^{-1}\cdot\nabla_{\X}\refphi_r(\X) \end{align*} $$

Integral transformation from physical to reference coordinates:

 
$$ \begin{equation} \int_{\Omega^{(e)}} a(\x)\nabla_{\x}\basphi_i\cdot\nabla_{\x}\basphi_j\dx = \int_{\tilde\Omega^r} a(\x(\X))(J^{-1}\cdot\nabla_{\X}\refphi_r)\cdot (J^{-1}\cdot\nabla\refphi_s)\det J\dX \end{equation} $$

Numerical integration

Numerical integration over reference cell triangles and tetrahedra:

 
$$ \int_{\tilde\Omega^r} g\dX = \sum_{j=0}^{n-1} w_j g(\bar\X_j)$$

Module numint.py contains different rules:

>>> import numint
>>> x, w = numint.quadrature_for_triangles(num_points=3)
>>> x
[(0.16666666666666666, 0.16666666666666666),
 (0.66666666666666666, 0.16666666666666666),
 (0.16666666666666666, 0.66666666666666666)]
>>> w
[0.16666666666666666, 0.16666666666666666, 0.16666666666666666]

• Triangle: rules with \( n=1,3,4,7 \) integrate exactly polynomials of degree \( 1,2,3,4 \), resp.
• Tetrahedron: rules with \( n=1,4,5,11 \) integrate exactly polynomials of degree \( 1,2,3,4 \), resp.

Time-dependent problems

• So far: used the finite element framework for discretizing in space
• What about \( u_t = u_{xx} + f \)?

1. Use finite differences in time to obtain a set of recursive spatial problems
2. Solve the spatial problems by the finite element method

Example: diffusion problem

 
$$ \begin{align} \frac{\partial u}{\partial t} &= \dfc\nabla^2 u + f(\x, t),\quad &\x\in\Omega, t\in (0,T] \tag{67}\\ u(\x, 0) & = I(\x),\quad &\x\in\Omega \tag{68}\\ \frac{\partial u}{\partial n} &= 0,\quad &\x\in\partial\Omega,\ t\in (0,T] \tag{69} \end{align} $$

A Forward Euler scheme; ideas

 
$$ \begin{equation} [D_t^+ u = \dfc\nabla^2 u + f]^n,\quad n=1,2,\ldots,N_t-1 \end{equation} $$

Solving wrt \( u^{n+1} \):

 
$$ \begin{equation} u^{n+1} = u^n + \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right) \tag{70} \end{equation} $$

• \( u^n = \sum_jc_j^n\baspsi_j\, \in V \), \( u^{n+1} = \sum_jc_j^{n+1}\baspsi_j\,\in V \)
• Compute \( u^0 \) from \( I \)
• Compute \( u^{n+1} \) from \( u^n \) by solving the PDE for \( u^{n+1} \) at each time level

A Forward Euler scheme; stages in the discretization

• \( \uex(\x,t) \): exact solution of the space-and time-continuous problem
• \( \uex^n(\x) \): exact solution of time-discrete problem (after applying a finite difference scheme in time)
• \( \uex^n(\x)\approx u^n = \sum_{j\in\If}c_j^n\baspsi_j = \) solution of the time- and space-discrete problem (after applying a Galerkin method in space)

 
$$ \begin{equation} \frac{\partial \uex}{\partial t} = \dfc\nabla^2 \uex + f(\x, t) \tag{71} \end{equation} $$

 
$$ \begin{equation} \uex^{n+1} = \uex^n + \Delta t \left( \dfc\nabla^2 \uex^n + f(\x, t_n)\right) \tag{72} \end{equation} $$

 
$$ \uex^n \approx u^n = \sum_{j=0}^{N} c_j^{n}\baspsi_j(\x),\quad \uex^{n+1} \approx u^{n+1} = \sum_{j=0}^{N} c_j^{n+1}\baspsi_j(\x) $$

 
$$ R = u^{n+1} - u^n - \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right)$$

A Forward Euler scheme; weighted residual (or Galerkin) principle

 
$$ R = u^{n+1} - u^n - \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right)$$

The weighted residual principle:

 
$$ \int_\Omega Rw\dx = 0,\quad \forall w\in W$$

results in

 
$$ \int_\Omega \left\lbrack u^{n+1} - u^n - \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right) \right\rbrack w \dx =0, \quad \forall w \in W $$

Galerkin: \( W=V \), \( w=v \)

A Forward Euler scheme; integration by parts

Isolating the unknown \( u^{n+1} \) on the left-hand side:

 
$$ \int_{\Omega} u^{n+1}\baspsi_i\dx = \int_{\Omega} \left\lbrack u^n - \Delta t \left( \dfc\nabla^2 u^n + f(\x, t_n)\right) \right\rbrack v\dx $$

Integration by parts of \( \int\dfc(\nabla^2 u^n) v\dx \):

 
$$ \int_{\Omega}\dfc(\nabla^2 u^n)v \dx = -\int_{\Omega}\dfc\nabla u^n\cdot\nabla v\dx + \underbrace{\int_{\partial\Omega}\dfc\frac{\partial u^n}{\partial n}v \dx}_{=0\quad\Leftarrow\quad\partial u^n/\partial n=0} $$

Variational form:

 
$$ \begin{equation} \int_{\Omega} u^{n+1} v\dx = \int_{\Omega} u^n v\dx - \Delta t \int_{\Omega}\dfc\nabla u^n\cdot\nabla v\dx + \Delta t\int_{\Omega}f^n v\dx,\quad\forall v\in V \tag{73} \end{equation} $$

New notation for the solution at the most recent time levels

• \( u \) and u: the spatial unknown function to be computed
• \( u_1 \) and u_1: the spatial function at the previous time level \( t-\Delta t \)
• \( u_2 \) and u_2: the spatial function at \( t-2\Delta t \)
• This new notation gives close correspondance between code and math

 
$$ \begin{equation} \int_{\Omega} u v\dx = \int_{\Omega} u_1 v\dx - \Delta t \int_{\Omega}\dfc\nabla u_1\cdot\nabla v\dx + \Delta t\int_{\Omega}f^n v\dx \tag{74} \end{equation} $$

or shorter

 
$$ \begin{equation} (u, \baspsi_i) = (u_1,v) - \Delta t (\dfc\nabla u_1,\nabla v) + (f^n, v) \tag{75} \end{equation} $$

Deriving the linear systems

• \( u = \sum_{j=0}^{N}c_j\baspsi_j(\x) \)
• \( u_1 = \sum_{j=0}^{N} c_{1,j}\baspsi_j(\x) \)
• \( \forall v\in V \): for \( v=\baspsi_i \), \( i=0,\ldots,N \)

 
$$

Insert these in

$$

 
(u, \baspsi_i) = (u_1,\baspsi_i) - \Delta t (\dfc\nabla u_1,\nabla\baspsi_i) + (f^n,\baspsi_i)

 
$$

and order terms as matrix-vector products:

$$

 
\begin{equation} \sum_{j=0}^{N} \underbrace{(\baspsi_i,\baspsi_j)}_{M_{i,j}} c_j = \sum_{j=0}^{N} \underbrace{(\baspsi_i,\baspsi_j)}_{M_{i,j}} c_{1,j} -\Delta t \sum_{j=0}^{N} \underbrace{(\nabla\baspsi_i,\dfc\nabla\baspsi_j)}_{K_{i,j}} c_{1,j} + (f^n,\baspsi_i),\quad i=0,\ldots,N \end{equation} $$

Structure of the linear systems

 
$$ \begin{equation} Mc = Mc_1 - \Delta t Kc_1 + f \end{equation} $$

 
$$ \begin{align*} M &= \{M_{i,j}\},\quad M_{i,j}=(\baspsi_i,\baspsi_j),\quad i,j\in\If\\ K &= \{K_{i,j}\},\quad K_{i,j}=(\nabla\baspsi_i,\dfc\nabla\baspsi_j),\quad i,j\in\If\\ f &= \{(f(\x,t_n),\baspsi_i)\}_{i\in\If}\\ c &= \{c_i\}_{i\in\If}\\ c_1 &= \{c_{1,i}\}_{i\in\If} \end{align*} $$

Computational algorithm

1. Compute \( M \) and \( K \).
2. Initialize \( u^0 \) by either interpolation or projection
3. For \( n=1,2,\ldots,N_t \):

1. compute \( b = Mc_1 - \Delta t Kc_1 + f \)
2. solve \( Mc = b \)
3. set \( c_1 = c \)

Initial condition:

• Either interpolation: \( c_{1,j} = I(\x_j) \) (finite elements)
• Or projection: solve \( \sum_j M_{i,j}c_{1,j} = (I,\baspsi_i) \), \( i\in\If \)

Comparing P1 elements with the finite difference method; ideas

• P1 elements in 1D
• Uniform mesh on \( [0,L] \) with cell length \( h \)
• No Dirichlet conditions: \( \baspsi_i=\basphi_i \), \( i=0,\ldots,N=N_n \)
• Have found formulas for \( M \) and \( K \) at the element level
• Have assembled the global matrices
• Have developed corresponding finite difference operator formulas
• \( M \): \( h[D_t^+(u + \frac{1}{6}h^2D_xD_x u)]^n_i \)
• \( K \): \( h[\dfc D_xD_x u]^n_i \)

Comparing P1 elements with the finite difference method; results

Diffusion equation with finite elements is equivalent to

 
$$ \begin{equation} [D_t^+(u + \frac{1}{6}h^2D_xD_x u) = \dfc D_xD_x u + f]^n_i \tag{76} \end{equation} $$

Can lump the mass matrix by Trapezoidal integration and get the standard finite difference scheme

 
$$ \begin{equation} [D_t^+u = \dfc D_xD_x u + f]^n_i \end{equation} $$

Discretization in time by a Backward Euler scheme

Backward Euler scheme in time:

 
$$ [D_t^- u = \dfc\nabla^2 u + f(\x, t)]^n \tp $$

 
$$ \begin{equation} \uex^{n} - \Delta t \left( \dfc\nabla^2 \uex^n + f(\x, t_{n})\right) = \uex^{n-1} \tag{77} \end{equation} $$

 
$$ \uex^n \approx u^n = \sum_{j=0}^{N} c_j^{n}\baspsi_j(\x),\quad \uex^{n+1} \approx u^{n+1} = \sum_{j=0}^{N} c_j^{n+1}\baspsi_j(\x)$$

The variational form of the time-discrete problem

 
$$ \begin{equation} \int_{\Omega} \left( u^{n}v + \Delta t \dfc\nabla u^n\cdot\nabla v\right)\dx = \int_{\Omega} u^{n-1} v\dx - \Delta t\int_{\Omega}f^n v\dx,\quad\forall v\in V \tag{78} \end{equation} $$

or

 
$$ \begin{equation} (u,v) + \Delta t (\dfc\nabla u,\nabla v) = (u_1,v) + \Delta t (f^n,\baspsi_i) \tag{79} \end{equation} $$

The linear system: insert \( u=\sum_j c_j\baspsi_i \) and \( u_1=\sum_j c_{1,j}\baspsi_i \),

 
$$ \begin{equation} (M + \Delta t \dfc K)c = Mc_1 + f \tag{80} \end{equation} $$

Calculations with P1 elements in 1D

Can interpret the resulting equation system as

 
$$ \begin{equation} [D_t^-(u + \frac{1}{6}h^2D_xD_x u) = \dfc D_xD_x u + f]^n_i \tag{81} \end{equation} $$

Lumped mass matrix (by Trapezoidal integration) gives a standard finite difference method:

 
$$ \begin{equation} [D_t^- u = \dfc D_xD_x u + f]^n_i \tag{82} \end{equation} $$

Dirichlet boundary conditions

Dirichlet condition at \( x=0 \) and Neumann condition at \( x=L \):

 
$$ \begin{align} u(\x,t) &= u_0(\x,t),\quad & \x\in\partial\Omega_D\\ -\dfc\frac{\partial}{\partial n} u(\x,t) &= g(\x,t),\quad & \x\in\partial{\Omega}_N \end{align} $$

Forward Euler in time, Galerkin's method, and integration by parts:

 
$$ \begin{equation} \int_\Omega u^{n+1}v\dx = \int_\Omega (u^n - \Delta t\dfc\nabla u^n\cdot\nabla v)\dx - \Delta t\int_{\partial\Omega_N} gv\ds,\quad \forall v\in V \end{equation} $$

Requirement: \( v=0 \) on \( \partial\Omega_D \)

Boundary function

 
$$ u^n(\x) = u_0(\x,t_n) + \sum_{j\in\If}c_j^n\baspsi_j(\x)$$

 
$$ \begin{align*} \sum_{j\in\If} \left(\int_\Omega \baspsi_i\baspsi_j\dx\right) c^{n+1}_j &= \sum_{j\in\If} \left(\int_\Omega\left( \baspsi_i\baspsi_j - \Delta t\dfc\nabla \baspsi_i\cdot\nabla\baspsi_j\right)\dx\right) c_j^n - \\ &\quad \int_\Omega\left( u_0(\x,t_{n+1}) - u_0(\x,t_n) + \Delta t\dfc\nabla u_0(\x,t_n)\cdot\nabla \baspsi_i\right)\dx \\ & \quad + \Delta t\int_\Omega f\baspsi_i\dx - \Delta t\int_{\partial\Omega_N} g\baspsi_i\ds, \quad i\in\If \end{align*} $$

Finite element basis functions

• \( B(\x,t_n)=\sum_{j\in\Ifb} U_j^n\basphi_j \)
• \( \baspsi_i = \basphi_{\nu(j)} \), \( j\in\If \)
• \( \nu(j) \), \( j\in\If \), are the node numbers corresponding to all nodes without a Dirichlet condition

 
$$ \begin{align*} u^n &= \sum_{j\in\Ifb} U_j^n\basphi_j + \sum_{j\in\If}c_{1,j}\basphi_{\nu(j)},\\ u^{n+1} &= \sum_{j\in\Ifb} U_j^{n+1}\basphi_j + \sum_{j\in\If}c_{j}\basphi_{\nu(j)} \end{align*} $$

 
$$ \begin{align*} \sum_{j\in\If} \left(\int_\Omega \basphi_i\basphi_j\dx\right) c_j &= \sum_{j\in\If} \left(\int_\Omega\left( \basphi_i\basphi_j - \Delta t\dfc\nabla \basphi_i\cdot\nabla\basphi_j\right)\dx\right) c_{1,j} - \\ &\quad \sum_{j\in\Ifb}\int_\Omega\left( \basphi_i\basphi_j(U_j^{n+1} - U_j^n) + \Delta t\dfc\nabla \basphi_i\cdot\nabla \basphi_jU_j^n\right)\dx \\ &\quad + \Delta t\int_\Omega f\basphi_i\dx - \Delta t\int_{\partial\Omega_N} g\basphi_i\ds, \quad i\in\If \end{align*} $$

Modification of the linear system; the raw system

• Drop boundary function
• Compute as if there are not Dirichlet conditions
• Modify the linear system to incorporate Dirichlet conditions
• \( \If \) holds the indices of all nodes \( \{0,1,\ldots,N=N_n\} \)

 
$$ \begin{align*} \sum_{j\in\If} \biggl(\underbrace{\int_\Omega \basphi_i\basphi_j\dx}_{M_{i,j}}\biggr) c_j &= \sum_{j\in\If} \biggl(\underbrace{\int_\Omega \basphi_i\basphi_j \dx}_{M_{i,j}} - \Delta t\underbrace{\int_\Omega \dfc\nabla \basphi_i\cdot\nabla\basphi_j\dx}_{K_{i,j}}\biggr) c_{1,j} \\ &\quad \underbrace{-\Delta t\int_\Omega f\basphi_i\dx - \Delta t\int_{\partial\Omega_N} g\basphi_i\ds}_{f_i},\quad i\in\If \end{align*} $$

Modification of the linear system; setting Dirichlet conditions

 
$$ \begin{equation} Mc = b,\quad b = Mc_1 - \Delta t Kc_1 + \Delta t f \end{equation} $$

For each \( k \) where a Dirichlet condition applies, \( u(\xno{k},t_{n+1})=U_k^{n+1} \),

• set row \( k \) in \( M \) to zero and 1 on the diagonal: \( M_{k,j}=0 \), \( j\in\If \), \( M_{k,k}=1 \)
• \( b_k = U_k^{n+1} \)

Or apply the slightly more complicated modification which preserves symmetry of \( M \)

Modification of the linear system; Backward Euler example

Backward Euler discretization in time gives a more complicated coefficient matrix:

 
$$ \begin{equation} Ac=b,\quad A = M + \Delta t K,\quad b = Mc_1 + \Delta t f\tp \end{equation} $$

• Set row \( k \) to zero and 1 on the diagonal: \( M_{k,j}=0 \), \( j\in\If \), \( M_{k,k}=1 \)
• Set row \( k \) to zero: \( K_{k,j}=0 \), \( j\in\If \)
• \( b_k = U_k^{n+1} \)

Observe: \( A_{k,k} = M_{k,k} + \Delta t K_{k,k} = 1 + 0 \), so \( c_k = U_k^{n+1} \)

Analysis of the discrete equations

The diffusion equation \( u_t = \dfc u_{xx} \) allows a (Fourier) wave component

 
$$ \begin{equation} u = \Aex^n e^{ikx},\quad \Aex = e^{-\dfc k^2\Delta t} \tag{83} \end{equation} $$

Numerical schemes often allow the similar solution

 
$$ \begin{equation} u^n_q = A^n e^{ikx} \tag{84} \end{equation} $$

• \( A \): amplification factor to be computed
• How good is this \( A \) compared to the exact one?

Handy formulas

 
$$ \begin{align*} [D_t^+ A^n e^{ikq\Delta x}]^n &= A^n e^{ikq\Delta x}\frac{A-1}{\Delta t},\\ [D_t^- A^n e^{ikq\Delta x}]^n &= A^n e^{ikq\Delta x}\frac{1-A^{-1}}{\Delta t},\\ [D_t A^n e^{ikq\Delta x}]^{n+\half} &= A^{n+\half} e^{ikq\Delta x}\frac{A^{\half}-A^{-\half}}{\Delta t} = A^ne^{ikq\Delta x}\frac{A-1}{\Delta t},\\ [D_xD_x A^ne^{ikq\Delta x}]_q &= -A^n \frac{4}{\Delta x^2}\sin^2\left(\frac{k\Delta x}{2}\right)\tp \end{align*} $$

Amplification factor for the Forward Euler method; results

Introduce \( p=k\Delta x/2 \) and \( C=\dfc\Delta t/\Delta x^2 \):

 
$$ A = 1 - 4C\frac{\sin^2 p}{1 - \underbrace{\frac{2}{3}\sin^2 p}_{\hbox{from }M}}$$

(See notes for details)

Stability: \( |A|\leq 1 \):

 
$$ \begin{equation} C\leq \frac{1}{6}\quad\Rightarrow\quad \Delta t\leq \frac{\Delta x^2}{6\dfc} \end{equation} $$

Finite differences: \( C\leq {\half} \), so finite elements give a stricter stability criterion for this PDE!

Amplification factor for the Backward Euler method; results

Coarse meshes:

 
$$ A = \left( 1 + 4C\frac{\sin^2 p}{1 + \frac{2}{3}\sin^2 p}\right)^{-1} \hbox{ (unconditionally stable)} $$

Amplification factors for smaller time steps; Forward Euler

Amplification factors for smaller time steps; Backward Euler