|
>
#So let me step you through the process here.
#First, you will want to read the software library
#as well as all the necessary maple stuff...
read "everything.txt":
with(plots):
with(linalg):
>
#Now let's make the table in which we will store all
#of our functions. The 7 refers to V7.
T:=gasket(7):
>
#Now, let's create a triharmonic function with which we will do all sorts
#of fun things with. Note that the values of the function will be stored in
#function 8, while the laplacian will be in 9, and the laplacian squared will be
#in 10. The first 9 lines are basically setting the boundary values which the
#algorithm uses to compute the rest of the function.
T[[0],8]:=0;
T[[0],9]:=50;
T[[0],10]:=-100;
T[[1],8]:=0;
T[[1],9]:=10;
T[[1],10]:=40;
T[[2],8]:=0;
T[[2],9]:=-70;
T[[2],10]:=-130;
#Here is the call to n_harm which calculates the rest of the function.
#T is where all the functions are stored, R is a table for calculating
#n-harmonic functions which is already stored when you run everything.txt.
#3 means tri-harmonic, 7 means calculate out to V7. The [0],[1],[2] should
#always remain that way, and the 8 is where the function is stored. Notice that
#after the function is run, the laplacian will be in function number 9, and L^2
#will be in function number 10.
n_harm(T,R,3,7,[0],[1],[2],8);
>
#Now suppose u was function number 8. Then f=-Lu+qu. Suppose q was the
#identically constant function 1. Then f=-Lu+u, and we know the exact value
#of Lu, because it is stored in function #9. Thus, let's calculate f and store
#it in function #11...
add_func_lin(T,7,8,9,1,-1,11);
#Here, T is the table where we store all the functions, 7 means we calculate out
#to V7, 8 is u, 9 is Lu, and the 1 and -1 basically implement 1*(fn8)+ (-1)*(fn9)
#The final result is stored in 11.
>
#Well, now that we know f, we can go ahead and start up the FEM (biharmonic).
#Step 1 is to calculate the biharmonic basis functions...
calc_biharms(T,R,7,20);
#Here, T means we will store it in table T. R is the table used to create n-
#harmonic functions, 7 means calculate to V7, and 20 means put the results in
#function number 20. It should be noted that calc_biharms uses a lot of
#functions, and the next safe place in T in which to assign functions is function
#number 50.
>
#Let's just set our spline basis out to V2...
m_basis:=2;
>
>
#Now, let's define our spline basis S. It is defined just like T, except
#that we will need the vertices of the basis (vm), and will need to find
#the neighbors of the vertices of S at level m_basis. Hence, the following...
S:=gasket(m_basis):
vm:=vertices(m_basis):
make_neighbors(S,m_basis,vm):
#IT IS VERY IMPORTANT TO RUN MAKE_NEIGHBORS HERE...
>
#The next step is to define the energy and inner product matrices. That
#is accomplished as follows...
EM:=gen_energy_matrix(S,m_basis,vm):
IP:=gen_iprod_matrix(S,m_basis,vm):
#Notice that I have to send the basis S and vm along every step of the way...
>
#So the matrix equation is now (EM+IP)x=y, since q=1.
#Let's add the matrices EM and IP together and store them
#in EM_IP.
EM_IP:=evalm(EM + IP):
>
#So what is y? Well, that can be generated using gen_f_vector...
f_v:=gen_f_vector(T,S,m_basis,7,vm,11,20):
#Here, the 7 means that f is defined out to V7. The 11 is the function in
#which f is stored. The 20 stands for the function number we used in
#calc_biharms, where all the basis elements are stored.
>
#Well, all that's left to do is solve the matrix equation.
#fortunately, we can just let Maple do all the work. The call to evalf
#just means to use floating point numbers; otherwise, the computations can
#take FOREVER...
C:=evalf(linsolve(EM_IP,f_v),20);
#ie, the final vector is stored in C; the 20 just means 20 decimal points
#worth of precision...
>
#So now we have the spline representation of our function stored in C.
#I guess that's cool, but it's nice to see it as well; hence the function
#call to make_pretty picture. The following call creates a function in
#function slot number 15, using the biharmonic basis elements in 20, and
#using the spline representation in C. The 7 refers to calculating the
#values of the function out to V7...
make_pretty_picture(T,S,m_basis,7,vm,15,20,C):
>
#The only problem thus far is that we haven't yet defined the x and y coordin-
#ates for the vertices of the gasket, which we will need if we want to display our
#function. Fortunately, we have a function known as make_points, which will do
#just that. It stores the x and y coordinates of the vertices in function slots
#5 and 6 respectively. The 7 means calculate coordinates out to V7...
make_points(T,7);
>
#Now we have to extract the actual values from the table T in such a format
#that maple can understand it and plot it. Hence SG_plot, which pulls
#function 15 out of T to level V7 and creates the list final_points which
#can be used by maple to plot...
final_points:=SG_plot(T,7,15):
>
#Let's store the plot of our solution...
p15:=pointplot3d(final_points, color=blue, symbol=POINT,title=`15(APPROX) is blue, 8(THEORY) is red`):
>
#Remember that fn #8 was our theoretical solution...
final_points2:=SG_plot(T,7,8):
>
p8:=pointplot3d(final_points2, color=red, symbol=POINT):
>
#And now let's display the final answers on top of each other...
display(p15,p8,axes=NORMAL);
#Came out pretty good, huh?
>
#Now let's determine the errors in our solutions...
>
#The following function takes the difference between 15 (our solution)
#and 8 (theoretical solution) and stores the difference in function
#number 55. It also uses up function 56 to store the difference squared...
error_fns(T,7,15,8,55);
#Note: the 7 in this function call, as well as all subsequent calls,
#refers to calculating out to V7.
>
#This function takes the L2 norm of the function in 55. However, it really
#uses the fact that function 56 contains the square of the difference...
L2_norm(T,7,55);
>
#This function returns the square root of the energy of the function stored
#in function 55. Notice that we have to run make_neighbors on T at level 7
#before we can calculate the energy (in this case, the energy error)...
vm7:=vertices(7):
make_neighbors(T,7,vm7);
get_root_energy(T,7,55);
>
#This function finds the maximum of a function (in this case, the maximum
#error).
find_max(T,7,55);
>
#We have used q = 1 throughout. Unfortunately, the posteriori error function
#only works with q defined as some function. So first, we have to define a
#function which is identically 1 and store it in function #14...
vm7:=vertices(7):
for i from 1 to nops(vm7) do
T[vm7[i],14]:=1 od:
q_funcno:=14; #Just to make it clear what goes where...
harm_funcno:=20;
p_funcno:=15; #as in the approximate (projective) solution
f_funcno:=11;
l_funcno:=58; #This is the slot in which we will store the laplacian of p
m_func:=7; #This basically means that all functions are defined to V7...
>
>
# This function determines the posteriori error (ie, -Lp+qp-f) where p is the approximate
#solution we just calculated.
posteriori_error_q_nonc(T,S,m_basis,m_func,C,p_funcno,q_funcno,f_funcno,harm_funcno,l_funcno);
>
#And that's just about it. The only thing which might change is if you use q nonconstant.
#How might you deal with that, you might ask? Well, let's try it...
>
#First, let's make q a harmonic function stored in function slot 14...
harmonic(T,7,1,0,0,14);
#This stores a harmonic function with boundary values 1, 0, and 0 into function 14 (defined
#out to V7).
>
#Now, instead of using the inner product matrix, we must generate the matrix of
#integral(q*basisi*basisj). To do this, we use gen_q_matrix...
QM:=gen_q_matrix(T,S,m_basis,m_func,vm,14,20):
>
>
#Now instead of EM+IP, we need EM+QM...
EM_QM:=evalm(EM + QM):
>
#Let's continue to suppose that our theoretical solution was 8. In order to determine
#f, we need to take -#9 + #14*#8...
mult(T,7,14,8,16); #multiply 14 by 8 out to level V7 and store in 16...
add_func_lin(T,7,9,16,-1,1,11); #store f in function 11, just as before...
>
>
#let's recalculate the f vector...
f_v:=gen_f_vector(T,S,m_basis,7,vm,11,20):
>
#Let's solve the equation again, shall we?
C:=evalf(linsolve(EM_QM,f_v),20);
>
#Let's look at the answer...
make_pretty_picture(T,S,m_basis,7,vm,15,20,C):
>
#Getting ready to plot...
final_points:=SG_plot(T,7,15):
>
#Let's store the plot of our solution...
p15:=pointplot3d(final_points, color=blue, symbol=POINT,title=`15(APPROX) is blue, 8(THEORY) is red`):
>
#And now let's display the final answer.
display(p15,p8,axes=NORMAL);
#Came out pretty good, huh?
|
