PHYCS 6730 Lab Exercise: Introduction to Finite Element Methods

These exercises experiment with a variety of approximate methods for solving a differential equation. Most of the exercise can be done by hand, but you may use Maple if you feel the need.

We will be working with the elementary ODE 

   R(x) = y'' - 6x = 0  on [0,1] with y(0) = 0 and y(1) = 1.
(See DeVries, First Course in Computational Physics.)

Exercise 1. The exact solution

For the exact solution, the residual function R(x) vanishes on the entire interval. Solve this equation by hand. (Note: this is trivial.)

Exercise 2. Approximating the solution

The finite element method approximates the solution as a linear combination of basis functions. Typically these basis functions have compact support - i.e. they are nonzero only over a small domain. To make it possible to do the calculations by hand, we use a polynomial basis here, but the methods described can be used with any basis set.

Approximate y(x) in the ODE with p(x) where 

  p(x) = a x2 + b x + c.
Here the basis functions are the monomials Phi2(x) = x2, Phi1(x) = x1, and Phi0(x) = x0 = 1. The coefficients a, b, c are to be determined from the boundary conditions and the requirement that the equation is approximately solved. Impose the boundary condition and show that 
  R(x) = 2a - 6x
So we have only one free parameter left to play with. This residual parameterizes the error in our approximation. We want to make it as close to zero on the inverval [0,1] as we can.

Exercise 3. Collocation method

This method picks a set of points on the interval and sets the residual to zero. That is, choose a set {xi} and require R(xi) = 0. Given the number of free parameters, how many points can we use in our case?

Choose the middle of the interval (x1 = 0.5). Show that you get a = 3/2. Plot the approximation and the exact solution side-by-side and discuss the difference.

Exercise 4. Subdomain method

This method divides the region into subdomains and sets the average of the residual on each domain to zero. The average is found by integrating the residual over the subdomain and dividing by the length (unnecessary, of course, if we are making it zero). Given the number of remaining free parameters, how many such subdomains can we use?

Choose the domain [0,1] and set the average to zero. Determine a.

Exercise 5. Least squares method

Minimize the square norm of the residual R: 
  I = int(R2(x),x=0..1)
Show that you also get a = 3/2. But please understand that there are no guarantees that the various methods will always agree.

Exercise 6. Galerkin method

We require that the inner product of the residual and a set of basis functions vanishes. The inner product of two functions is defined as the integral of the product over the interval. That is, we require 

  int( R(x)*Phii(x),x=0..1) = 0
for some choice of values for i. Given the number of free parameters we have, how many such conditions can we impose?

Try Phi1(x) and show that you get a = 2.

Exercise 7. Rayleigh Ritz method

Many differential equations in science and engineering applications can be derived from a variational principle. Minimizing a variational integral gives a function that satisfies the equation. For our simple ODE the variational integral 
  I(y(x)) = int(1/2 (dy/dx)2 + 6 x y, x = 0..1)
is to be minimized for all functions y(x) that satisfy the boundary conditions. It can be shown that the function that minimizes it satisfies the differential equation.

For optimizing our approximation, we plug in the approximation and find the parameters that minimize the variational integral.

Show that for our approximation, 

  I = a2/6 - a/2 + 5/2
and find the value a that minimizes it. Then calculate the variational integral for the exact solution and show that it lowers the minimum value of I.

Discussion

The example here is so simple, we were left with only one parameter to play with, and with each method we had only one linear constraint equation to solve. For more complex problems in any number of dimensions we would use typical finite element basis sets of compact support. Depending on how well we wanted to approximate the solution and the power of our computer, the number of parameters unconstrained by boundary conditions could be extremely large. The constraint equations would generate a large linear system in the parameters. But with a basis of compact support, the linear system would be sparse and easily solved with iterative methods.