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

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

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 is taken to be close to 1/2; a taylor expansion around 1/2 up to order 3 will be performed.

ap:=1/2+T;

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 16 do for j from 1 to 16 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):
OP[i,j]:=mtaylor(OP[i,j],T,Order)
od:od:

The operator (here OP) is then stored in R, a 16*16*6*3 tensor, the following way, corresponding to
the asymptotic development explained above:
for all i and j,

OP[i,j]= R[i,j,1,1]+T*R[i,j,1,2]+T^2*R[i,j,1,3]
+(R[i,j,2,1]+T*R[i,j,2,2]+T^2*R[i,j,2,3])*X
+(R[i,j,3,1]+T*R[i,j,3,2]+T^2*R[i,j,3,3])*Y
+(R[i,j,4,1]+T*R[i,j,4,2]+T^2*R[i,j,4,3])*X^2
+(R[i,j,5,1]+T*R[i,j,5,2]+T^2*R[i,j,5,3])*Y^2
+(R[i,j,6,1]+T*R[i,j,6,2]+T^2*R[i,j,6,3])*X*Y
In fact, the singularity around 1/2 with respect to T is of order 2
because of the symmetry of the problem, there is no need to keep
track of that many terms, only the following are necessary

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

R:=array(1..16,1..16,1..6,1..3);

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

save R,"FullLegendre.m";

Gaussian Elimination

We then perform a gaussian elimination on the operator LM1, together with a third order taylor development in X and Y, and T.
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 (in X, Y and T) is therefore given by R[i,j,1,1]:=R[i,j,1,1] - R[i,k,1,1]*R[k,j,1,1]/R[k,k,1,1]

The procedure MultiWaveletTaylor performs these operations, for each given step.
The procedure MultiWaveletLoop is the gaussian elimination for the first 12 steps.
For the last four steps (the reduction within the scaling vectors),
the taylor developpment in T is omitted (all terms are gathered together), and the procedure used for the Haar Reduction procedure is used.
HaarLoop returns the resulting reduced operator, after all eliminations have been performed, that is

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

MultiWaveletTaylor := proc(R, i, j, n)
local t30, t31, t48, t49, t50, t51, t42, t40, t41, t38, t39,
t67, t68, t69, t66, t161, t61, t62, t54, t52, t53, t70, t81,
t77, t78, t79, t80, t76, t82, t96, t89, t90, t99, t97, t102,
t103, t104, t158, t159, t141, t164, t165, t166, t5, t6, t4,
t3, t2, result, t220, t217, t218, t221, t225, t227, t9, t12,
t11, t226, t10, t13, t17;
result := array(1 .. 6, 1 .. 3);
t5 := R[n, n, 1, 1];
t6 := 1/t5;
t2 := R[n, j, 1, 1];
t3 := R[i, n, 1, 1];
t4 := t2*t3;
result[1, 1] := R[i, j, 1, 1] - t4*t6;
t9 := R[i, n, 1, 2];
t10 := t2*t9;
t11 := R[n, j, 1, 2];
t12 := t11*t3;
t13 := R[n, n, 1, 2];
t17 := (t10 + t12 - t4*t6*t13)*t6;
result[1, 2] := R[i, j, 1, 2] - t17;
t30 := t5^2;
t31 := 1/t30;
result[1, 3] := -(t2*R[i, n, 1, 3] + t11*t9
+ R[n, j, 1, 3]*t3 - t4*t6*R[n, n, 1, 3]
- (t10*t5 + t12*t5 - t4*t13)*t31*t13)*t6
+ R[i, j, 1, 3];
t38 := R[i, n, 2, 1];
t39 := t2*t38;
t40 := R[n, j, 2, 1];
t41 := t40*t3;
t42 := R[n, n, 2, 1];
result[2, 1] :=
R[i, j, 2, 1] - (t39 + t41 - t4*t6*t42)*t6;
t48 := R[i, n, 2, 2];
t49 := t2*t48;
t50 := t11*t38;
t51 := t40*t9;
t52 := R[n, j, 2, 2];
t53 := t52*t3;
t54 := R[n, n, 2, 2];
t61 := t39*t5 + t41*t5 - t4*t42;
t62 := t61*t31;
result[2, 2] := R[i, j, 2, 2] - (
t49 + t50 + t51 + t53 - t4*t6*t54 - t17*t42 - t62*t13
)*t6;
result[2, 3] := 0;
t66 := R[i, n, 3, 1];
t67 := t2*t66;
t68 := R[n, j, 3, 1];
t69 := t68*t3;
t70 := R[n, n, 3, 1];
result[3, 1] :=
-(t67 + t69 - t4*t6*t70)*t6 + R[i, j, 3, 1];
t76 := R[i, n, 3, 2];
t77 := t2*t76;
t78 := t11*t66;
t79 := t68*t9;
t80 := R[n, j, 3, 2];
t81 := t80*t3;
t82 := R[n, n, 3, 2];
t89 := t67*t5 + t69*t5 - t4*t70;
t90 := t89*t31;
result[3, 2] := -(
t77 + t78 + t79 + t81 - t4*t6*t82 - t17*t70 - t90*t13
)*t6 + R[i, j, 3, 2];
result[3, 3] := 0;
t96 := R[i, n, 4, 1];
t97 := t2*t96;
t99 := R[n, n, 4, 1];
t102 := R[n, j, 4, 1];
t103 := t102*t3;
t104 := t40*t38;
result[4, 1] := R[i, j, 4, 1]
- (t97 - t62*t42 - t4*t6*t99 + t103 + t104)*t6;
result[4, 2] := 0;
result[4, 3] := 0;
t158 := R[i, n, 5, 1];
t159 := t2*t158;
t161 := R[n, n, 5, 1];
t164 := R[n, j, 5, 1];
t165 := t164*t3;
t166 := t68*t66;
result[5, 1] := R[i, j, 5, 1]
- (t159 - t90*t70 - t4*t6*t161 + t165 + t166)*t6;
result[5, 2] := 0;
result[5, 3] := 0;
t217 := R[i, n, 6, 1];
t218 := t2*t217;
t220 := t40*t66;
t221 := R[n, n, 6, 1];
t225 := R[n, j, 6, 1];
t226 := t225*t3;
t227 := t68*t38;
result[6, 1] := R[i, j, 6, 1] - (t218 - t90*t42 + t220
- t4*t6*t221 - t62*t70 + t226 + t227)*t6;
result[6, 2] := 0;
result[6, 3] := 0;
result
end proc;

> MultiWaveletLoop:=proc(R)
local nr,i,j,k,l,inter, SaveFile;
description "The Main Loop";
for nr to 12 do
for i from nr+1 to 16 do
printf(" nr is %d, i is %d \n",nr,i);
inter:=array(1..6,1..3):
for j from nr+1 to 16 do
inter:=MultiWaveletTaylor(R,i,j,nr):
for k from 1 to 6 do
for l from 1 to 3 do
R[i,j,k,l]:=normal(inter[k,l]):
end do
end do;
end do
end do;
for i from nr+1 to 16 do
for k to 6 do
for l to 3 do
R[nr,i,k,l] := 0:
R[nr,i,k,l] := 0:
end do end do
end do;
SaveFile := cat("Save",nr+1,".m");
# save R, SaveFile
end do;
return R
end proc;

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;

> MultiWaveletLoop(R);

> R2:=array(1..6,1..4,1..4):

> for i from 13 to 16 do
for j from 13 to 16 do
for k from 1 to 6 do
R2[k,i-12,j-12]:=normal(R[i,j,k,1])+T*normal(R[i,j,k,2]) +T^2*normal(R[i,j,k,3]):
end do;
end do;
end do;

> HaarLoop(R2);

The result is stored in the matrix Hom

Hom:=matrix(2,2,[16*R2[4,4,4],8*R2[6,4,4],8*R2[6,4,4],16*R2[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:={a4=Z+a2*a3/a1};
CU2:={a3*a3=a1^2-a1*X};
CU3:={a2*a2=a1^2-a1*Y};

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],[a1,X,Y])
od od;
T1:=X+simplify(Hom[1,1]-Hom[1,2]-X):
T2:=Y+simplify(Hom[2,2]-Hom[1,2]-Y):
Hom[1,2]:=simplify(Hom[1,2]*(T1+T2))/(T1+T2):
Hom[2,1]:=Hom[1,2]:
Hom[1,1]:=T1+Hom[1,2]:
Hom[2,2]:=T2+Hom[1,2]:
print(Hom);

>