Reduction of a diffusion operator from the fine scale to the coarse scale,
in two dimensions for the Haar Basis.

restart:with(linalg):
LePath:="./":
LeName:=Haar:
DataFile:=cat(LePath,"Operator",LeName,".m"):
SaveFile:=cat(LePath,"HomMat",LeName,".m"):
read DataFile:

Simplification hypothesis and simplified representation:

1. Consider only symmetric differences schemes:

am:=1-ap;

we compute the result for a true forward/backward finite difference scheme.
The parameter ap can be maintained unknown, the final result is parameter independent.

ap:=1;

We will compute taylor expansions only up to 3rd order. The simplification assumptions HS0 to HS3||3 reflect this choice. The variable X stand for the differenciation in X, whereas the variable Y stands for the differentiation in Y.

Order:=3;
HS0:={Dx=X,Dy=Y}:
HS1:={Dxx=X^2/4,Dyy=Y^2/4}:
HS3||1:={X^3=0}:
HS3||2:={Y^3=0}:
HS3||3:={Y^2*X=0}:
HS3||4:={X^2*Y=0}:

We simplify the operator according to the simplification rules we have defined.

for i from 1 to 4 do for j from 1 to 4 do
OP[i,j]:=simplify(OP[i,j],HS0):
OP[i,j]:=simplify(OP[i,j],HS1):
OP[i,j]:=mtaylor(OP[i,j],[X,Y],Order)
od:od:

The operator (here OP) is then stored in R, a 6*4*4 tensor, the following way:
for all i and j,

OP[i,j]=R[1,i,j]+R[2,i,j]*X+R[3,i,j]*Y+R[4,i,j]*X^2+R[5,i,j]*Y^2+R[6,i,j]*X*Y

R:=array(1..6,1..4,1..4);

for i from 1 to 4 do for j from 1 to 4 do
R[1,i,j]:=simplify(OP[i,j],{X=0,Y=0}):
R[2,i,j]:=simplify(coeff(OP[i,j],X),{Y=0}):
R[3,i,j]:=simplify(coeff(OP[i,j],Y),{X=0}):
R[4,i,j]:=coeff(OP[i,j],X^2):
R[5,i,j]:=coeff(OP[i,j],Y^2):
R[6,i,j]:=coeff(coeff(OP[i,j],X),Y)
od od;

save R,"FullHaar.m";

Gaussian Elimination

We then perform a gaussian elimination on the operator LM1, together with a third order taylor development in X and Y.
The resulting coefficients are stored in R, according to the same rule as before. This is of course in fact only an algebraic operation
on the coefficients of R.

For example at step k of the gaussian elimination, the (i,j) coefficient is given by

OP[i,j]:=OP[i,j]-OP[i,k]*OP[k,j]/OP[k,k]

The zeroth order coefficient is therefore given by R[1,i,j]:=R[1,i,j] - R[1,i,k]*R[1,k,j]/R[1,k,k]

and the first order coefficient in X is given by R[2,i,j]:=R[2,i,j] - (R[2,i,k]*R[1,k,j]+R[1,i,k]*R[2,k,j])/R[1,k,k]
+ R[1,i,k]*R[1,k,j]*R[2,k,k]/R[1,k,k]^2

The procedure HaarLoop performs these operations, and returns the resulting operator, after all eliminations have been performed, that is

OP[4,4]:=R[1,4,4]+R[2,4,4]*X+R[3,4,4]*Y+R[4,4,4]*X^2+R[5,4,4]*Y^2+R[6,4,4]*X*Y

HaarLoop := proc (R)
local nr, t5, t6, t48, t14, t24, t148, t146, i, t2, t40, t22, t12, j, t3, t4, t10, t20, t42, t45, t65, k, SaveFile; description "The Main Loop";
for nr to 3 do
t5 := R[1,nr,nr];
t6 := 1/t5;
t48 := t6^2;
t14 := R[2,nr,nr];
t24 := R[3,nr,nr];
t148 := t14*t48;
t146 := t6*t14;
for i from nr+1 to 4 do
printf(" nr is %d, i is %d \n",nr,i);
t2 := R[1,i,nr];
t40 := t2*t5;
t22 := R[3,i,nr];
t12 := R[2,i,nr];
for j from nr+1 to 4 do
t3 := R[1,nr,j];
t4 := t3*t2;
R[1,i,j] := normal(R[1,i,j]-t4*t6);
t10 := R[2,nr,j];
R[2,i,j] := normal(R[2,i,j]-(t2*t10+t12*t3-t4*t146)*t6);
t20 := R[3,nr,j];
R[3,i,j] := normal(-(t2*t20+t22*t3-t4*t24*t6)*t6+R[3,i,j]);
t42 := t5*t3;
t45 := -t40*t10-t42*t12+t4*t14;
R[4,i,j] := normal(R[4,i,j]-(t2*R[4,nr,j]-t4*R[4,nr,nr]*t6+R[4,i,nr]*t3+t10*t12+t45*t148)*t6);
t65 := -t40*t20-t42*t22+t24*t4;
R[5,i,j] := normal(R[5,i,j]-(t2*R[5,nr,j]-t4*R[5,nr,nr]*t6+R[5,i,nr]*t3+t22*t20+t65*t24*t48)*t6);
R[6,i,j] := normal(R[6,i,j]-(t2*R[6,nr,j]-t4*R[6,nr,nr]*t6
+R[6,i,nr]*t3+t10*t22+t12*t20+t65*t148+t45*t24*t48)*t6)
end do
end do;
for i from nr+1 to 4 do
for k to 6 do
R[k,nr,i] := 0;
R[k,i,nr] := 0
end do end do;
SaveFile := cat("Save",nr+1,".m");
# save R, SaveFile
end do;
return R[1,4,4], R[2,4,4], R[3,4,4], R[4,4,4], R[5,4,4], R[6,4,4]
end proc;

> HaarLoop(R);

The result is stored in the matrix Hom

Hom:=matrix(2,2,[16*R[4,4,4],8*R[6,4,4],8*R[6,4,4],16*R[5,4,4]]):
save Hom,SaveFile;

Simplified form of the coefficients

The simplified formulas introduced in the article refer to the case of isotropic or diagonal coefficients.
The changes of unknown CU1, CU2 and CU3 introduce the simplified variables.
CU1:={a4y=Zy+a2y*a3y/a1y,a4x=Zx+a2x*a3x/a1x};
CU2:={a3y*a3y=a1y^2-a1y*Xy,a3x*a3x=a1x^2-a1x*Xx};
CU3:={a2y*a2y=a1y^2-a1y*Yy,a2x*a2x=a1x^2-a1x*Yx};

b1:=0:b2:=0:b3:=0:b4:=0:
for i from 1 to 2 do
for j from 1 to 2 do
Hom[i,j]:=simplify(Hom[i,j]);
Hom[i,j]:=simplify(subs(CU1,Hom[i,j]));
Hom[i,j]:=simplify(Hom[i,j],CU2);
Hom[i,j]:=simplify(Hom[i,j],CU3);
Hom[i,j]:=collect(Hom[i,j],[a1x,a1y,Xx,Yy])
od od;
print(Hom);

>