Overview of FEM programs
by Mike Gibbons and Arjun Raj
Note: All programs written in Maple V on unix

Implementing the Finite Element Method (as described in Mike Usher and Prof. Strichartz's paper Splines on Fractals) required several parts. First, we had to implement numerical integration on the Sierpinski gasket (referred to as SG from now on). Then, we had to generate the appropriate FEM matrices and vectors (ie, the energy and inner product matrices, as well as the vector of integral(f*basisi).

The very first step in creating the programs was to create an addressing system for the gasket. We chose the following: each triangle on the gasket can be addressed in an obvious way by a string of 0's, 1's, and 2's. As you read through the string, each digit corresponds to a triangle contraction: 0 means "top" triangle, 1 means bottom right triangle, and the 2 is the bottom left triangle. In other parts of the program, we use this system to address triangles. However, here we want to address vertices. So after the string defining the triangle, we then tack on a 01 or a 02 or a 12. This basically specifies that the vertex in question is the midpoint of the 0th and the 1st vertices of the triangle addressed in the first part of the string (or the 0th and the 2nd, or the 1st and 2nd, respectively). Thus, every vertex can be addressed in this way. An exception: the vertices of the SG are arbitrarily labelled 0, 1, and 2 (they cannot be put into this scheme since they aren't midpoints of any vertices). Also, the points in V1-V0 are labeled 01, 02, 12, since there is no triangle address, as they are midpoints of the entire SG. One advantage to this numbering scheme is that all the vertices have a unique label. Also, the length of the address string to a vertex tells you which level of the gasket the vertex is at; the vertex is in Vm if the length of the string is m+1. However, one disadvantage is that it is quite difficult to determine the neighboring points to any particular vertex. We were able to create an algorithm to do such a thing which was not that complicated, however, and this algorithm is in the procedure make_neighbors in SG.mws. It works by basically hard coding the neighbors for the first few levels; then, for all the levels up to m-1, it uses the fact that if you are at the ith level, then the neighbors are all in Vm, so you can take the triangle address for the point in question, and then tack on the rest of the digits by taking the first few digits of the neighbors of a point in V1-V0 (notice that this all depends on the last 2 digits: ie, if the vertex in question ends in a 12, then you need to tack on digits from the 12 vertex). This orientation (meaning the last 2 digits) is necessary because we want to maintain some sort of canonical ordering to the neighbors. We decided that the easiest and most useful way to do this is rotationally. If you consider any vertex, the 4 neighboring vertices all fall on one side of a circle around the vertex. Then, if you were to traverse the circle clockwise, the first vertex you encounter on that side of the circle would be the first neighbor, and the 2nd would be the 2nd neighbor, etc. For instance, the 1st neighbor of 01 in V3 is 1101, the second is 1102, 3rd is 0012 and the 4th is 0001. One reason this is handy is that you can then easily determine the reverse neighbor; ie, if vertex2 is the 1st neighbor of vertex1, then vertex 1 is the 4th neighbor of vertex2. Thus we can assign the neighbors of the neighbors (at least some of them). The reason for that is that when we want to assign neighbors at the deepest level (Vm), there are only 2 neighbors which are still unassigned, and determining which they are is quite easy (just tack on a different ending).

Anyway, now that we had addressing, we created a maple table to store functions in. We addressed the table as follows: if the table was T, then we would have functions stored in as follows: T[[0,1,2,1,1,2],10] would be the value of function number 10 at the vertex defined by 012112. It should be noted that we used the first 6 function numbers for internal purposes: function 0 simply returns the vertex, 1->4 are the 4 neighbors (in order), and 5 and 6 store the x and y coordinates of each point (used for plotting purposes).

To create the spline basis, one basically just creates another gasket table, which we usually denoted as S. Now for the biharmonic splines, there are 2 types of basis vectors, (we call them f and g, as following the splines paper). Thus, we defined each spline basis element to be a list whose first element was the vertex, and the second element was the type. For instance, [[0,1,2,0,1],f]. It should be noted that there are no f-type vectors on the boundaries (set to 0 on the boundaries). There is the matter of using the fact that the normal derivatives at all vertices should be continuous. The question is if we have say C[i] amount of basis element b[i] (which is a g-type vector), then in which direction is the normal derivative positive, and in which direction is it negative? Here, we used the neighbors: if we are considering the first triangle (as determined by the first 2 neighbors), then the normal derivative is positive. If we are considering the 2nd triangle (3rd and 4th neighbors), then the normal derivative is negative. On the boundaries of the SG, we arbitrarily set the normal derivative of a basis vector to be positive. It should also be noted that we needed to create the actual basis functions and store them somewhere in order to integrate against them and the such. These functions are created by calc_biharms and calc_harms, and are usually stored in T. Whenever a function needs to know where these functions have been stored, it usually has a parameter called harm_funcno.

Now, if we are solving -L(u) + qu = f, we need to have the energy matrix and (if q is constant) the inner product matrix. These can be computed theoretically, and this is done in gen_energy_matrix and gen_iprod_matrix. Note that we must always send our spline basis S along for the ride as well as a structure containing Vm_basis, so that we know all about the neighbors. We also can create the vector for f by using gen_f_vector. That function creates the vector by using Simpson's rule for integrating, which we use everywhere, because it works significantly better than the trapezoid rule. If q is not constant, then we have to perform numerical integrations to determine the matrix integral(q*basisi*basisj). This is performed in gen_q_matrix. Notice that here we are talking about biharmonic splines. To perform the same calculations with harmonic splines, the functions are the same, except that you tack an h on the end (gen_q_matrix -> gen_q_matrix_h, gen_energy_matrix->gen_energy_matrix_h). Once the matrices are generated, most of the computations are basically elementary matrix operations like adding matrices and solving systems of linear equations. Then you can use the function make_pretty_picture (make_pretty_picture_h) to create a function out of the spline representation (the vector), and then plot it using the appropriate plotting functions.

One can also use the matrices to solve eigenvector equations simply by solving the generalized eigenvalue equation Eu=lambda*Gu, where E is the energy matrix, and G is the grammian (inner product) matrix. Thus, all one needs to do is tell the computer to find eigenvalues of inv(G)*E. One can also solve the wave and heat equation by exponentiating the appropriate matrix. A computational note: while these calculations are quite simple in theory, they are quite long in practice, because biharmonic splines in level 5 generate 729x729 matrices. While Maple V was useful in that it allowed for quick coding, it was not very powerful as a numerical computation engine. Hence, for most of these kind of calculations, we worked in Matlab by generating our matrices in Maple and then reading them into Matlab, performing the calculations, and then reading the files back into Maple for plotting and the such. While this sounds complicated, it really isn't, and I have included several sample worksheets to document the process. It should be noted that all the matrices calculated are sparse because only a few basis elements will ever overlap with one another. We used this fact when storing the matrices in maple (using the sparse indexing function), and then used that fact when transferring to matlab. Matlab reads sparse matrices by reading 3 column vectors of the same length, the first two of which are the ith and jth component, and the last of which is the value of the matrix at the ith and jth coordinate (all other values are assumed to be zero). The use of sparse matrices significantly reduces memory usage as well as computation time. The program which stores matrices into files (from maple) in save_matrices.mws.

Anyway, I hope that these programs will be helpful to whomever uses them, either for simple modification or extension, or as a guide to performing the same calculations on the pentagasket or other such things.

A listing of main program files in Maple, along with a short description:

SG.mws: This file contains all the basic functionality of generating gaskets, making neighbors, and plotting them on screen. Also, it has functions for making animations of functions based on a bunch of spline representations (biharmonic) for looking at time dependent solutions...

FEM.mws: This file contains the functions required to generate all the matrices for biharmonic splines. It also contains the function to create a function out of a spline projective vector representation. Note that this file calls a lot of functions in splines.mws...

splines.mws: This file contains a lot of the meat of the functionality for generating energy matrices, as well as inner product matrices and q matrices. It also is used for making the f vector (all for biharmonic splines). Many of these programs are quite complicated...

harm_FEM.mws: This file contains everything for doing the FEM with harmonic splines. It will generate all the matrices and vectors, as well as making a picture out of the final spline representation. All the functions are similarly named to those in FEM.mws, but end with _h.

functions.mws: This file contains a bunch of misc. functions which are used in several other programs. It has numerical integration in it, as well as all sorts of methods for generating function.

harmonic.mws: This file contains programs to create harmonic as well as n-harmonic functions. It is used in other programs when creating basis vectors or creating the final pictures.

save_matrices.mws: This file contains a few rudimentary programs to save enery and inner product matrices to text files. I have already included several of the data files for many of the matrices, so chances are these programs will never be used, but here it is anyways.

error_funcs.mws: This file contains a bunch of functions you might use to determine errors of solutions. Included are methods of calculating the error given a theoretical solution, as well as a method of calculating the error when you don't know what the theoretical solution is.

everything.mws: This file contains all the functions from above in one big file, so it's easier to manage. Note that generally speaking, all of the above program files must be run in order for anything to work.

There are also plain text versions of all the above worksheets. They are named the same thing, but instead of ending .mws, they end .txt.

The following are sample worksheets which illustrate the usage of the functions in the above programs. They also come in .txt format...

harm_FEM_sample.mws: Shows how to run the harmonic Finite Element Method...

FEM_sample.mws: Shows how to run the biharmonic Finite Element method...

wave_eqn_sample: Shows how to create animations of the wave eqn using maple. It should be noted that all the calculations are done using matlab, as documented in matlab_sample.txt. This also requires the document wave_time3.dat, which is a solution to the wave eqn with time steps of .01 and a biharmonic spline basis out to V3. This is a solution to the wave eqn with intial values 0 everywhere, and giving a derivative of 20 at the point 01.

The following are matlab .m files which are used for generating wave and heat solutions:

loadIP2.m -> loadIP5.m: These m files return the inner product matrix for the biharmonic splines at level 2 to 5 respectively.

loadEM2.m -> loadEM5.m: same as above, but returns the energy matrix.

Note: both of the above files require EM2_x.dat, IP4_i.dat, and so forth to be in the home directory, otherwise they won't work. Also, the documentation for them is only in loadEM2.m and loadIP2.m, since they all work quite similarly.

generate_heat_time_series.m: solves the heat eqn using finite differences

generate_wave_time_series.m: solves the wave eqn using finite differences

generate_wave_time_series2.m: solves the wave eqn using matrix exponential

HSnorm.m: returns the Hilbert Schmidt norm of a matrix (used in calculating the spectral radius)

spec_radius.m: returns the spectral radius of a matrix using the HS norm

heat_eqn.m: solves the heat eqn given initial values and values of k and h
heat_eqn_t.m: solves the heat eqn given inital values with a timestep and number of steps

wave_eqn_t.m: like heat_eqn_t, except solves the wave eqn.
wave_eqn_t2.m: solves the wave eqn at a given value of t, but uses matrix exponential.

The following are the text files which are used by the various loads in order to load the matrices (biharmonic)...

EM2_i.dat EM2_j.dat EM2_x.dat
IP2_i.dat IP2_j.dat IP2_x.dat
EM3_i.dat EM3_j.dat EM3_x.dat
IP3_i.dat IP3_j.dat IP3_x.dat
EM4_i.dat EM4_j.dat EM4_x.dat
IP4_i.dat IP4_j.dat IP4_x.dat
EM5_i.dat EM5_j.dat EM5_x.dat
IP5_i.dat IP5_j.dat IP5_x.dat

The following file shows how to use the matlab files: matlab sample

Go to top