Construction of the representation of difference operators from a fine scale to a coarse scale

in one dimension for the Legendre Multiwavelet Basis.

restart:with(linalg):
LePath:="./":
LeName:=Legendre:
SaveFile:=cat(LePath,"Data",LeName,".m");

Definition of the Lagrange Multiwavelets Y0 and Y1 and of the corresponding shape functions F0 and F1

psi||0:=(x->piecewise(0<=x and x<1/2, (1-4*x)*sqrt(3),1/2<=x and x<1, (4*x-3)*sqrt(3),0)):
simplify(psi0(x));
psi||1:=(x->piecewise(0<=x and x<1/2, 6*x-1,1/2<=x and x<1, 6*x-5,0)):
simplify(psi1(x));
phi||0:=(x->piecewise(0<=x and x<1,1,0)):
simplify(phi0(x));
phi||1:=(x->piecewise(0<=x and x<1,sqrt(3)*(2*x-1),0)):
simplify(phi1(x));

Verification of the orthogonality, and computation of norms.
for i from 0 to 1 do
for j from 0 to 1 do
f||i||j:=simplify(int(phi||i(x)*phi||j(x),x=0..1));
g||i||j:=simplify(int(psi||i(x)*psi||j(x),x=0..1));
h||i||j:=simplify(int(phi||i(x)*psi||j(x),x=0..1));
v1:=cat('Phi',i," x Phi",j, " = ",f||i||j,", Psi",i," x Psi",j, " = ",g||i||j, "Phi",i," x Phi",j, " = ",h||i||j);
print(v1);

od od;

Verification of the interscale relations

H||0:=matrix(2,2,[1,0,-sqrt(3)/2,1/2]):
H||1:=matrix(2,2,[1,0,sqrt(3)/2,1/2]):
G||0:=matrix(2,2,[0,-1,1/2,sqrt(3)/2]):
G||1:=matrix(2,2,[0,1,-1/2,sqrt(3)/2]):
MatJ:=matrix(4,4):
MatJ[1,1]:=H0[1,1]:MatJ[1,2]:=H0[1,2]:MatJ[1,3]:=H1[1,1]:MatJ[1,4]:=H1[1,2]:
MatJ[2,1]:=H0[2,1]:MatJ[2,2]:=H0[2,2]:MatJ[2,3]:=H1[2,1]:MatJ[2,4]:=H1[2,2]:
MatJ[3,1]:=G0[1,1]:MatJ[3,2]:=G0[1,2]:MatJ[3,3]:=G1[1,1]:MatJ[3,4]:=G1[1,2]:
MatJ[4,1]:=G0[2,1]:MatJ[4,2]:=G0[2,2]:MatJ[4,3]:=G1[2,1]:MatJ[4,4]:=G1[2,2]:

Wfs is the fine scale data
Wfs:=(x->matrix(4,1,[phi||0(2*x),phi||1(2*x),phi||0(2*x-1),phi||1(2*x-1)])):
Wcs is the coarse scale shape function and wavelet, given by MatJ*Wfs.
Wcs:=simplify(multiply(MatJ,Wfs(x))):
Checking...
simplify(Wcs[1,1]-phi||0(x));
simplify(Wcs[2,1]-phi||1(x));
simplify(Wcs[3,1]-psi||0(x));
simplify(Wcs[4,1]-psi||1(x));

Definition of the finite difference operator
The parameters "a" and "b" corresponds to the relative weights of the forward and backward difference.

R||0:=matrix(2,2,[(b-a),(2-a-b)*sqrt(3),(-a-b)*sqrt(3),(b-a)*3]):
R||1:=matrix(2,2,[-b,-sqrt(3)*b,sqrt(3)*b,3*b]):
R||(-1):=matrix(2,2,[a,-sqrt(3)*a,sqrt(3)*a,-a*3]):
print(R||(-1),R||0,R||1);

The functions
F0,F1 and Y0,Y1 on the next finer scale (J+1) to compute the action of R, the difference operator.
M0 stands for the f0(x) and f1(x) component (where f||k(x) is truly
F||k ) ,
M||(-1) stands for the f0(x-1) and f1(x-1) component

Ftest:=matrix(4,1,[f||0(x),f||1(x),f||0(x-1),f||1(x-1)]):
PreF:=simplify(multiply(MatJ,Ftest)):
M0:=transpose(matrix(4,2,[coeff(PreF[1,1],f0(x)),coeff(PreF[1,1],f1(x)),coeff(PreF[2,1],f0(x)),
coeff(PreF[2,1],f1(x)),coeff(PreF[3,1],f0(x)),coeff(PreF[3,1],f1(x)),coeff(PreF[4,1],f0(x)),
coeff(PreF[4,1],f1(x))])):
M||(-1):=transpose(matrix(4,2,[coeff(PreF[1,1],f0(x-1)),coeff(PreF[1,1],f1(x-1)),coeff(PreF[2,1],f0(x-1)),
coeff(PreF[2,1],f1(x-1)),coeff(PreF[3,1],f0(x-1)),coeff(PreF[3,1],f1(x-1)),coeff(PreF[4,1],f0(x-1)),
coeff(PreF[4,1],f1(x-1))])):

M||a||b is the difference operator in x+b applied to the x+a term

M||0||0:=multiply(R||0,M||0):
M||0||1:=multiply(R||(-1),M0):
M||0||(-1):=multiply(R||1,M0):

M||(-1)||0:=multiply(R||0,M||(-1)):
M||(-1)||1:=multiply(R||(-1),M||(-1)):
M||(-1)||(-1):=multiply(R||1,M||(-1)):

Md1:=multiply(matrix(1,2,[f0(x+1),f1(x+1)]),M01):
Md0:=multiply(matrix(1,2,[f0(x),f1(x)]),evalm(M00+M||(-1)||1)):
Md||(-1):=multiply(matrix(1,2,[f0(x-1),f1(x-1)]),evalm(M||(-1)||0+M||0||(-1))):
Md||(-2):=multiply(matrix(1,2,[f0(x-2),f1(x-2)]),M||(-1)||(-1)):

The total difference operator is given by

MD:=evalm(Md1+Md0+Md||(-1)+Md||(-2)):

F0,F1 and Y0, Y1 (noted S and P)
Invsc:=simplify(multiply(inverse(MatJ),
matrix(4,1,[S0(x),S1(x),P0(x),P1(x)]))):
f0(x):=Invsc[1,1]:
f1(x):=Invsc[2,1]:
f0(x-1):=Invsc[3,1]:
f1(x-1):=Invsc[4,1]:
Invsc:=simplify(multiply(inverse(MatJ),
matrix(4,1,[S0(x+1),S1(x+1),P0(x+1),P1(x+1)]))):
f0(x+2):=Invsc[1,1]:
f1(x+2):=Invsc[2,1]:
f0(x+1):=Invsc[3,1]:
f1(x+1):=Invsc[4,1]:
Invsc:=simplify(multiply(inverse(MatJ),
matrix(4,1,[S0(x-1),S1(x-1),P0(x-1),P1(x-1)]))):
f0(x-2):=Invsc[1,1]:
f1(x-2):=Invsc[2,1]:
f0(x-3):=Invsc[3,1]:
f1(x-3):=Invsc[4,1]:

for k from 1 to 1 do
for i from 1 to 4 do MD[1,i]:=collect(simplify(MD[1,i]),[S0(x),S1(x),P0(x),P1(x),S0(x+1),S1(x+1),P0(x+1),P1(x+1),S0(x-1),S1(x-1),P0(x-1),P1(x-1)]) od; od;

The resulting difference operators are stored in three matrices,
A0 for the "x " part,

A1 for the "x-1 " part,

A||(-1) for the "x+1" part.

A0:=matrix(4,4):A||(-1):=matrix(4,4):A||1:=matrix(4,4):
for k from 1 to 1 do for i from 1 to 4 do
A0[1,i]:=coeff(MD[1,i],S0(x)):
A0[2,i]:=coeff(MD[1,i],S1(x)):
A0[3,i]:=coeff(MD[1,i],P0(x)):
A0[4,i]:=coeff(MD[1,i],P1(x)):
A||(-1)[1,i]:=coeff(MD[1,i],S0(x+1)):
A||(-1)[2,i]:=coeff(MD[1,i],S1(x+1)):
A||(-1)[3,i]:=coeff(MD[1,i],P0(x+1)):
A||(-1)[4,i]:=coeff(MD[1,i],P1(x+1)):
A1[1,i]:=coeff(MD[1,i],S0(x-1)):
A1[2,i]:=coeff(MD[1,i],S1(x-1)):
A1[3,i]:=coeff(MD[1,i],P0(x-1)):
A1[4,i]:=coeff(MD[1,i],P1(x-1)):
od od;
for i from 1 to 4 do
for j from 1 to 4 do
for k from (-1) to 1 do
A||k[i,j]:= factor(simplify(A||k[i,j])) od od od;
print("1,0,-1 ", A1,A0,A||(-1));
print("J 1,0,-1 ", R||1,R0,R||(-1));

The results are stored in SaveFile (defined in the introduction).

save psi||0,psi||1,phi||0,phi||1,R0,R1,R||(-1),G0,G1,H0,H1,A0,A1,A||(-1),SaveFile;

>