problem boundary value ODE solved by finite element method

Consider the two-point boundary value problem -u'' = pi^2 cos(pi x) for x in (0,1), with u(0) = 1 and u(1)=-1.

• Find the analytical solution of this problem.
• Program the finite element method for this problem, using continuous piecewise linear elements.
• Plot the log of the error in the solution at the mesh points versus the log of the number of basis functions, for 2^n elements with n between 1 and 10. What is the slope of this curve (ie, the order of convergence)?
• Plot the log of the error in the derivative of the solution at the mesh points versus the log of the number of basis functions, for 2^n elements with n between 1 and 10. What is the slope of this curve?
• Suppose that we change the right-hand boundary condition to u'(1) = 0. Modify the finite element method to handle this boundary condition. What is the order of convergence of u(1) and u'(1) for this finite element approximation?

Solution:

The independent and dependent variables in the ODE are automatically separated. So the general solution to the ODE is

u(x) = cos (pi x) + C1 x + C2

And using the two Dirichlet boundary conditions, we found C1 = C2 = 0. So the solution to the BVP is

u(x) = cos (pi x)

Code:

1DFEM.c is written in C++ and can be compiled by "g++ -o FEM 1DFEM.c ". Run the executable file "./FEM". The program asks for boundary conditions at x = 1 on screen: "1 for Dirichlet; other for Newmann". The program then asks "n" on screen, for which there are "L = 2^n" elements in total.

The program output on screen in the last three lines

"n, L, totalerror" , where "totalerror" is the total error of all the 1024+1 dyadic mesh points.

"n, L, maxerror" , where "maxerror" is the maximum error of the 2^n+1 mesh points.

"n, L, maxderiverror" , where "maxderiverror" is the maximum derivative error of the 2^n+1 mesh points.

It also output

"A.dat" : stiffness matrix A

"A1.dat","A2dat" and "A3.dat": the three diagonals of A

"f.dat": vector f

"u.dat": coefficient vector u. (Au = f)

"solution.dat": FEM function values of a scanning of x points.

The FEM function values for Dirichlet and Newmann BCs compared with analytical solution are shown

The errors for different L are listed in this table

The log of total error, max error and max derivative error changes with log L like

The total error, max error and max derivative error converges with slope -2.04, -0.2 and -0.98.