Construction of the coarse scale representation of a diffusion operator with coefficients defined on the next finer scale
in two dimension for the Legendre Multiwavelet Basis.

restart:with(linalg):

LePath:="./":
LeName:=Legendre:
Datafile:=cat(LePath,"Data",LeName,".m"):
CoefFile:=cat(LePath,"Coeffs",LeName,".m"):
SaveFile:=cat(LePath,"Operator",LeName,".m"):

Remark concerning the coefficients.

If the computation concerns a full or diagonal coefficient matrix, the number of unknown coefficients should be increased accordingly, as it is shown below
The scripts here are presented for an isotropic diffusion matrix, so that the final computation (the interscale reduction in the next script) can be performed in

a reasonable amount of time.

NTMX:=matrix(4,4):NTMY:=matrix(4,4):NTMZ:=matrix(4,4):
NTMX:=simplify(NTM,{a1=a1x,a2=a2x,a3=a3x,a4=a4x}):
NTMY:=simplify(NTM,{a1=a1y,a2=a2y,a3=a3y,a4=a4y}):
NTMZ:=simplify(NTM,{a1=b1,a2=b2,a3=b3,a4=b4}):

NTMX:=matrix(16,16):NTMY:=matrix(16,16):NTMZ:=matrix(16,16):
NTMX:=NTM:
NTMY:=NTM:
NTMZ:=0:

Definition of the (forward and backward) difference operators, from the one dimensional operator obtained in the first script.

The simplification hypothesis HS and HS2 are the same as that described in the paper: Di corresponds to a centered difference whereas Delta is the centered second order operator.
HS:={DeltaP=Di/2+Delta/2,DeltaM=Di/2-Delta/2}:
HS2:={Di*Di=Delta*4+Delta*Delta}:

The differentiation operator is DF; MDF is its one dimensional equivalent.

DF:=matrix(4,4):
for i from 1 to 4 do for j from 1 to 4 do
DF[i,j]:=A||(-1)[i,j]*(DeltaP+1)+A||1[i,j]*(-DeltaM+1)+A0[i,j];
DF[i,j]:=simplify(DF[i,j],HS);
DF[i,j]:=simplify(DF[i,j],HS2);
od od;
MDF:=matrix(2,2):
for i from 1 to 2 do for j from 1 to 2 do
MDF[i,j]:=R||(-1)[i,j]*(DeltaP+1)+R||1[i,j]*(-DeltaM+1)+R0[i,j];
MDF[i,j]:=simplify(MDF[i,j],HS);
MDF[i,j]:=simplify(MDF[i,j],HS2);
od od;
print(DF,MDF);

The 2 Dimensional differentiation operators are constructed by tensor product.

DF2DX:=matrix(16,16):DF2DY:=matrix(16,16):
DFX:=simplify(DF,{Di=Dx,Delta=Dxx}):
DFY:=simplify(DF,{Di=Dy,Delta=Dyy}):

X: Columns 1 and 2
DF2DX[1,1] :=DFX[1,1]:DF2DX[1,2] :=0 :
DF2DX[2,1] :=0 :DF2DX[2,2] :=DFX[1,1]:
DF2DX[3,1] :=DFX[2,1]:DF2DX[3,2] :=0 :
DF2DX[4,1] :=0 :DF2DX[4,2] :=DFX[2,1]:
DF2DX[5,1] :=0 :DF2DX[5,2] :=0 :
DF2DX[6,1] :=0 :DF2DX[6,2] :=0 :
DF2DX[7,1] :=0 :DF2DX[7,2] :=0 :
DF2DX[8,1] :=0 :DF2DX[8,2] :=0 :
DF2DX[9,1] :=DFX[3,1]:DF2DX[9,2] :=0 :
DF2DX[10,1] :=0 :DF2DX[10,2] :=DFX[3,1]:
DF2DX[11,1] :=DFX[4,1]:DF2DX[11,2] :=0 :
DF2DX[12,1] :=0 :DF2DX[12,2] :=DFX[4,1]:
DF2DX[13,1] :=0 :DF2DX[13,2] :=0 :
DF2DX[14,1] :=0 :DF2DX[14,2] :=0 :
DF2DX[15,1] :=0 :DF2DX[15,2] :=0 :
DF2DX[16,1] :=0 :DF2DX[16,2] :=0 :

X: Columns 3 and 4
DF2DX[1,3] :=DFX[1,2]:DF2DX[1,4] :=0 :
DF2DX[2,3] :=0 :DF2DX[2,4] :=DFX[1,2]:
DF2DX[3,3] :=DFX[2,2]:DF2DX[3,4] :=0 :
DF2DX[4,3] :=0 :DF2DX[4,4] :=DFX[2,2]:
DF2DX[5,3] :=0 :DF2DX[5,4] :=0 :
DF2DX[6,3] :=0 :DF2DX[6,4] :=0 :
DF2DX[7,3] :=0 :DF2DX[7,4] :=0 :
DF2DX[8,3] :=0 :DF2DX[8,4] :=0 :
DF2DX[9,3] :=DFX[3,2]:DF2DX[9,4] :=0 :
DF2DX[10,3] :=0 :DF2DX[10,4] :=DFX[3,2]:
DF2DX[11,3] :=DFX[4,2]:DF2DX[11,4] :=0 :
DF2DX[12,3] :=0 :DF2DX[12,4] :=DFX[4,2]:
DF2DX[13,3] :=0 :DF2DX[13,4] :=0 :
DF2DX[14,3] :=0 :DF2DX[14,4] :=0 :
DF2DX[15,3] :=0 :DF2DX[15,4] :=0 :
DF2DX[16,3] :=0 :DF2DX[16,4] :=0 :

X: Columns 5 and 6
DF2DX[1,5] :=0 :DF2DX[1,6] :=0 :
DF2DX[2,5] :=0 :DF2DX[2,6] :=0 :
DF2DX[3,5] :=0 :DF2DX[3,6] :=0 :
DF2DX[4,5] :=0 :DF2DX[4,6] :=0 :
DF2DX[5,5] :=DFX[1,1]:DF2DX[5,6] :=0 :
DF2DX[6,5] :=0 :DF2DX[6,6] :=DFX[1,1]:
DF2DX[7,5] :=DFX[2,1]:DF2DX[7,6] :=0 :
DF2DX[8,5] :=0 :DF2DX[8,6] :=DFX[2,1]:
DF2DX[9,5] :=0 :DF2DX[9,6] :=0 :
DF2DX[10,5] :=0 :DF2DX[10,6] :=0 :
DF2DX[11,5] :=0 :DF2DX[11,6] :=0 :
DF2DX[12,5] :=0 :DF2DX[12,6] :=0 :
DF2DX[13,5] :=DFX[3,1]:DF2DX[13,6] :=0 :
DF2DX[14,5] :=0 :DF2DX[14,6] :=DFX[3,1]:
DF2DX[15,5] :=DFX[4,1]:DF2DX[15,6] :=0 :
DF2DX[16,5] :=0 :DF2DX[16,6] :=DFX[4,1]:

X: Columns 7 and 8
DF2DX[1,7] :=0 :DF2DX[1,8] :=0 :
DF2DX[2,7] :=0 :DF2DX[2,8] :=0 :
DF2DX[3,7] :=0 :DF2DX[3,8] :=0 :
DF2DX[4,7] :=0 :DF2DX[4,8] :=0 :
DF2DX[5,7] :=DFX[1,2]:DF2DX[5,8] :=0 :
DF2DX[6,7] :=0 :DF2DX[6,8] :=DFX[1,2]:
DF2DX[7,7] :=DFX[2,2]:DF2DX[7,8] :=0 :
DF2DX[8,7] :=0 :DF2DX[8,8] :=DFX[2,2]:
DF2DX[9,7] :=0 :DF2DX[9,8] :=0 :
DF2DX[10,7] :=0 :DF2DX[10,8] :=0 :
DF2DX[11,7] :=0 :DF2DX[11,8] :=0 :
DF2DX[12,7] :=0 :DF2DX[12,8] :=0 :
DF2DX[13,7] :=DFX[3,2]:DF2DX[13,8] :=0 :
DF2DX[14,7] :=0 :DF2DX[14,8] :=DFX[3,2]:
DF2DX[15,7] :=DFX[4,2]:DF2DX[15,8] :=0 :
DF2DX[16,7] :=0 :DF2DX[16,8] :=DFX[4,2]:

X: Columns 9 and 10
DF2DX[1,9] :=DFX[1,3]:DF2DX[1,10] :=0 :
DF2DX[2,9] :=0 :DF2DX[2,10] :=DFX[1,3]:
DF2DX[3,9] :=DFX[2,3]:DF2DX[3,10] :=0 :
DF2DX[4,9] :=0 :DF2DX[4,10] :=DFX[2,3]:
DF2DX[5,9] :=0 :DF2DX[5,10] :=0 :
DF2DX[6,9] :=0 :DF2DX[6,10] :=0 :
DF2DX[7,9] :=0 :DF2DX[7,10] :=0 :
DF2DX[8,9] :=0 :DF2DX[8,10] :=0 :
DF2DX[9,9] :=DFX[3,3]:DF2DX[9,10] :=0 :
DF2DX[10,9] :=0 :DF2DX[10,10]:=DFX[3,3]:
DF2DX[11,9] :=DFX[4,3]:DF2DX[11,10]:=0 :
DF2DX[12,9] :=0 :DF2DX[12,10]:=DFX[4,3]:
DF2DX[13,9] :=0 :DF2DX[13,10]:=0 :
DF2DX[14,9] :=0 :DF2DX[14,10]:=0 :
DF2DX[15,9] :=0 :DF2DX[15,10]:=0 :
DF2DX[16,9] :=0 :DF2DX[16,10]:=0 :

X: Columns 11 and 12
DF2DX[1,11] :=DFX[1,4]:DF2DX[1,12] :=0 :
DF2DX[2,11] :=0 :DF2DX[2,12] :=DFX[1,4]:
DF2DX[3,11] :=DFX[2,4]:DF2DX[3,12] :=0 :
DF2DX[4,11] :=0 :DF2DX[4,12] :=DFX[2,4]:
DF2DX[5,11] :=0 :DF2DX[5,12] :=0 :
DF2DX[6,11] :=0 :DF2DX[6,12] :=0 :
DF2DX[7,11] :=0 :DF2DX[7,12] :=0 :
DF2DX[8,11] :=0 :DF2DX[8,12] :=0 :
DF2DX[9,11] :=DFX[3,4]:DF2DX[9,12] :=0 :
DF2DX[10,11]:=0 :DF2DX[10,12]:=DFX[3,4]:
DF2DX[11,11]:=DFX[4,4]:DF2DX[11,12]:=0 :
DF2DX[12,11]:=0 :DF2DX[12,12]:=DFX[4,4]:
DF2DX[13,11]:=0 :DF2DX[13,12]:=0 :
DF2DX[14,11]:=0 :DF2DX[14,12]:=0 :
DF2DX[15,11]:=0 :DF2DX[15,12]:=0 :
DF2DX[16,11]:=0 :DF2DX[16,12]:=0 :

X: Columns 13 and 14
DF2DX[1,13] :=0 :DF2DX[1,14] :=0 :
DF2DX[2,13] :=0 :DF2DX[2,14] :=0 :
DF2DX[3,13] :=0 :DF2DX[3,14] :=0 :
DF2DX[4,13] :=0 :DF2DX[4,14] :=0 :
DF2DX[5,13] :=DFX[1,3]:DF2DX[5,14] :=0 :
DF2DX[6,13] :=0 :DF2DX[6,14] :=DFX[1,3]:
DF2DX[7,13] :=DFX[2,3]:DF2DX[7,14] :=0 :
DF2DX[8,13] :=0 :DF2DX[8,14] :=DFX[2,3]:
DF2DX[9,13] :=0 :DF2DX[9,14] :=0 :
DF2DX[10,13]:=0 :DF2DX[10,14]:=0 :
DF2DX[11,13]:=0 :DF2DX[11,14]:=0 :
DF2DX[12,13]:=0 :DF2DX[12,14]:=0 :
DF2DX[13,13]:=DFX[3,3]:DF2DX[13,14]:=0 :
DF2DX[14,13]:=0 :DF2DX[14,14]:=DFX[3,3]:
DF2DX[15,13]:=DFX[4,3]:DF2DX[15,14]:=0 :
DF2DX[16,13]:=0 :DF2DX[16,14]:=DFX[4,3]:

X: Columns 15 and 16
DF2DX[1,15] :=0 :DF2DX[1,16] :=0 :
DF2DX[2,15] :=0 :DF2DX[2,16] :=0 :
DF2DX[3,15] :=0 :DF2DX[3,16] :=0 :
DF2DX[4,15] :=0 :DF2DX[4,16] :=0 :
DF2DX[5,15] :=DFX[1,4]:DF2DX[5,16] :=0 :
DF2DX[6,15] :=0 :DF2DX[6,16] :=DFX[1,4]:
DF2DX[7,15] :=DFX[2,4]:DF2DX[7,16] :=0 :
DF2DX[8,15] :=0 :DF2DX[8,16] :=DFX[2,4]:
DF2DX[9,15] :=0 :DF2DX[9,16] :=0 :
DF2DX[10,15]:=0 :DF2DX[10,16]:=0 :
DF2DX[11,15]:=0 :DF2DX[11,16]:=0 :
DF2DX[12,15]:=0 :DF2DX[12,16]:=0 :
DF2DX[13,15]:=DFX[3,4]:DF2DX[13,16]:=0 :
DF2DX[14,15]:=0 :DF2DX[14,16]:=DFX[3,4]:
DF2DX[15,15]:=DFX[4,4]:DF2DX[15,16]:=0 :
DF2DX[16,15]:=0 :DF2DX[16,16]:=DFX[4,4]:

Y: Columns 1 and 3
DF2DY[1,1] :=DFY[1,1]:DF2DY[1,3] :=0 :
DF2DY[2,1] :=DFY[2,1]:DF2DY[2,3] :=0 :
DF2DY[3,1] :=0 :DF2DY[3,3] :=DFY[1,1]:
DF2DY[4,1] :=0 :DF2DY[4,3] :=DFY[2,1]:
DF2DY[5,1] :=DFY[3,1]:DF2DY[5,3] :=0 :
DF2DY[6,1] :=DFY[4,1]:DF2DY[6,3] :=0 :
DF2DY[7,1] :=0 :DF2DY[7,3] :=DFY[3,1]:
DF2DY[8,1] :=0 :DF2DY[8,3] :=DFY[4,1]:
DF2DY[9,1] :=0 :DF2DY[9,3] :=0 :
DF2DY[10,1] :=0 :DF2DY[10,3] :=0 :
DF2DY[11,1] :=0 :DF2DY[11,3] :=0 :
DF2DY[12,1] :=0 :DF2DY[12,3] :=0 :
DF2DY[13,1] :=0 :DF2DY[13,3] :=0 :
DF2DY[14,1] :=0 :DF2DY[14,3] :=0 :
DF2DY[15,1] :=0 :DF2DY[15,3] :=0 :
DF2DY[16,1] :=0 :DF2DY[16,3] :=0 :

Y: Columns 2 and 4
DF2DY[1,2] :=DFY[1,2]:DF2DY[1,4] :=0 :
DF2DY[2,2] :=DFY[2,2]:DF2DY[2,4] :=0 :
DF2DY[3,2] :=0 :DF2DY[3,4] :=DFY[1,2]:
DF2DY[4,2] :=0 :DF2DY[4,4] :=DFY[2,2]:
DF2DY[5,2] :=DFY[3,2]:DF2DY[5,4] :=0 :
DF2DY[6,2] :=DFY[4,2]:DF2DY[6,4] :=0 :
DF2DY[7,2] :=0 :DF2DY[7,4] :=DFY[3,2]:
DF2DY[8,2] :=0 :DF2DY[8,4] :=DFY[4,2]:
DF2DY[9,2] :=0 :DF2DY[9,4] :=0 :
DF2DY[10,2] :=0 :DF2DY[10,4] :=0 :
DF2DY[11,2] :=0 :DF2DY[11,4] :=0 :
DF2DY[12,2] :=0 :DF2DY[12,4] :=0 :
DF2DY[13,2] :=0 :DF2DY[13,4] :=0 :
DF2DY[14,2] :=0 :DF2DY[14,4] :=0 :
DF2DY[15,2] :=0 :DF2DY[15,4] :=0 :
DF2DY[16,2] :=0 :DF2DY[16,4] :=0 :

Y: Columns 5 and 7
DF2DY[1,5] :=DFY[1,3]:DF2DY[1,7] :=0 :
DF2DY[2,5] :=DFY[2,3]:DF2DY[2,7] :=0 :
DF2DY[3,5] :=0 :DF2DY[3,7] :=DFY[1,3]:
DF2DY[4,5] :=0 :DF2DY[4,7] :=DFY[2,3]:
DF2DY[5,5] :=DFY[3,3]:DF2DY[5,7] :=0 :
DF2DY[6,5] :=DFY[4,3]:DF2DY[6,7] :=0 :
DF2DY[7,5] :=0 :DF2DY[7,7] :=DFY[3,3]:
DF2DY[8,5] :=0 :DF2DY[8,7] :=DFY[4,3]:
DF2DY[9,5] :=0 :DF2DY[9,7] :=0 :
DF2DY[10,5] :=0 :DF2DY[10,7] :=0 :
DF2DY[11,5] :=0 :DF2DY[11,7] :=0 :
DF2DY[12,5] :=0 :DF2DY[12,7] :=0 :
DF2DY[13,5] :=0 :DF2DY[13,7] :=0 :
DF2DY[14,5] :=0 :DF2DY[14,7] :=0 :
DF2DY[15,5] :=0 :DF2DY[15,7] :=0 :
DF2DY[16,5] :=0 :DF2DY[16,7] :=0 :

Y: Columns 6 and 8
DF2DY[1,6] :=DFY[1,4]:DF2DY[1,8] :=0 :
DF2DY[2,6] :=DFY[2,4]:DF2DY[2,8] :=0 :
DF2DY[3,6] :=0 :DF2DY[3,8] :=DFY[1,4]:
DF2DY[4,6] :=0 :DF2DY[4,8] :=DFY[2,4]:
DF2DY[5,6] :=DFY[3,4]:DF2DY[5,8] :=0 :
DF2DY[6,6] :=DFY[4,4]:DF2DY[6,8] :=0 :
DF2DY[7,6] :=0 :DF2DY[7,8] :=DFY[3,4]:
DF2DY[8,6] :=0 :DF2DY[8,8] :=DFY[4,4]:
DF2DY[9,6] :=0 :DF2DY[9,8] :=0 :
DF2DY[10,6] :=0 :DF2DY[10,8] :=0 :
DF2DY[11,6] :=0 :DF2DY[11,8] :=0 :
DF2DY[12,6] :=0 :DF2DY[12,8] :=0 :
DF2DY[13,6] :=0 :DF2DY[13,8] :=0 :
DF2DY[14,6] :=0 :DF2DY[14,8] :=0 :
DF2DY[15,6] :=0 :DF2DY[15,8] :=0 :
DF2DY[16,6] :=0 :DF2DY[16,8] :=0 :

Y: Columns 9 and 11
DF2DY[1,9] :=0 :DF2DY[1,11] :=0 :
DF2DY[2,9] :=0 :DF2DY[2,11] :=0 :
DF2DY[3,9] :=0 :DF2DY[3,11] :=0 :
DF2DY[4,9] :=0 :DF2DY[4,11] :=0 :
DF2DY[5,9] :=0 :DF2DY[5,11] :=0 :
DF2DY[6,9] :=0 :DF2DY[6,11] :=0 :
DF2DY[7,9] :=0 :DF2DY[7,11] :=0 :
DF2DY[8,9] :=0 :DF2DY[8,11] :=0 :
DF2DY[9,9] :=DFY[1,1]:DF2DY[9,11] :=0 :
DF2DY[10,9] :=DFY[2,1]:DF2DY[10,11]:=0 :
DF2DY[11,9] :=0 :DF2DY[11,11]:=DFY[1,1]:
DF2DY[12,9] :=0 :DF2DY[12,11]:=DFY[2,1]:
DF2DY[13,9] :=DFY[3,1]:DF2DY[13,11]:=0 :
DF2DY[14,9] :=DFY[4,1]:DF2DY[14,11]:=0 :
DF2DY[15,9] :=0 :DF2DY[15,11]:=DFY[3,1]:
DF2DY[16,9] :=0 :DF2DY[16,11]:=DFY[4,1]:

Y: Columns 10 and 12
DF2DY[1,10] :=0 :DF2DY[1,12] :=0 :
DF2DY[2,10] :=0 :DF2DY[2,12] :=0 :
DF2DY[3,10] :=0 :DF2DY[3,12] :=0 :
DF2DY[4,10] :=0 :DF2DY[4,12] :=0 :
DF2DY[5,10] :=0 :DF2DY[5,12] :=0 :
DF2DY[6,10] :=0 :DF2DY[6,12] :=0 :
DF2DY[7,10] :=0 :DF2DY[7,12] :=0 :
DF2DY[8,10] :=0 :DF2DY[8,12] :=0 :
DF2DY[9,10] :=DFY[1,2]:DF2DY[9,12] :=0 :
DF2DY[10,10]:=DFY[2,2]:DF2DY[10,12]:=0 :
DF2DY[11,10]:=0 :DF2DY[11,12]:=DFY[1,2]:
DF2DY[12,10]:=0 :DF2DY[12,12]:=DFY[2,2]:
DF2DY[13,10]:=DFY[3,2]:DF2DY[13,12]:=0 :
DF2DY[14,10]:=DFY[4,2]:DF2DY[14,12]:=0 :
DF2DY[15,10]:=0 :DF2DY[15,12]:=DFY[3,2]:
DF2DY[16,10]:=0 :DF2DY[16,12]:=DFY[4,2]:

Y: Columns 13 and 15
DF2DY[1,13] :=0 :DF2DY[1,15] :=0 :
DF2DY[2,13] :=0 :DF2DY[2,15] :=0 :
DF2DY[3,13] :=0 :DF2DY[3,15] :=0 :
DF2DY[4,13] :=0 :DF2DY[4,15] :=0 :
DF2DY[5,13] :=0 :DF2DY[5,15] :=0 :
DF2DY[6,13] :=0 :DF2DY[6,15] :=0 :
DF2DY[7,13] :=0 :DF2DY[7,15] :=0 :
DF2DY[8,13] :=0 :DF2DY[8,15] :=0 :
DF2DY[9,13] :=DFY[1,3]:DF2DY[9,15] :=0 :
DF2DY[10,13]:=DFY[2,3]:DF2DY[10,15]:=0 :
DF2DY[11,13]:=0 :DF2DY[11,15]:=DFY[1,3]:
DF2DY[12,13]:=0 :DF2DY[12,15]:=DFY[2,3]:
DF2DY[13,13]:=DFY[3,3]:DF2DY[13,15]:=0 :
DF2DY[14,13]:=DFY[4,3]:DF2DY[14,15]:=0 :
DF2DY[15,13]:=0 :DF2DY[15,15]:=DFY[3,3]:
DF2DY[16,13]:=0 :DF2DY[16,15]:=DFY[4,3]:

Y: Columns 14 and 16
DF2DY[1,14] :=0 :DF2DY[1,16] :=0 :
DF2DY[2,14] :=0 :DF2DY[2,16] :=0 :
DF2DY[3,14] :=0 :DF2DY[3,16] :=0 :
DF2DY[4,14] :=0 :DF2DY[4,16] :=0 :
DF2DY[5,14] :=0 :DF2DY[5,16] :=0 :
DF2DY[6,14] :=0 :DF2DY[6,16] :=0 :
DF2DY[7,14] :=0 :DF2DY[7,16] :=0 :
DF2DY[8,14] :=0 :DF2DY[8,16] :=0 :
DF2DY[9,14] :=DFY[1,4]:DF2DY[9,16] :=0 :
DF2DY[10,14]:=DFY[2,4]:DF2DY[10,16]:=0 :
DF2DY[11,14]:=0 :DF2DY[11,16]:=DFY[1,4]:
DF2DY[12,14]:=0 :DF2DY[12,16]:=DFY[2,4]:
DF2DY[13,14]:=DFY[3,4]:DF2DY[13,16]:=0 :
DF2DY[14,14]:=DFY[4,4]:DF2DY[14,16]:=0 :
DF2DY[15,14]:=0 :DF2DY[15,16]:=DFY[3,4]:
DF2DY[16,14]:=0 :DF2DY[16,16]:=DFY[4,4]:

Computation of the representation of the elliptic operator, including the coefficients

The following identities will be used to simplify the formulas

HSX:={Dx*Dx=Dxx*4+Dxx*Dxx}:
HSY:={Dy*Dy=Dyy*4+Dyy*Dyy}:

The forward and backward operator have a total weight of 1, so that they truly correspond to a differentiation.
The forward difference weight is "ap", the backward difference weight is "am"

DPX:=simplify(DF2DX,{a=ap,b=1-ap}):DMX:=simplify(DF2DX,{a=am,b=1-am}):
DPY:=simplify(DF2DY,{a=ap,b=1-ap}):DMY:=simplify(DF2DY,{a=am,b=1-am}):

The final operator is then obtained by matrix multiplications and additions.

OPX:=simplify(evalm(&*(DPX,NTMX,DMX)),HSX):
OPY:=simplify(evalm(&*(DPY,NTMY,DMY)),HSY):
OPXY:=simplify(simplify(evalm(&*(DPX,NTMZ,DMY)),HSY),HSX):
OPYX:=simplify(simplify(evalm(&*(DPY,NTMZ,DMX)),HSX),HSY):

OP:=simplify(evalm(OPX+OPY+OPXY+OPYX)):

For the gaussian elimination to come, we want the coarse scale shape functions to be located at the bottom left of the matrix, whereas the wavelet component should be a the top right. Since this is the opposite of the convention adopte so far, we flip the matrix.

chgt:=matrix(16,16):
for i from 1 to 16 do for j from 1 to 16 do
chgt[i,j]:=0:
if j=17-i then chgt[i,j]:=1 fi;
od od;
#print(chgt);
OP:=multiply(chgt,OP,chgt):
# print(OP);

The result is saved in SaveFile, defined in the introduction

save MDF,OP,SaveFile:

>