23
Electrostatics via Finite Elements

---

In Chapter 19, we discussed the simple, but powerful, finite-differences method for solving Poisson’s and Laplace’s equations on a lattice in space. In this chapter, we provide a basic outline of the finite-element method (FEM) for solving PDEs. Our usual approach to solving PDEs uses the finite-difference method to approximate various derivatives in terms of the finite differences of a function evaluated upon a fixed grid. The FEM, in contrast, breaks space up into multiple geometric objects (elements), determines an approximate form for the solution appropriate to each element, and then matches the solutions up at the elements′edges. The FEM is ultimately more efficient and powerful than the finite-differences method; however, much more work is required to derive the algorithm. In practice, it is rare to solve a PDE from scratch by deriving the FEM for a particular problem. Rather, and for good reasons, many FEM applications use highly developed FEM packages that get customized for an individual problem. (Python’s finite element library is FiPy.) Our aim is to give the reader some basic understanding of the FEM, not to develop a practitioner. Accordingly, we examine a 1D problem in some detail, and then outline the similar steps followed for the same equation extended to 2D.

---

23.1 Finite-Element Method

The theory and practice of FEM as a numerical method for solving partial differential equations have been developed over the last 30 years and is still an active field of research. One of the theoretical strengths of FEM is that its mathematical foundations allow for elegant proofs of the convergence of its solutions. One of the practical strengths of FEM is that it offers great flexibility for problems on irregular domains, or for problems with highly varying conditions or even singularities. Although finite-difference methods are simpler to implement than FEM, they are less robust mathematically and for big problems less efficient in terms of computer time. Finite elements, in turn, are more complicated to implement, but more appropriate and precise for complicated equations and complicated geometries. In addition, the same basic finite-element technique can be applied to many problems with only minor modifications, and yields solutions that may be evaluated throughout all space, not just on a grid. In fact, the FEM with various preprogrammed multigrid packages has very much become the standard for large-scale engineering applications. Our discussion is based upon Shaw (1992); Li (2014); Otto (2011).

23.2 Electric Field from Charge Density (Problem)

As shown in Figure 23.1, you are given two conducting plates a distance b − a apart, with the lower one kept at potential U_a, the upper plate at potential U_b, and a uniform charge density ρ(x) placed between them. Your problem is to compute the electric potential between the plates.

23.3 Analytic Solution

The relation between charge density ρ(x) and potential U(x) is given by Poisson’s equation (19.6). For our problem, the potential U changes only in the x direction, and so the PDE becomes the ODE:

(23.1)

where we have set ρ(x) = 1/4π to simplify the programming. The solution we want is subject to the Dirichlet boundary conditions:

(23.2)

(23.3)

Although, we know the analytic solution (23.3), we shall develop the FEM for solving the ODE as if it was a PDE (it would be in 2D), and as if we did not know the solution. Although we will not demonstrate it, this method works equally well for any charge density ρ(x).

Figure 23.1 A finite element solution to Laplace’s equation for two metal plates with a charge density between them. The dots are the nodes x_i, and the lines connecting the nodes are the finite elements.

23.4 Finite-Element (Not Difference) Methods, 1D

In an FEM, the domain in which the PDE is solved is split into finite subdomains, called elements, and a trial solution to the PDE in each subdomain is hypothesized. Then the parameters of the trial solution are adjusted to obtain a best fit (in the sense of Chapter 7) to the exact solution. Essentially, this approach converts a given PDE into an integral equation known as the weak or variational form (“weak” because there is no longer the requirement that the second derivative of the solution exists). A trial solution on each element is then postulated, and this leads to the numerically intensive work of finding the best values for the parameters in the trial solution, and in matching up the various trial solutions from different subdomains.

In general, an FEM solution follows six steps (Li, 2014):

1. Derivation of a weak form of the PDE. This is equivalent to a least-squares minimization of the integral of the difference between the approximate and exact solutions.
2. Discretization of the computational domains.
3. Generation of interpolating or trial functions.
4. Conversion of the “weak form” integral equation into a set of linear equations.
5. Implementation of the boundary conditions.
6. Solution of the resulting linear system of equations.

23.4.1 Weak Form of PDE

Finite-difference methods yield an approximate solution of an approximate PDE. Finite-element methods yield the best possible global agreement between an approximate solution and the exact solution. We start the FEM with the differential equation we want to solve,

(23.4)

We form an integral of the product of the exact solution U(x) and the approximate solution or trial solution Ф(x) over the solution domain. This will be used as a measure of overall agreement between the two solutions. We assume that the trial vanishes at the endpoints, Ф(a) = Ф (b) = 0 (we satisfy general boundary conditions later). We next multiply both sides of the differential equation (23.1) by Ф and integrate by parts from a to b:

(23.5)

(23.6)

(23.7)

(23.8)

Equation 23.8 is the weak form of the PDE, “weak” in the sense that it does not require the existence of the second derivative of U, or the continuity of ρ. Because the approximate and exact solutions are related by the integral of their difference over the entire domain, the algorithm provides the global best agreement between the two.

23.4.2 Galerkin Spectral Decomposition

The approximate solution to the weak PDE proceeds via three steps. First, we split the full domain of the PDE into subdomains called elements, then we find approximate solutions within each element, and finally we match the elemental solutions onto each other. For our 1D problem, the subdomain elements are straight lines of equal length, while for the 2D problem to be considered soon, the elements are triangles (Figure 23.4).

The critical step in the FEM is the expansion of the solution U in terms of a set of basis functions φ_i:

(23.9)

Even when the basis functions are not sines or cosines, this expansion is still called a spectral decomposition. We will choose φ_i’s that are convenient for computation, and so the solution reduces to determining the unknown expansion coefficients α_j. Later, in order to satisfy the boundary conditions, we will add an additional term to this expansion.

Considerable study has gone into determining the effectiveness of different basis functions φ_i that are used to represent the solution on the finite elements. If the sizes of the finite elements are made sufficiently small, then good accuracy is obtained with simple piecewise-continuous basis functions φ_i. For our 1D problem, we use finite elements that are line segments between x_i and x_i+1, and we use basis functions, representing the solution on each line segment, that have the form of triangle or “hats” between x_i−1 and x_i+1 (Figure 23.2). We also require that each basis function equals 1 at the x_i vertex, φ_i(x_i) = 1:

(23.10)

Figure 23.2 Basis functions used in finite-elements solution of the 1D Laplace equation. (a) A set of overlapping basis functions φ_i. Each function is a triangle from x_i−1 to x_i+1. (b) A Piecewise-linear function. (c) A piecewise-quadratic function.

Because this choice means that each basis function equals 0 or 1 at the nodes,

(23.11)

the values of the expansion coefficients α_i must equal the values of the (still unknown) solution at the nodes:

(23.12)

(23.13)

Equation 23.13 makes it clear that the expansion in terms of basis functions is essentially an interpolation between the actual solution at the nodes.

23.4.2.1 Solution via Linear Equations

Because the basis functions φ_i in (23.9) are known, solving for U(x) involves determining the coefficients α_j, which, as we just said, are the unknown values of the true solution U(x) on the nodes. We determine those values by substituting the expansions for U(x) and Ф(x) into the weak form of the PDE (23.8). This converts the integral equation into a set of simultaneous linear equations. As discussed in Chapter 6, there is a standard matrix form for a set of linear equations,

(23.14)

Our equations fit that, with y a vector of unknowns, and where we still need to specify the known stiffness matrix A and the known load matrix b. To that end, we substitute the expansion into the weak form (23.8) to obtain:

By successively selecting we obtain N simultaneous linear equations for the unknown α_j’s:

(23.15)

We factor out the unknown α_j’s and write the equations out explicitly:

Because we have chosen the φ_i’s to be the simple hat functions, the derivatives are easy to evaluate analytically (for other bases they can be carried out numerically):

(23.16)

The integrations are now fairly simple:

(23.17)

(23.18)

(23.19)

We rewrite these equations in the standard matrix form (23.14) with y constructed from the unknown α_j’s, and the tridiagonal stiffness matrix A constructed from the integrals over the derivatives:

(23.20)

(23.21)

The elements in A are just combinations of inverse step sizes and so do not change for different charge densities ρ(x). This is part of what makes FEM so efficient once set up. The elements in b do change for different ρ’s, but the required integrals can be performed analytically or with Gaussian quadrature (Chapter 5). Once A and b are computed, highly efficient methods from a linear algebra package are used to solve the matrix equations for y, and thus the expansion coefficients α_j.

23.4.2.2 Dirichlet Boundary Conditions

Because the basis functions vanish at the endpoints, a solution expanded in them must also vanishes there. This will not do in general, and so we must add to our general solution U(x), a particular solution U_aφ₀(x) that satisfies the boundary conditions (Li, 2014):

(23.22)

where U_a = U(x_a). We substitute U(x) − U_aφ₀(x) into the weak form of the PDE to obtain (N + 1) simultaneous equations, still of the form Ay = b, but now with

(23.23)

This is equivalent to adding a unit element to A and adding a new load vector element:

(23.24)

To impose the boundary condition at x = b, we again add a particular solution U_bφ_N−1(x) and substitute it into the weak form to obtain

(23.25)

So now we need to solve the linear equations Ay = b'. For 1D problems, 100–1000 equations are common, while for 3D problems there may be millions. Because the number of calculations varies approximately as N², it is important to utilize an efficient and accurate algorithm, because otherwise the round-off error can easily dominate the solution.

23.5 1D FEM Implementation and Exercises

In Listing 23.1, we give our program LaplaceFEM_1D.py that determines the 1D FEM solution, and in Figure 23.3 we show that solution. We see on the left of the figure that three elements do not provide even visual agreement with the analytic result, whereas N = 11 elements do.

1. Examine the FEM solution for the choice of parameters
   
   (23.26)
2. Generate your own triangulation by assigning explicit x values at the nodes over the interval [0, 1].
3. Start with N = 3 and solve the equations for N values up to 1000.
4. Examine the stiffness matrix A and ensure that it is triangular.
5. Verify that the integrations used to compute the load vector b are accurate.
6. Verify that the solution of the linear equation Ay = b is correct.
7. Plot the numerical solution for U(x) for N = 10, 100, and 1000, and compare with the analytic solution.

   

   Figure 23.3 Exact (line) vs. FEM solution (points) for the two-plate problem for N = 3 and N = 11 finite elements (N = 3 displaced upwards for clarity). On this scale, the N = 11 solution IS identical to the exact one.

8. The log of the relative global error (number of significant figures) is
   
   (23.27)
   
   Plot the global error vs. x for N = 10, 100, and 1000.

Listing 23.1 LaplaceFEM_1D.py provides an FEM solution of the 1D Laplace equation via a Galerkin spectral decomposition. The resulting matrix equations are solved with Matplotlib.

23.5.1 1D Exploration

1. Modify your program to use piecewise-quadratic functions for interpolation, and compare the results obtained to those obtained with the linear functions.
2. Explore the resulting electric potential and check that the charge distribution between the plates has the explicit x dependence

(23.28)

23.6 Extension to 2D Finite Elements

The steps followed to derive the 2D finite elements method are similar to those for the 1D method, with the big difference being that the finite elements are now 2D triangles as opposed to 1D lines. Figure 23.4 shows how an arbitrarily shaped domain might be decomposed into triangles. Although life is simpler if all the finite elements are of the same size and shape, this is not necessary, and, indeed, as we have shown in the figure, higher precision and faster run times maybe obtained by picking smaller domains in regions where the solution is known to vary rapidly, and picking larger domains in regions of slow variation. As you can imagine, 2D and 3D FEMs can get rather complicated, but not to worry, we will just outline the 2D method and refer the interested reader to Polycarpou (2006) and Reddy (1993) for fuller discussions.

Figure 23.4 (a) Decomposition of a domain into triangular elements. Smaller triangles are used in regions of rapid variation and larger triangles are used in regions of slow variation. Discretization errors occur at boundaries. (b) A decomposition of a rectangular domain into 32 right triangles on a mesh with 25 nodes.

We extend our previous 1D method to solve the 2D version of the Laplace equation (19.4),

(23.29)

There are now 2D Dirichlet boundary conditions:

(23.30)

Here, h is the height and w is the width of the rectangular domain in which we desire a solution. Because our problem domain is now rectangular, it is easy to divide it into right triangles, as we have shown in Figure 23.4b.

23.6.1 Weak Form of PDE

For the 2D problem, the weak form of the PDE again follows from multiplying both sides of the PDE by the trial solution Φ, and then integrating (Polycarpou, 2006):

(23.31)

Here, Ω is a surface boundary of the domain in which we seek a solution, Γ is a perimeter around the surface, U is the solution of the PDE, and n_x and n_y are outward-facing unit normal to Γ. For Dirichlet boundary conditions the contribution of the line integral on the RHS vanishes.

23.6.2 Galerkin’s Spectral Decomposition

As in the 1D method, the approximate solution U(x, y) is expanded in a set φ_i(x, y) of basis functions, in this case 2D functions:

(23.32)

After setting U = φ_j for j = 1, 2, … , N − 1, the weak form of the PDE becomes a set of linear equations:

(23.33)

We rewrite these equations in the standard matrix form (23.14) for linear equations:

(23.34)

Here the U_i’s are the vectors of unknowns, the stiffness matrix elements are

(23.35)

and the load vector b is given by (23.20).

23.6.3 Triangular Elements

Triangular elements are often used in 2D FEM because they can be fit into many arbitrary geometries with little overlap and with little discretization error at the boundary edges (see Figure 23.4). As we see in Figure 23.5a, we take these elements to be triangles of arbitrary shape, with the CCW arrow indicating the direction in which the nodes are numbered. While it is easier to fit arbitrary shaped triangles into a general region, it is easier to do mathematics with right triangles, such as the master triangle shown in Figure 23.5b. The latter have their orthogonal sides lying along the and η axes, with the (x, y) and (, η) coordinates related by a linear coordinate transformation.

These master triangles are the linear interpolation functions in the and η variables that are used in 2D FEM. For example, the linear function for node 1 has the form

(23.36)

The constants are determined by evaluating the functions at each node, for example,

(23.37)

Figure 23.5 (a) Linear triangular elements in the x–y plane. (b) Linear triangular element (master element) in the η plane.

(23.38)

(23.39)

Similar evaluations at the other nodes yield (Polycarpou, 2006):

(23.40)

With these interpolation functions in hand, it is possible to express the x and y coordinates of any point inside an element in terms of the master coordinates:

(23.41)

(23.42)

(23.43)

Next, we take these discrete forms for the interpolation functions, return to the Galerkin spectral decomposition, and use (23.35) to evaluate the A matrix. The required derivatives are evaluated using the chain rule:

(23.44)

(23.45)

We write these equations in the matrix form as

(23.46)

where the 2 × 2 matrix that defines the coordinate transformation between the (x, y) and the (, η) derivatives is called the Jacobian matrix J. After substitution of the explicit forms for the φ’s, the Jacobian takes the simple form:

(23.47)

Likewise, the derivatives in the A matrix can be expressed in terms of the x and y derivatives by using the inverse of the Jacobain matrix:

(23.48)

(23.49)

Continued evaluation of the Galerkin matrix elements yields

(23.50)

(23.51)

After similar evaluations for φ₂ and φ₃, we obtain the six needed derivatives:

(23.52)

(23.53)

(23.54)

23.6.4 Solution as Linear Equations

The final evaluations of the stiffness matrix elements are made using the and η coordinates, for example,

(23.55)

These elements are found to form a symmetric matrix with values:

(23.56)

(23.57)

(23.58)

Next we evaluates the coordinate transformations:

(23.59)

(23.60)

After substituting for and η, we are left with the desired interpolation functions expressed in terms of just x and y.

23.6.5 Imposing Boundary Conditions

The procedure to impose Dirichlet’s boundary conditions for the 2D case is essentially the same as that for the 1D case (Section 23.4.2.2), with it now applied to all nodes that lie on the boundary Γ.

Listing 23.2 The code LaplaceFEM_2D.py solves the 2D Laplace equation using a finite elements method.

23.6.6 FEM 2D Implementation and Exercise

As shown in Figure 23.4b, our application of 2D FEM has the solution domain covered by a collection of triangular elements. Each triangle in the mesh is numbered, in this case from 1 to 32. In addition, the three vertices of each triangle are numbered in a counter-clockwise direction from 1 to 3. Next, each node in the mesh (the dark circles in Figure 23.4 where lines intersect) is numbered, in this case from 1 to 25. Accordingly, the stiffness matrix A in (23.34) has dimension 25 × 25, while the load vector b has dimension 25 × 1.

Listing 23.2 presents our implementation of the 2D FEM solution to the 2D Laplace equation, based on the Matlab code of Polycarpou (2006). It utilizes 800 elements and 441 nodes. The output of this code is essentially the same as our solution to the same problem using the finite-differences method.

23.6.7 FEM 2D Exercises

1. Examine the effect of varying the domain height and width, as well as the number of elements.
2. Compare this numerical solution to the analytic one (the Fourier series in Section 19.3) and determine how the precision changes as the number of elements is varied.
3. Modify the program so that it solves the parallel plate capacitor problem and compare to the finite-difference solution.