LINE3: Continue your 1D finite element program. Write "program3", which adds new stuff to program2. Add a function called BASIS, to evaluate any finite element basis function.

            function value = basis ( t, x, n, j )
            T is a position between BCLEFT and BCRIGHT.
            X is the list of nodes you defined earlier.
            N is the number of nodes.
            J is an index between 1 and N (or 0 and N-1).

Depending on the language you use, the details will differ, but the idea will be roughly that you have to handle the first and last values of J separately; you have to determine which two values of X surround the value of T, and then evaluate the appropriate linear function (up, down, or zero).

            value = 0
            if J = 1, then:
              if T < X(2),            value = ( X(2) - T ) / ( X(2) - X(1) );
            if 1 < J < N then:
              if X(J-1) <  T <= X(J), value = ( T - X(J-1) ) / ( X(J) - X(J-1) ),
              if X(J)   <= T < X(J+1) value = (X(J+1) - T ) / ( X(J+1) - X(J) );
            if J = N then:
              if X(N-1) <  T <= X(N), value = ( T - X(N-1) ) / ( X(N) - X(N-1) ).

Verify that your function is working correctly by having program 3 creating a plot of the values of basis function number J=4, using N = 11 points for the mesh, and 101 equally spaced points T in [BCLEFT,BCRIGHT] for the plot.

LINE4: Continue your 1D finite element program. Write "program4", which adds new stuff to program 3. Evaluate the interpolant to sin(2*pi*t). We still have BCLEFT=0, BCRIGHT=1, N=11 and our points X. Define a vector C of length N, and set C(J) to be the value of sin(2*pi*X(J)). Our interpolant is

            s(t) = c(1) * b(1)(t) + c(2) * b(2)(t) + ...+ c(n) * b(n)(t).

where b(1)(t) is the first basis function, evaluated at t, and so on. Compare s(t) and sin(2*pi*t) at 101 equally spaced points in [0,1], with a plot or a table.

LINE5: Continue your 1D finite element program. Make a copy of program4 and call it program5. We are going to add another tool to the program, which can compute the derivative of any basis function. In program5, make a copy of the basis function code, and call it BASISP. It will look almost the same as basis:

            function value = basisp ( t, x, n, j )

but now we want it to evaluate the derivative of a basis function. Because we are only using piecewise linear functions so far, the derivatives will be very simple to compute:

            if J = 1, then:
              if T < X(2),            value = -1 / ( X(2) - X(1) );
            if 1 < J < N then:
              if X(J-1) <  T <= X(J), value = +1 / ( X(J) - X(J-1) ),
              if X(J)   <= T < X(J+1) value = -1 / ( X(J+1) - X(J) );
            if J = N then:
              if X(N-1) <  T <= X(N), value = +1 / ( X(N) - X(N-1) ).

Verify that your function is working correctly by creating a plot of the derivative of basis function number J=4, using N = 11 points for the mesh, and 101 equally spaced points T in [BCLEFT,BCRIGHT] for the plot.

LINE6: Numerical integration for 1D. Make a copy of program5 and call it program 6 The linear system that will define our solution involves integrals. We will approximate those integrals using a 2 point quadrature rule:

            integral(a<=x<=b) f(x) dx = (b-a) * ( w1 * f(x1) + w2 * f(x2) ).

Here,

w1 = w2 = 1/2;
x1 = ((1+s)*a+(1-s)*b)/2;
x2 = ((1-s)*a+(1+s)*b)/2;
s = 1/sqrt(3).

Using this rule, and assuming we use N=11 equally spaced nodes on [0,1] as usual, estimate the following integrals, printing your answer each time:

integral ( 0 <= x <= 1 ) cos(2*pi*x) dx;
integral ( ELEMENT #3 ) cos(2*pi*x) dx;
integral ( ELEMENT #3 ) basis(#4) dx;
integral ( ELEMENT #3 ) basisp(#3) dx.
integral ( ELEMENT #3 ) basis(#3) * basis(#4) dx;

Except for the first integral, we only want to consider the integral over the region that "belongs" to element 3.

LINE7: The Stiffness Matrix. Make a copy of program6 and call it program 7. Set up a matrix A of size N by N, where N is the number of nodes. We want A(i,j) to be the integral over the whole region of the product of the derivatives of basis functions i and j:

            A(i,j) = integral basisp(i,x) * basisp(j,x) dx

It can be tricky to handle the cases when I=1 and I=N. The cleanest way is to do the following:

            For each ELEMENT
              For I = E(1:2,ELEMENT)
                For J = E(1:2,ELEMENT)
                  A(i,j) = A(i,j) 
                    + integral ( x in ELEMENT ) basisp(i,x) * basisp(j,x) dx

In other words, do the integral one element at a time. Most of the entries of A(i,j) will be zero. Most of the other entries will be simple, especially if you use N = 11.

LINE8: The Right Hand Side. Make a new program called program 8, which is a copy of program 7, plus some new stuff. Supposing the right hand side of our differential equation is f(x), we now need to compute a vector called RHS, again using integration. Define a vector RHS of size (N) or, in MATLAB, of size (N,1). We want to compute

            RHS(i) = integral ( 0 <= x <= 1 ) f(x) * basis(i,x) dx

Again, we fill in the values by considering each element:

            For each ELEMENT
              For I = E(1:2,ELEMENT)
                RHS(i) = RHS(i) 
                  + integral ( x in ELEMENT ) f(x) * basisp(i,x) dx

Use F(X) = -exp(x), use N = 11, compute RHS, and print it.

LINE9: Boundary conditions. You now have a matrix A and right hand side vector RHS. However, we need to adjust A and RHS to account for boundary conditions at the first node, X(1), and the last node, X(N). All the entries in row 1 of A should be zeroed out, except that A(1,1) should be set to 1. RHS(1) should be set to the boundary value. For this problem, that value is 0. Similarly, set all the entries in row N of the matrix A to zero, except that A(N,N) is 1, and set RHS(N) to 0. Print A and RHS.

LINE10: Solve the system. Make a new program #10 which is a copy of program #9. We want to solve the linear system

            A*C=RHS

for C, the coefficients of the finite element basis functions. If you are working in MATLAB, this solution step is easy. If you are working in C/C++ or Fortran, you need to locate a linear solver you can use. Print out the 11 entries of your solution vector C.

LINE11: Tabulate the solution. Copy your previous program and call it program 11. We want to display the solution as a table, and as a graph. Given the coefficients C, the finite element estimate of the solution at any point X is

            UHAT(X) = sum ( 1 <= I <= N ) C(I) * BASIS(I,X)

Use 25 equally spaced points to print a table of UHAT(X). Use 101 equally spaced points to create a plot.

LINE12: Compare computed and exact solution at nodes. Copy your previous program and call it program 12. Our finite element program was solving the following problem:

            -u'' = - exp(x), u(0) = u(1) = 0.

The exact solution to this problem is

            u(x) = exp(x) + x - x exp(1) - 1

Use 101 equally spaced points to create a "line" plot of your solution, AND a "dot" plot of the exact solution. Expect the line to go right through the dots. The error at the nodes is easy to compute, because at those special places, we have that uhat(x(i)) = c(i), so:

            error = u(x(i)) - uhat(x(i)) = u(x(i)) - c(i)

Print the error at the N=11 nodes.

LINE13: Compare computed and exact solution as an integral. Copy your previous program and call it program 13. Estimate the following integral:

            L2 = sqrt ( integral ( 0 <= X <= 1 ) ( U(X) - UHAT(X) )^2 dx )

This is a complicated operation. Think about it in small steps, starting with the idea of integrating over a single element.

LINE14: We want to solve a new problem of the form

            - u'' + 2 * u = 6x^2 - 8x + 4, for 0 <= x <= 1

with boundary conditions

            u(0) = 5, u(1) = 4.

In order to do this, we have to modify our finite element system:

            integral [ u'*b' + 2*u*b ] dx = integral [ (6x^2-8x+4)*b ] dx

We have new boundary conditions and a different right hand side, but the big change is that our matrix A now includes something called the "mass matrix". To set up this new matrix, do the following:

            For each ELEMENT
              For I = E(1:2,ELEMENT)
                For J = E(1:2,ELEMENT)
                  A(i,j) = A(i,j) 
                    + integral ( x in ELEMENT )     basisp(i,x) * basisp(j,x) dx
                    + integral ( x in ELEMENT ) 2 * basis(i,x)  * basis(j,x)  dx

Using N=11 nodes on [0,1] as usual, set up and solve this new linear system. Compare your computed answer UHAT to the exact answer U(X)=3x^2-4x+5.

LINE15: Ordinarily, in MATLAB, you solve a linear system A*x=b by writing x=A\b. In other computer languages, you have to write out the instructions to do the solution, or find a prewritten library routine. For the finite element method, one alternative is to use a simple iterative method. The Jacobi iterative method starts with an approximate solution x0 to the linear system, and tries to return an improved estimate x1. The formula is:

            x1(i) = ( b(i) - sum ( j =/= i ) A(i,j) * x0(j) ) / A(i,i)

Write a function like this: function x1 = jacobi ( A, x0, b ) which does one step of the Jacobi iteration. Let A be the 11x11 matrix which looks like this: [ 2 -1 0 0 0 ...0 0 ] [ -1 2 -1 0 0 ...0 0 ] [ 0 -1 2 -1 0 ...0 0 ] ... [ 0 0 0 0 0 ..-1 2 ] let b be the vector [ 0, 0, 0, ..., 0, 11 ] and let the starting vector x0 be [ 1, 1, 1, ..., 1 ]. The norm of the residual r is the maximum value of |A * x - b|. What is the residual for x0? Compute the improved Jacobi estimate x1. What is the residual for x1? Use x1 as the input to jacobi to get x2, and repeat this process 10 times, so that you get x10. Print the residuals for x0 through x10.

LINE16: Write a program which starts with the matrix A, right hand side b, and initial guess x0 that we used in the previous example. But now the program must repeatedly call jacobi until the norm of the residual is less than 0.00001. When you satisfy this tolerance, print out the number of iterations you took, and the estimated solution.

LINE17: Make a copy of PROGRAM14, but now use your Jacobi function to solve the linear system, instead of Matlab's solver. Again, try to get an answer for which the residual is less than 0.00001. How many times did you have to call jacobi?