The Angiogenesis Project

For this project, your mentors are Professor Howard Levine and Professor Janet Peterson.

Your initial task is to learn about the computer solution of a set of partial differential equations using the finite element method.

To do this, start by copying the FEM1D MATLAB routines. There are about 10 files, which together, constitute a program to solve a simple problem.

To begin with:

print out the files.
figure out how to run the program.
The program solves a problem over a fixed one dimensional interval, and computes an approximate solution U^h. Figure out how to make a plot of the solution.
The program allows you to choose the number of elements, and can report two measures of error between the approximate solution U^h and the exact solution U. If we let h be the size of an element, then we expect that, as h goes to zero, the error will go to zero. We are interested in how fast the error goes to zero; compute the error for a series of values of h and try to determine the rate of convergence. (I can explain what you are looking for here).

See if the code will run faster using a tridiagonal solver. There are several things we can do to improve the handling of the linear system, but to start with, we'll try replacing the line

    x = A\f

by the line

    x = tri_sol ( sub, diag, sup, f )

You'll need the MATLAB DIAG command to make SUB, DIAG and SUP from A, and also the following routines:

TRI_SOL.M, solves a tridiagonal system.
TRI_FAC.M, computes LU factors of tridiagonal system.
UBI_SOL.M, solves the U part of the system.
LBI_SOL.M, solves the L part of the system.
TRI_MUL.M, computes A*X for a tridiagonal system, useful for checking results of TRI_SOL.

Let us try the code on a different problem. Again, we are "cheating", because we are going to start out knowing the exact answer. The function we are interested in is

    u(x) = cos(2x)

on the interval [0,pi]. The corresponding problem to be solved is

    -u''(x) = 4*cos(2x)
    u(0) = 1
    u(pi) = 1

To get the program to solve this problem correctly, you have to:

Modify the computation of the node locations;
Change the computation of the coefficients associated with the boundary values;
Modify the routines RHSFUN and EXACT;

Change the code to solve this problem, and verify that the error is "small", and gets smaller as you increase the number of nodes.

We have seen that the accuracy of the solution improves if we increase the number of elements. But in some cases, it can be just as important to use the elements "wisely", that is, to use many elements in areas where the function is changing rapidly, and not so many where it is linear. The hard part is deciding where the nodes go. After that, there are no other changes to the code! For this task, we are going to have an exact solution:

    u(x) = 1 / ( x^2 + 0.1 )

on the interval [-5,5]. Work out the corresponding right hand side:

    -u''(x) = ??
    u(-5) = 1 / 25.1
    u(5) = 1 / 25.1

First, try using a uniform mesh of nodes, and a sequence of values of N like 11, 21, 41, 81, 161. Record the values of the error for each value of N. The errors should still behave like h^2, but they should be larger than we saw for our first problem. Graph the exact solution, and notice that the solution changes rapidly between -1 and 1. We're going to guess that we should put most of our nodes here. For a given value of N, try putting the nodes as follows:

    xx = linspace ( N, -5, 5 )
    x(1) = -5
    x(N) = +5
    for i = 2 : N - 1
      x(i) = (10 / pi) * atan( xx )
    end

Plot XX versus X, and verify that the X values lie between -5 and 5, rise monotonically, and that there are more points near 0. Once you get your variable mesh set up, rerun the program with the same set of values of N, and record the error. For each value of N, is the error smaller than for the uniform mesh?

First, we'd like to see if we can't speed up the FE_ERROR calculation. My first stab at this is to make the FE_LINEAR routine more efficient. Try copying FE_LINEAR2.M and using it in FE_ERROR.

Secondly, we'd like to look at boundary conditions that involve derivatives. Based on what we discuss, try to set up and solve the following problem:

    u(x) = cos(2x)

on the interval [0,pi]. The corresponding problem to be solved is

    -u''(x) = 4*cos(2x)
    u'(0) = 0
    u(pi) = 1

We've already solved a problem with this exact solution. The only change is in our specification of the left boundary condition. It's easy to verify that a boundary condition of the form u(a)=b is correctly computed, but how can we check u'(a)=b?

Thirdly, we're going to look at making the equation we solve a little more complicated. Specifically, let's try to solve a problem of the form

    -u'' + u = f

Our specific example will have be

    -u''(x) + u = x^3 - 9x - 2
    u(-2) = -4
    u(3) = +16

whose exact solution is

    u(x) = x^3 -3x -2

It's time to start sketching your report. You can download a suggested outline for your report.

Last modified on 08 August 2001.