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

in one dimension for the Haar Basis.

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

Definition of the Haar Wavelet Y , and of the corresponding shape function F
psi:=x->piecewise(0<=x and x<1/2, 1,1/2<=x and x<1, -1,0);
phi:=x->piecewise(0<=x and x<1,1,0);

Verification of the orthogonality of F and Y , and computation of their norms.
simplify(int(phi(x)*phi(x),x=0..1));
simplify(int(psi(x)*psi(x),x=0..1));
simplify(int(phi(x)*psi(x),x=0..1));

Verification of the interscale relations
H||0:=1:H||1:=1:G||0:=1:G||1:=-1:
MatJ:=matrix(2,2,[H0,H1,G0,G1]);
Wfs is the fine scale data
Wfs:=(x->matrix(2,1,[phi(2*x),phi(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(x));simplify(Wcs[2,1]-psi(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(1,1,[b-a]):
R||(-1):=matrix(1,1,[a]):
R||1:=matrix(1,1,[-b]):
print(R||0,R||1,R||(-1));

The functions
F and Y on the next finer scale (J+1) to compute the action of R, the difference operator.

M0 stands for the f(x) component (where f(x) is truly F ) ,
M||(-1) stands for the f(x-1) component

Ftest:=matrix(2,1,[f(x),f(x-1)]):
PreF:=simplify(multiply(MatJ,Ftest)):
M0:=transpose(matrix(2,1,[coeff(PreF[1,1],f(x)),coeff(PreF[2,1],f(x))])):
M||(-1):=transpose(matrix(2,1,[coeff(PreF[1,1],f(x-1)),coeff(PreF[2,1],f(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)):

Terms are gathered by position and re-assigned their functions.
M||a||b is applied to f(x+a+b)

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

The total difference operator is given by

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

We now return to coarse scale notations,
F and Y (noted S and P)

Invsc:=simplify(multiply(inverse(MatJ),matrix(2,1,[S(x),P(x)]))):
f(x):=Invsc[1,1]:
f(x-1):=Invsc[2,1]:
Invsc:=simplify(multiply(inverse(MatJ),matrix(2,1,[S(x+1),P(x+1)]))):
f(x+2):=Invsc[1,1]:
f(x+1):=Invsc[2,1]:
Invsc:=simplify(multiply(inverse(MatJ),matrix(2,1,[S(x-1),P(x-1)]))):
f(x-2):=Invsc[1,1]:
f(x-3):=Invsc[2,1]:
for k from 1 to 1 do
for i from 1 to 2 do
MD[1,i]:=collect(simplify(MD[1,i]),[S(x),P(x),S(x+1),P(x+1),S(x-1),P(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(2,2):A||(-1):=matrix(2,2):A||1:=matrix(2,2):
for k from 1 to 1 do for i from 1 to 2 do
A0[1,i]:=coeff(MD[1,i],S(x)):
A0[2,i]:=coeff(MD[1,i],P(x)):
A||(-1)[1,i]:=coeff(MD[1,i],S(x+1)):
A||(-1)[2,i]:=coeff(MD[1,i],P(x+1)):
A1[1,i]:=coeff(MD[1,i],S(x-1)):
A1[2,i]:=coeff(MD[1,i],P(x-1)):
od od;
for i from 1 to 2 do
for j from 1 to 2 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,phi,R0,R1,R||(-1),G0,G1,H0,H1,A0,A1,A||(-1),SaveFile;

>