|
|
LOTS OF PROGRAMS - Outline
|
 |
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
|
|
|
