Step 1: Loading The Code

Try to not overlook this step, please.  The SG code can be found on Gibbons website here, and the SG irregular decomposition can be found on Coletta's here.  Start up your favorite version of Maple (7 or higher please), and open everything.mws, new_stuff.mws, and product_fem.mws.  Respectively, this is the code for the Sierpinski Gasket, irregular decomposition, and my new code for the products.  For each of these worksheets, select Edit -> Execute -> Worksheet from the menu, and viola! the code is loaded.

Step 2: Creating the Gasket

The very first thing to do is to decide how the gasket will be decomposed.  For more information about how our decomposition works, or why we decompose follow this link.  For this tutorial, we want to have the singularity of our function at [[0,1],[0,1]], so we want to create a decomposition centered around this point.  To this effect, we must first generate a wordlist of cells to decompose:

> W:=reg_wordlist(1);
W:=zoom_wordlist(W,[0,1],1);
W:=zoom_wordlist(W,[0,1],1);

reg_wordlist creates an empty wordlist, and zoom_wordlist decomposes the cells that contain  [[0,1],[0,1]].  The next thing to do is to make the wordlist for the SG into a wordlist for SG2.

> W:=wordlist_makeprod(W);

We now have our decomposition, and can create our gasket and the vertex list.

>  T := pd_gasket_irreg(W);
    vm := pd_vertices_irreg(W);

pd_gasket_irreg(W) creates a gasket based on the wordlist W, and pd_vertices_irreg generates the corresponding vertex list.  For more information on what the gasket looks like to Maple, click here.  Next, we generate the listing of the neighbors, which are very useful later on.

>  pd_make_neighbors_irreg(T,W);

pd_make_neighbors needs the gasket and the wordlist as its input. 

Our final step is to generate the pluribiharmonic basis elements that will be used by our spline.

> pd_calc_biharms(T,R,vm,4,20);

This saves the 27 biharmonic basis functions in slots 20 - 47 (argument 5).  These can be put anywhere you desire; just remember that it will take 27 slots!  T is our gasket, and R should always be passed in the second argument (as a throwback to the old Gibbons code).  No code for generating the biharmonic basis elements with an irregular decomposition exists (perhaps you would like to write it?), so we still have to pass the vertices list in argument 3, and argument 4 is the depth to calculate the biharmonic basis to.

Congratulations!  You have successfully generated a SG²!

Step 3:  Create the Spline

The spline is created in several steps.  First, we create the outline of it; it is held in Maple just like the gasket is.  Note that the spline, of course, is not as detailed as the gasket is.  As before, we create the wordlist, then the gasket, then the vertices, etc.

> W2 := reg_wordlist(1);
   W2 := zoom_wordlist(W2,[0,1],1);
   S := pd_gasket_irreg(W2);
   vm2 := pd_vertices_irreg(W2);
   pd_make_neighbors_irreg(S,W2)

And yeah... that's it.  Rather anti-climatic after doing the gasket, eh?

Step 4:  Generating the Energy Matrix and the F vector

This is another easy step.  The energy matrix is generated by a simple call to pd_gen_energy_matrix, passing the spline and the corresponding wordlist.

>  EM := pd_gen_energy_matrix_irreg(S,W);

Now, make the F-vector.  It will be the same dimension as EM, but consist of all zeros except for the one corresponding to the delta mass - the singularity in the Green's Function.  For this example, let's use [[0,1],[0,1]] as the singularity.

>  f_vec :=  Vector(RowDimension(EM),0);
    f_vec[1] := 1;

Step 5:  Solving the problem

Here's where you get to make your own decisions!  If you have gnu-octave correctly installed, and in your PATH, you can get the solution really quick... if not, you will have to wait (forever) on the Maple solver.  To install octave, hop on over to http://www.octave.org and download and install the appropriate software, then make sure that octave is included in your path (In Linux/Unix, type 'octave' in a command prompt, and see if it starts up.  In windows, simply type 'octave.exe.lnk'.  Email me when you get confused).  The Maple code to do this is

>  C:=evalf(linsolve(EM,f_vec),20);

This will be EXTREMELY slow for any decomposition more than level 0.  My recommendation is to use the octave code; first, save the matrix and vector, execute octave, and load the solution:

Windows:

>  save_EM_f();
    system("octave.exe.lnk matrix_solve_octave");
    load_sol();

Linux:

>  save_EM_f();
    system("octave matrix_solve_octave");
    load_sol();

Step 6:  Getting something out of it...

So, now you have a vector, sol, that gives you the solution in terms of the basis functions, which is... not... so useful.  We need to take this and convert it back to point-values on the gasket.  My favorite way to do this is to take the solution, and look at it only on one of the factors;  make_SG_plot will do this for you:

>  make_SG_plot(W,3,sol,[0,1]);

W is the wordlist for the spline, 3 is the detail to compute to, sol is the solution vector, and [0,1] is the point you want to look at (this will plot [0,1] x SG'').  You should at least get a pretty picture out of it.

Step 7:  Take a Coffee Break

Congrats!  You've just solved a fractal differential equation on the product of two Sierpinski Gaskets.  Go take a coffee break and brag to your coworkers.  Then look up Greens_Func in the help files of this site - you just done in 20 or so lines what you can do in one.  Ha!  You should have read more of the documentation before you jumped straight into the tutorial....

>  Greens_Func([[0,1],[0,1]],2,[[1,2],[1,2]],0,1);