 |
>
#This function returns a table R which contains the values of p and q as defined
#by Mike Usher for determining the values of n-harmonic extensions. The values
#a and b are values of inner products. The function basically implements the
#formula for determining extensions based on an iterative method involving inner
#products. See "Splines on Fractals" from more information (ie, don't ask me,
#I just work here...).
ps_n_qs:=proc(n)
local i,l,j,R,m;
R:=table();
R[p,0]:=2/5;
R[q,0]:=1/5;
R[a,0]:=2;
R[b,0]:=-1;
R[p,1]:=-3/25;
R[q,1]:=-7/75;
R[a,1]:=7/45;
R[b,1]:=4/45;
for i from 1 to n-1 do
R[a,i+1]:=(2/15)*(1/(5^(2*i)-6*5^(i-1) + 1/5)) * sum( (((4*5^i-4/15)*R[p,i+1-j]+(5^i+13/15)*R[q,i+1-j])*(R[a,j]+2*R[b,j])), j=1..i);
R[b,i+1]:=(2/15)*(1/(5^(2*i)-6*5^(i-1) + 1/5))*sum( (((3*5^i-13/15) * R[p,i+1-j] + (2*5^i-14/15) * R[q,i+1-j]) * (R[a,j] + 2*R[b,j])),j=1..i);
R[p,i+1]:=(-1/5)*((9/5)*(R[a,i+1] + 2*R[b,i+1]) + sum( ((4*R[p,i+1-m] + R[q,i+1-m])*R[a,m] + (3 * R[p,i+1-m] + 2*R[q,i+1-m]) * R[b,m]),m=1..i));
R[q,i+1]:=(-1/5)*((7/5)*(R[a,i+1] + 2*R[b,i+1]) + sum( ((2*R[p,i+1-m] + 3*R[q,i+1-m])*R[a,m] + (4 * R[p,i+1-m] + R[q,i+1-m]) * R[b,m]),m=1..i));
od;
R;
end;
>
#As long as I'm here, why don't we just calculate this thing out to 5
#so that we can do sex-harmonic functions (Yeah, baby!)
R:=ps_n_qs(5);
>
eval(R):
>
#This function creates a n-harmonic function, using the values for harmonic extensions calculated
#in R. It uses the recursive method of calculating the n-harmonic extension based on Mike
#Usher's work. In order to understand how the recursion addressing works, I recommend you look
#at the documentation of make_points and gen_points in SG.mws, where I fully explain how the
#recursion works (in terms of the addressing). Basically, after that, this function
#is quite straightforward; it simply calculates the next level based on the values of the previous
#level, as well as the necessary laplacians. This procedure is called by the wrapper function
#n_harm.
nharm:=proc(T,R,n,m,add0,add1,add2,funcno)
local i,j;
if m <> 0 then
#print(funcno);
#print(add0,add1,add2);
if add1[nops(add1)-1..nops(add1)] = [0,1] then
#This loop is because we need to calculate up to n-1 laplacians of the function
#at this level as well...
for i from 0 to n-1 do
#the sums here determine the value at the next level based on the values of
#various laplacians and the such before. Again, see Mike Usher's work for a
#derivation of this formula.
T[[op(add1[1..nops(add1)-2]),0,0,1],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add1,funcno+n-1-i+j]) + R[q,j] * T[add2,funcno+n-1-i+j] ),j=0..i);
T[[op(add1[1..nops(add1)-2]),0,1,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add1,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add0,funcno+n-1-i+j] ),j=0..i);
T[[op(add1[1..nops(add1)-2]),0,0,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add1,funcno+n-1-i+j] ),j=0..i);
nharm(T,R,n,m-1,add0,[op(add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2],funcno);
nharm(T,R,n,m-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add1[1..nops(add1)-2]),0,0,1],funcno);
nharm(T,R,n,m-1,add2,[op(add1[1..nops(add1)-2]),0,0,2],[op(add1[1..nops(add1)-2]),0,1,2],funcno);
fi;
if add1[nops(add1)-1..nops(add1)] = [1,2] then
for i from 0 to n-1 do
T[[op(add1[1..nops(add1)-2]),1,1,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add1,funcno+n-1-i+j]) + R[q,j] * T[add2,funcno+n-1-i+j] ),j=0..i);
T[[op(add1[1..nops(add1)-2]),1,0,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add1,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add0,funcno+n-1-i+j] ),j=0..i);
T[[op(add1[1..nops(add1)-2]),1,0,1],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add1,funcno+n-1-i+j] ),j=0..i);
od;
nharm(T,R,n,m-1,add0,[op(add1[1..nops(add1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1],funcno);
nharm(T,R,n,m-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)-2]),1,1,2],funcno);
nharm(T,R,n,m-1,add2,[op(add1[1..nops(add1)-2]),1,0,1],[op(add1[1..nops(add1)-2]),1,0,2],funcno);
fi;
if add1[nops(add1)-1..nops(add1)] = [0,2] then
for i from 0 to n-1 do
T[[op(add1[1..nops(add1)-2]),2,0,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add1,funcno+n-1-i+j]) + R[q,j] * T[add2,funcno+n-1-i+j] ),j=0..i);
T[[op(add1[1..nops(add1)-2]),2,0,1],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add1,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add0,funcno+n-1-i+j] ),j=0..i);
T[[op(add1[1..nops(add1)-2]),2,1,2],funcno+n-1-i] := sum((1/(5^(j*(nops(add1)))))*(R[p,j] * (T[add0,funcno+n-1-i+j] + T[add2,funcno+n-1-i+j]) + R[q,j] * T[add1,funcno+n-1-i+j] ),j=0..i);
od;
nharm(T,R,n,m-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[op(add1[1..nops(add1)-2]),2,1,2],funcno);
nharm(T,R,n,m-1,add1,[op(add1[1..nops(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2],funcno);
nharm(T,R,n,m-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(add1[1..nops(add1)-2]),2,0,1],funcno);
fi;
fi;
end;
>
#So this is basically a wrapper function around the function nharm. This function is
#necessary because nharm won't work until V1, due to idiosyncracies in the labeling system.
#This function basically computes the first n-harmonic extension, and then sends the rest out to
#nharm to compute. n stands for n-harmonic, and R is the table of harmonic extension values
#calculated using ps_and_qs. In order to make this function work, you have to put in the values
#of the function at the boundaries in funcno, as well as the values of the laplacians of the
#function in funcno+1, L^2 of the function in funcno+2... and so forth. This uniquely determines
#an n-harmonic function (ie, you need to specify boundary values of L^(0->n-1)). Notice that
#the algorithm developed required that the function as well as all of it's laplacians be
#defined out to all levels. This has the added bonus that while calculating the function, you have
#calculated the laplacian of it in funcno+1, L^2(f) in funcno+2... Fun, isn't it?
n_harm:=proc(T,R,n,m,add0,add1,add2,funcno)
local i,j;
for i from 0 to n-1 do
T[[0,1],funcno+n-1-i]:=sum((1/(5^(j*nops(add1))))*(R[p,j]*(T[add0,funcno+n-1-i+j]+T[add1,funcno+n-1-i+j])+R[q,j]*T[add2,funcno+n-1-i+j]),j=0..i);
T[[1,2],funcno+n-1-i]:=sum((1/(5^(j*1)))*(R[p,j]*(T[add1,funcno+n-1-i+j]+T[add2,funcno+n-1-i+j])+R[q,j]*T[add0,funcno+n-1-i+j]),j=0..i);
T[[0,2],funcno+n-1-i]:=sum((1/(5^(j*1)))*(R[p,j]*(T[add0,funcno+n-1-i+j]+T[add2,funcno+n-1-i+j])+R[q,j]*T[add1,funcno+n-1-i+j]),j=0..i);
od;
nharm(T,R,n,m-1,[0],[0,1],[0,2],funcno);
nharm(T,R,n,m-1,[1],[1,2],[0,1],funcno);
nharm(T,R,n,m-1,[2],[0,2],[1,2],funcno);
end;
>
#This function recursively generates a harmonic function. It is a much simpler algorithm than
#the n-harmonic algorithm, so the addressing business is much more transparent here. Basically,
#it implements the 2/5, 2/5, 1/5 rule for generating the harmonic extension, and recursively
#generates values all the way down into the gasket. The recursion is a little bit complicated,
#because each step of the way, you need to know the orientation of the triangle in order to know
#what type of points are assigned what values. Basically, we used a rotational method of
#determining what the orientation is: we made it so that the 0th address always was a vertex of the
#previous triangle (ie, the 0th address of the first triangle is [1]...). This way, we always
#which address to assign to, based on the last digits of address 1 (if address 1 is [1,2], then
#we know the triangle is a 1 triangle. Thus, you know what addresses to assign what values.
#Anyway, it works. If you want a somewhat more transparent example, see make_points and gen_points
#in SG.mws. The recursion method used there is exactly the same as used here, and should be easily
#generalizable to other fractals (ie, pentagasket...).
harm:=proc(T,m,add0,add1,add2,funcno)
if m <> 0 then
#print(funcno);
#print(add0,add1,add2);
if add1[nops(add1)-1..nops(add1)] = [0,1] then
T[[op(add1[1..nops(add1)-2]),0,0,1],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add1,funcno] + 1/5 * T[add2,funcno];
T[[op(add1[1..nops(add1)-2]),0,1,2],funcno] := 2/5 * T[add1,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add0,funcno];
T[[op(add1[1..nops(add1)-2]),0,0,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add1,funcno];
harm(T,m-1,add0,[op(add1[1..nops(add1)-2]),0,0,1],[op(add1[1..nops(add1)-2]),0,0,2],funcno);
harm(T,m-1,add1,[op(add1[1..nops(add1)-2]),0,1,2],[op(add1[1..nops(add1)-2]),0,0,1],funcno);
harm(T,m-1,add2,[op(add1[1..nops(add1)-2]),0,0,2],[op(add1[1..nops(add1)-2]),0,1,2],funcno);
fi;
if add1[nops(add1)-1..nops(add1)] = [1,2] then
T[[op(add1[1..nops(add1)-2]),1,1,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add1,funcno] + 1/5 * T[add2,funcno];
T[[op(add1[1..nops(add1)-2]),1,0,2],funcno] := 2/5 * T[add1,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add0,funcno];
T[[op(add1[1..nops(add1)-2]),1,0,1],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add1,funcno];
harm(T,m-1,add0,[op(add1[1..nops(add1)-2]),1,1,2],[op(add1[1..nops(add1)-2]),1,0,1],funcno);
harm(T,m-1,add1,[op(add1[1..nops(add1)-2]),1,0,2],[op(add1[1..nops(add1)-2]),1,1,2],funcno);
harm(T,m-1,add2,[op(add1[1..nops(add1)-2]),1,0,1],[op(add1[1..nops(add1)-2]),1,0,2],funcno);
fi;
if add1[nops(add1)-1..nops(add1)] = [0,2] then
T[[op(add1[1..nops(add1)-2]),2,0,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add1,funcno] + 1/5 * T[add2,funcno];
T[[op(add1[1..nops(add1)-2]),2,0,1],funcno] := 2/5 * T[add1,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add0,funcno];
T[[op(add1[1..nops(add1)-2]),2,1,2],funcno] := 2/5 * T[add0,funcno] + 2/5 * T[add2,funcno] + 1/5 * T[add1,funcno];
harm(T,m-1,add0,[op(add1[1..nops(add1)-2]),2,0,2],[op(add1[1..nops(add1)-2]),2,1,2],funcno);
harm(T,m-1,add1,[op(add1[1..nops(add1)-2]),2,0,1],[op(add1[1..nops(add1)-2]),2,0,2],funcno);
harm(T,m-1,add2,[op(add1[1..nops(add1)-2]),2,1,2],[op(add1[1..nops(add1)-2]),2,0,1],funcno);
fi;
fi;
end;
>
#This function is basically a wrapper for harm
#It generates a harmonic function with val0,val1,val2 on the boundaries.
#This function calculates the first harmonic extension, and then calls harm
#on the 3 parts to generate the rest of the harmonic function. It can't do it directly due
#to idiosyncracies in the nomenclature of the V0 points.
harmonic:=proc(T,m,val0,val1,val2,funcno)
T[[0],funcno]:=val0;
T[[1],funcno]:=val1;
T[[2],funcno]:=val2;
T[[0,1],funcno]:=2/5 * val0 + 2/5 * val1 + 1/5 * val2;
T[[0,2],funcno]:=2/5 * val0 + 2/5 * val2 + 1/5 * val1;
T[[1,2],funcno]:=2/5 * val1 + 2/5 * val2 + 1/5 * val0;
harm(T,m-1,[2],[0,2], [1,2], funcno);
harm(T,m-1,[0],[0,1], [0,2], funcno);
harm(T,m-1,[1],[1,2], [0,1], funcno);
end;
|
