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

restart:with(linalg):

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

The suffix can be replaced by "ISO" for example if the computation is for an isotropic diffusion coefficient.

Introduction of the coefficients.

If the computation concerns only a isotropic, or diagonal coefficient, the number of unknown coefficients should be reduced accordingly.
Here the script is presented in the case of a non-diagonal diffusion matrix.

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}):
print(NTMX,NTMY,NTMZ);

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(2,2):
for i from 1 to 2 do
for j from 1 to 2 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(1,1):
for i from 1 to 1 do for j from 1 to 1 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(4,4):DF2DY:=matrix(4,4):
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]:

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]:

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]:

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]:

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(multiply(multiply(DPX,NTMX),DMX),HSX):
OPY:=simplify(multiply(multiply(DPY,NTMY),DMY),HSY):
OPXY:=simplify(simplify(multiply(multiply(DPX,NTMZ),DMY),HSY),HSX):
OPYX:=simplify(simplify(multiply(multiply(DPY,NTMZ),DMX),HSX),HSY):

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(4,4):
for i from 1 to 4 do for j from 1 to 4 do
chgt[i,j]:=0:
if j=5-i then chgt[i,j]:=1 fi;
od od;
#print(chgt);

OP:=multiply(chgt,OP,chgt):
# print(OP,MDF);

The result is saved in SaveFile, defined in the introduction

save MDF,OP,SaveFile:

>