Symbolic local Fourier analysis for determining an approximation error estimate

This Mathematica notebook accompanies the paper

	“Using cylindrical algebraic decomposition and local Fourier analysis to analyze numerical methods: two examples” by S. Takacs

 A preprint version is available at

	http://www.numa.uni-linz.ac.at/~stefant/J3362/slfa/

Case 1: The Courant element

First, we define the basis functions:

(Subscript[φ, k_,i_])[x_]:=If[(Subscript[h, k]) (i-1)<x≤(Subscript[h, k]) i,x/(Subscript[h, k]) -(i-1),0]+If[(Subscript[h, k]) i<x≤(Subscript[h, k]) (i+1),-x/(Subscript[h, k]) +(i+1),0]

Plot[(Subscript[φ, 2,2])[x]/.(Subscript[h, k_]):>(2^(-k)),{x,0,1.5}]

The next step is to compute the integrals that determine the mass matrix. As the support of the basis functions consists of two elements, we obtain a tri-diagonal matrix, where the diagonal entries and the off-diagonal entries have the following values:

Integrate[(((Subscript[φ, k,i])[x])^2),{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[(Subscript[φ, k,i])[x]*(Subscript[φ, k,i+1])[x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

So, we obtain that the mass matrix (Subscript[M, k]) has the following tri-diagonal form:

(Subscript[M, k])=((Subscript[h, k])/6)(({{4, 1,  ,  ,   ,  ,  ,  ,  ,  }, {1, 4, 1,  ,  ,  ,  ,  ,  ,  }, { , 1, 4, 1,  ,  ,  ,  ,  ,  }, { ,  , 1, 4, 1,  ,  ,  ,  ,  }, { ,  ,  , 1, 4, 1,  ,  ,  ,  }, { ,  ,  ,  , 1, 4, 1,  ,  ,  }, { ,  ,  ,  ,  , 1, 4, 1,  ,  }, { ,  ,  ,  ,  ,  , 1, 4, 1,  }, { ,  ,  ,  ,  ,  ,  , 1, 4, 1}, { ,  ,  ,  ,  ,  ,  ,  , 1, 4}}))

By multiplying (Subscript[M, k]) with the complex exponentials (Subscript[(UnderBar[ϕ]), k])(θ) =(Subscript[((e^(i θ j))), j]) obtain in the j-th row:
	((Subscript[h, k])/6)(1*(e^(i θ (j-1)))+4*(e^(i θ j))+1*(e^(i θ (j+1))))=((Subscript[h, k])/6)(4+(e^(-i θ))+(e^(i θ )))(e^(i θ j)).

So, (Subscript[M, k])(Subscript[(UnderBar[ϕ]), k])(θ)  =((Subscript[h, k])/6)(4+(e^(-i θ))+(e^(i θ )))(Subscript[(UnderBar[ϕ]), k])(θ)  .

So, the symbol reads as follows:

(Subscript[M, k_])[θ_]:=((Subscript[h, k])/6)(4+(E^(-I θ))+(E^(I θ)))

Now we compute the integrals that determine the stiffness matrix. As the support of the basis functions consists of two elements, we obtain a tri-diagonal matrix, where the diagonal entries and the off-diagonal entries have the following values:

Integrate[((D[(Subscript[φ, k,i])[x],x])^2),{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[D[(Subscript[φ, k,i])[x],x]D[(Subscript[φ, k,i+1])[x],x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

So, we obtain that the stiffness matrix (Subscript[K, k]) has the following tri-diagonal form:

 (Subscript[K, k])=(1/(Subscript[h, k]))(({{2, -1,  ,  ,   ,  ,  ,  ,  ,  }, {-1, 2, -1,  ,  ,  ,  ,  ,  ,  }, { , -1, 2, -1,  ,  ,  ,  ,  ,  }, { ,  , -1, 2, -1,  ,  ,  ,  ,  }, { ,  ,  , -1, 2, -1,  ,  ,  ,  }, { ,  ,  ,  , -1, 2, -1,  ,  ,  }, { ,  ,  ,  ,  , -1, 2, -1,  ,  }, { ,  ,  ,  ,  ,  , -1, 2, -1,  }, { ,  ,  ,  ,  ,  ,  , -1, 2, -1}, { ,  ,  ,  ,  ,  ,  ,  , -1, 2}}))

The symbol reads as follows:

(Subscript[K, k_])[θ_]:=(1/(Subscript[h, k]))(2-(E^(-I θ))-(E^(I θ)))

The next step is to determine the intergrid-taransfer matrix, as outlined in the paper.

(Subscript[ϕ, k_])[θ_,x_]:=(∑_(i=-∞)^∞)(Subscript[φ, k,i])[x]Exp[i θ I]

Here, we have to solve the equations (11) and (13) from the paper:

Simplify[(Subscript[ϕ, k-1])[2θ,0] == A (Subscript[ϕ, k])[θ,0]+B (Subscript[ϕ, k])[θ+π,0],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

Simplify[(Subscript[ϕ, k-1])[2θ,(Subscript[h, k])] == A (Subscript[ϕ, k])[θ,(Subscript[h, k])]+B (Subscript[ϕ, k])[θ+π,(Subscript[h, k])],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

Solve[{%,%%}]

{A,B}/.(%[[1]])//Expand

P[θ_]:={{(1/2)+((E^(-I θ))/4)+((E^(I θ))/4), (1/2)-((E^(-I θ))/4)-((E^(I θ))/4)}}

As outlined in the paper, we have to set up also the symbols for the mass matrix and the stiffnes matrix for the two-dimensional basis:

(Subscript[M2, k_])[θ_]:=DiagonalMatrix[{(Subscript[M, k])[θ],(Subscript[M, k])[θ+π]}]

(Subscript[K2, k_])[θ_]:=DiagonalMatrix[{(Subscript[K, k])[θ],(Subscript[K, k])[θ+π]}]

Now, we check if the Galerkin identity is satisfied (up to scaling):

FullSimplify[P[θ].(Subscript[M2, k])[θ].Transpose[P[θ]]==(1/2) (Subscript[M, k-1])[2θ],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

FullSimplify[P[θ].(Subscript[K2, k])[θ].Transpose[P[θ]]==(1/2) (Subscript[K, k-1])[2θ],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

The symbol of the whole operator (Subscript[G, k]) =(1/((Subscript[h, k])^2))( I - (Subscript[Π, k]) ) ((Subscript[K, k])^(-1))(Subscript[M, k])  reads as follows:

(Subscript[G, k_])[θ_]=FullSimplify[(IdentityMatrix[2]-Transpose[P[θ]].Inverse[(1/2){{(Subscript[K, k-1])[2θ]}}].P[θ].(Subscript[K2, k])[θ]).Inverse[(Subscript[K2, k])[θ]].(Subscript[M2, k])[θ]/((Subscript[h, k])^2),Assumptions->{(Subscript[h, k-1])==2(Subscript[h, k])}]

As this symbol has rank 1, the spectral radius is equal to to the sum of the eigenvalues, which coincides with the trace of the matrix:

Simplify[(Subscript[G, k])[θ][[1,1]]+(Subscript[G, k])[θ][[2,2]]]

In this case, we do not have to determine the supremum anymore, because the spectral radius takes the value (1/3) for all frequencies θ.

Case 2: A (P^2)-spline discretization

First, we define the basis functions:

(Subscript[φ, k_,i_])[x_]:=If[(Subscript[h, k]) (i-1)<x≤(Subscript[h, k]) i,(1/(2((Subscript[h, k])^2)))(((x-(Subscript[h, k])(i-1)))^2),0]+If[(Subscript[h, k]) i<x≤(Subscript[h, k]) (i+1),(3/4)-(1/(4 ((Subscript[h, k])^2)))(((2x-(Subscript[h, k]) i-(Subscript[h, k])(i+1)))^2),0]+If[(Subscript[h, k]) (i+1)<x≤(Subscript[h, k]) (i+2),(1/(2((Subscript[h, k])^2)))(((x-(Subscript[h, k])(i+2)))^2),0]

Plot[(Subscript[φ, 2,2])[x]/.(Subscript[h, k_]):>(2^(-k)),{x,0,2},PlotRange->{0,1}]

The next step is to compute the integrals that determine the mass matrix. As the support of the basis functions consists of three elements, we obtain a multi-diagonal matrix, where the diagonal entries and the off-diagonal entries have the following values:

Integrate[(((Subscript[φ, k,i])[x])^2),{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[(Subscript[φ, k,i])[x]*(Subscript[φ, k,i+1])[x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[(Subscript[φ, k,i])[x]*(Subscript[φ, k,i+2])[x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

So, we obtain that the mass matrix (Subscript[M, k]) has the following form:

(Subscript[M, k])=((Subscript[h, k])/120)(({{66, 26, 1,  ,   ,  ,  ,  ,  ,  }, {26, 66, 26, 1 ,  ,  ,  ,  ,  ,  }, { 1, 26, 66, 26, 1 ,  ,  ,  ,  ,  }, { ,  1, 26, 66, 26, 1 ,  ,  ,  ,  }, { ,  ,  1, 26, 66, 26, 1 ,  ,  ,  }, { ,  ,  ,  1, 26, 66, 26, 1 ,  ,  }, { ,  ,  ,  ,  1, 26, 66, 26, 1 ,  }, { ,  ,  ,  ,  ,  1, 26, 66, 26, 1 }, { ,  ,  ,  ,  ,  ,  1, 26, 66, 26}, { ,  ,  ,  ,  ,  ,  , 1 , 26, 66}}))

The symbol reads as follows:

(Subscript[M, k_])[θ_]:=((Subscript[h, k])/120)(66+26(E^(-I θ))+26(E^(I θ))+(E^(-2I θ))+(E^(2I θ)))

Now we compute the integrals that determine the stiffness matrix. As the support of the basis functions consists of three elements, we obtain a multi-diagonal matrix, where the diagonal entries and the off-diagonal entries have the following values:

Integrate[((D[(Subscript[φ, k,i])[x],x])^2),{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[D[(Subscript[φ, k,i])[x],x]D[(Subscript[φ, k,i+1])[x],x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[D[(Subscript[φ, k,i])[x],x]D[(Subscript[φ, k,i+2])[x],x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

So, we obtain that the stiffness matrix (Subscript[K, k]) has the following form:

 (Subscript[K, k])=(1/(6(Subscript[h, k])))(({{6, -2,  -1,  ,   ,  ,  ,  ,  ,  }, {-2, 6, -2,  -1,  ,  ,  ,  ,  ,  }, { -1, -2, 6, -2, -1 ,  ,  ,  ,  ,  }, { , -1 , -2, 6, -2, -1 ,  ,  ,  ,  }, { ,  , -1 , -2, 6, -2, -1 ,  ,  ,  }, { ,  ,  , -1 , -2, 6, -2, -1 ,  ,  }, { ,  ,  ,  , -1 , -2, 6, -2, -1 ,  }, { ,  ,  ,  ,  , -1 , -2, 6, -2, -1 }, { ,  ,  ,  ,  ,  , -1 , -2, 6, -2}, { ,  ,  ,  ,  ,  ,  , -1 , -2, 6}}))

The symbol reads as follows:

(Subscript[K, k_])[θ_]:=(1/(6(Subscript[h, k])))(6-2(E^(-I θ))-2(E^(I θ))-(E^(-2I θ))-(E^(2I θ)))

The next step is to determine the intergrid-taransfer matrix, as outlined in the paper:

L=10;

(Subscript[ϕ, k_])[θ_,x_]:=(∑_(i=-L)^L)(Subscript[φ, k,i])[x]Exp[i θ I]

Here, we have to solve the equations (11) and (13) from the paper:

Simplify[(Subscript[ϕ, k-1])[2θ,0] == A (Subscript[ϕ, k])[θ,0]+B (Subscript[ϕ, k])[θ+π,0],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

Simplify[(Subscript[ϕ, k-1])[2θ,(Subscript[h, k])] == A  (Subscript[ϕ, k])[θ,(Subscript[h, k])]+B (Subscript[ϕ, k])[θ+π,(Subscript[h, k])],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

Solve[{%,%%},{A,B}]

{A,B}/.(%[[1]])//Expand

P[θ_]:={{(3/8)+((3 (E^(-I θ)))/8)+((E^(I θ))/8)+(1/8) (E^(-2 I θ)),(3/8)-((3 (E^(-I θ)))/8)-((E^(I θ))/8)+(1/8) (E^(-2 I θ))}}

As outlined in the paper, we have to set up also the symbols for the mass matrix and the stiffnes matrix for the two-dimensional basis:

(Subscript[M2, k_])[θ_]:=DiagonalMatrix[{(Subscript[M, k])[θ],(Subscript[M, k])[θ+π]}]

(Subscript[K2, k_])[θ_]:=DiagonalMatrix[{(Subscript[K, k])[θ],(Subscript[K, k])[θ+π]}]

Now we check if the Galerkin identity is satisfied (up to scaling):

FullSimplify[P[θ].(Subscript[M2, k])[θ].Conjugate[Transpose[P[θ]]]=={{(1/2) (Subscript[M, k-1])[2θ]}},Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k]),θ∈Reals}]

FullSimplify[P[θ].(Subscript[K2, k])[θ].Conjugate[Transpose[P[θ]]]=={{(1/2) (Subscript[K, k-1])[2θ]}},Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k]),θ∈Reals}]

The symbol of the whole operator (Subscript[G, k]) =(1/((Subscript[h, k])^2))( I - (Subscript[Π, k]) ) ((Subscript[K, k])^(-1))(Subscript[M, k])  reads as follows:

(Subscript[G, k_])[θ_]=FullSimplify[(IdentityMatrix[2]-Conjugate[Transpose[P[θ]]].Inverse[(1/2){{(Subscript[K, k-1])[2θ]}}].P[θ].(Subscript[K2, k])[θ]).Inverse[(Subscript[K2, k])[θ]].(Subscript[M2, k])[θ]/((Subscript[h, k])^2),Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k]),θ∈Reals}]

As this symbol has rank 1, the spectral radius is equal to to the sum of the eigenvalues, which coincides with the trace of the matrix:

spectralradius=Simplify[(Subscript[G, k])[θ][[1,1]]+(Subscript[G, k])[θ][[2,2]]]

Now, we rewrite this term as rational function by replacing cos(θ) by c:

spectralradius=spectralradius/.Cos[x_]:>ChebyshevT[x/θ,c]

Here, we can use CAD to determine the supremum:

Resolve[ForAll[c,-1≤c≤1,-λ≤spectralradius≤λ]]

So, we obtain that the supremum is (2/5).

Case 3: A stanadard (P^2) discretization

First, we define the basis functions:

(Subscript[φ, k_,i_,1])[x_]:=If[(Subscript[h, k]) (i-1)<x≤(Subscript[h, k]) i,(2/((Subscript[h, k])^2))(x-(Subscript[h, k])(i-1))(x-(Subscript[h, k])(i-1/2)),0]+If[(Subscript[h, k]) i<x≤(Subscript[h, k]) (i+1),(2/((Subscript[h, k])^2))(x-(Subscript[h, k])(i+1))(x-(Subscript[h, k])(i+1/2)),0]

(Subscript[φ, k_,i_,2])[x_]:=If[(Subscript[h, k]) i<x≤(Subscript[h, k]) (i+1),-(4/((Subscript[h, k])^2))(x-(Subscript[h, k]) i)(x-(Subscript[h, k])(i+1)),0]

(Subscript[φ, k_,i_])[x_]:=If[Mod[i,2]==0,(Subscript[φ, k,i/2,1])[x],(Subscript[φ, k,(i-1)/2,2])[x]]

Plot[(Subscript[φ, 2,3])[x]/.(Subscript[h, k_]):>(2^(-k)),{x,0,1.5},PlotRange->{-.25,1}]

Plot[(Subscript[φ, 2,4])[x]/.(Subscript[h, k_]):>(2^(-k)),{x,0,1.5},PlotRange->{-.25,1}]

The next step is to compute the integrals that determine the mass matrix. Also here, we obtain a multi-diagonal matrix. However, as there are two kinds of basis functions, the coefficients are alternating:

Integrate[(((Subscript[φ, k,i])[x])^2),{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[(Subscript[φ, k,i])[x]*(Subscript[φ, k,i+1])[x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[(Subscript[φ, k,i])[x]*(Subscript[φ, k,i+2])[x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

We obtain that the mass matrix (Subscript[M, k]) has the following form:

(Subscript[M, k])=((Subscript[h, k])/30)(({{8, 2, -1,  ,   ,  ,  ,  ,  ,  }, {2, 16, 2, 0 ,  ,  ,  ,  ,  ,  }, { -1, 2, 8, 2, -1 ,  ,  ,  ,  ,  }, { ,  0, 2, 16, 2, 0 ,  ,  ,  ,  }, { ,  ,  -1, 2, 8, 2, -1 ,  ,  ,  }, { ,  ,  ,  0, 2, 16, 2, 0,  ,  }, { ,  ,  ,  ,  -1, 2, 8, 2, -1 ,  }, { ,  ,  ,  ,  ,  0, 2, 16, 2, 0 }, { ,  ,  ,  ,  ,  ,  -1, 2, 8, 2}, { ,  ,  ,  ,  ,  ,  , 0 , 2, 16}}))=((Subscript[h, k])/60)(({{24, 4, -1,  ,   ,  ,  ,  ,  ,  }, {4, 24, 4, -1 ,  ,  ,  ,  ,  ,  }, { -1, 4, 24, 4, -1 ,  ,  ,  ,  ,  }, { ,  -1, 4, 24, 4, -1 ,  ,  ,  ,  }, { ,  ,  -1, 4, 24, 4, -1 ,  ,  ,  }, { ,  ,  ,  -1, 4, 24, 4, -1,  ,  }, { ,  ,  ,  ,  -1, 4, 24, 4, -1 ,  }, { ,  ,  ,  ,  , -1, 4, 24, 4, -1 }, { ,  ,  ,  ,  ,  ,  -1, 4, 24, 4}, { ,  ,  ,  ,  ,  ,  , -1 , 4, 24}}))+((Subscript[h, k])/60)(({{-8, 0, -1,  ,   ,  ,  ,  ,  ,  }, {0, 8, 0, 1 ,  ,  ,  ,  ,  ,  }, { -1, 0, -8, 0, -1 ,  ,  ,  ,  ,  }, { ,  1, 0, 8, 0, 1 ,  ,  ,  ,  }, { ,  ,  -1, 0, -8, 0, -1 ,  ,  ,  }, { ,  ,  ,  1, 0, 8, 0, 1,  ,  }, { ,  ,  ,  ,  -1, 0, -8, 0, -1 ,  }, { ,  ,  ,  ,  , 1, 0, 8, 0, 1 }, { ,  ,  ,  ,  ,  ,  -1, 0, -8, 0}, { ,  ,  ,  ,  ,  ,  , 1 , 0, 8}}))=:(Subscript[A, k])+(Subscript[B, k])

By multiplying (Subscript[A, k]) with the complex exponentials (Subscript[(UnderBar[ϕ]), k])(θ) =(Subscript[((e^(i θ j))), j]), we obtain in the j-th row:
	((Subscript[h, k])/60)(-1*(e^(i θ (j-2)))+4*(e^(i θ (j-1)))+24*(e^(i θ j))+4*(e^(i θ (j+1)))-1*(e^(i θ (j+2))))=((Subscript[h, k])/60)(24+4(e^(-i θ))+4(e^(i θ ))-(e^(-2i θ))-4(e^(2i θ )))(e^(i θ j)).

So, (Subscript[A, k])(Subscript[(UnderBar[ϕ]), k])(θ)  =(Underscript[((Subscript[h, k])/60)(24+4(e^(-i θ))+4(e^(i θ ))-(e^(-2i θ))-4(e^(2i θ ))), = :(Overscript[(Subscript[A, k]), ∧])(θ)]) (Subscript[(UnderBar[ϕ]), k])(θ)

By multiplying (Subscript[B, k]) with the complex exponentials (Subscript[(UnderBar[ϕ]), k])(θ) =(Subscript[((e^(i θ j))), j]), we obtain in the j-th row:
	(h/60)((((-1))^(j+1))*(e^(i θ (j-2)))+(((-1))^(j+1))*8*(e^(i θ j))+(((-1))^(j+1))*(e^(i θ (j+2))))=(h/60)(-(e^(-2i θ))-8-(e^(2i θ )))(e^(i (θ +π)j)).

So, (Subscript[B, k])(Subscript[(UnderBar[ϕ]), k])(θ) =(Underscript[((Subscript[h, k])/60)(-(e^(-2i θ))-8-(e^(2i θ ))), =:(Overscript[(Subscript[B, k]), ∧])(θ)])(Subscript[(UnderBar[ϕ]), k])(θ+π)

Concluding, we obtain

	(Subscript[M, k])(Subscript[(UnderBar[ϕ]), k])(θ) = (Overscript[(Subscript[A, k]), ∧])(θ)(Subscript[ϕ, k])(θ) + (Overscript[(Subscript[B, k]), ∧])(θ)(Subscript[(UnderBar[ϕ]), k])(θ+π)

So, ob observe that (Subscript[M, k])(Subscript[(UnderBar[ϕ]), k])(θ) is not an element of span{(Subscript[(UnderBar[ϕ]), k])(θ)}, but of span{(Subscript[(UnderBar[ϕ]), k])(θ),(Subscript[(UnderBar[ϕ]), k])(θ+π)}.

So, the symbol has to be formulated in the basis ((Subscript[(UnderBar[ϕ]), k])(θ),(Subscript[(UnderBar[ϕ]), k])(θ+π)).

As we obviously also have

	(Subscript[M, k])(Subscript[(UnderBar[ϕ]), k])(θ++π) = (Overscript[(Subscript[B, k]), ∧])(θ++π)(Subscript[(UnderBar[ϕ]), k])(θ)+ (Overscript[(Subscript[A, k]), ∧])(θ++π)(Subscript[(UnderBar[ϕ]), k])(θ++π) ,

 we see immediately that

	(Overscript[(Subscript[M, k]), ∧])(θ)=(({{(Overscript[(Subscript[A, k]), ∧])(θ), (Overscript[(Subscript[B, k]), ∧])(θ)}, {(Overscript[(Subscript[B, k]), ∧])(θ+π), (Overscript[(Subscript[A, k]), ∧])(θ+π)}})).
	
in the basis ((Subscript[(UnderBar[ϕ]), k])(θ),(Subscript[(UnderBar[ϕ]), k])(θ+π)). 

So, we define as follows:

(Subscript[M, k_])[θ_]:=((Subscript[h, k])/60)(({{24+4 (E^(-I θ))+4(E^(I θ))-(E^(-2I θ))-(E^(2I θ)), -8-(E^(-2I θ))-(E^(2I θ))}, {-8-(E^(-2I (θ+π)))-(E^(2I (θ+π))), 24+4 (E^(-I (θ+π)))+4(E^(I (θ+π)))-(E^(-2I (θ+π)))-(E^(2I (θ+π)))}}))

Now we compute the integrals that determine the stiffness matrix. As the support of the basis functions consists of three elements, we obtain a multi-diagonal matrix, where the diagonal entries and the off-diagonal entries have the following values:

Integrate[((D[(Subscript[φ, k,i])[x],x])^2),{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[D[(Subscript[φ, k,i])[x],x]D[(Subscript[φ, k,i+1])[x],x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

Integrate[D[(Subscript[φ, k,i])[x],x]D[(Subscript[φ, k,i+2])[x],x],{x,-∞,∞},Assumptions->(Subscript[h, k])>0]

So, we obtain that the stiffness matrix (Subscript[K, k]) has the following form:

 (Subscript[K, k])=(1/(3(Subscript[h, k])))(({{14, -8,  1,  ,   ,  ,  ,  ,  ,  }, {-8, 16, -8,  0,  ,  ,  ,  ,  ,  }, { 1, -8, 14, -8, 1 ,  ,  ,  ,  ,  }, { , 0 , -8, 16, -8, 0 ,  ,  ,  ,  }, { ,  , 1 , -8, 14, -8, 1 ,  ,  ,  }, { ,  ,  , 0 , -8, 16, -8, 0 ,  ,  }, { ,  ,  ,  , 1 , -8, 14, -8, 1 ,  }, { ,  ,  ,  ,  , 0, -8, 16, -8, 0 }, { ,  ,  ,  ,  ,  , 1 , -8, 14, -8}, { ,  ,  ,  ,  ,  ,  , 0 , -8, 16}}))=(1/(6(Subscript[h, k])))(({{30, -16, 1,  ,   ,  ,  ,  ,  ,  }, {-16, 30, -16, 1 ,  ,  ,  ,  ,  ,  }, { 1, -16, 30, -16, 1 ,  ,  ,  ,  ,  }, { ,  1, -16, 30, -16, 1 ,  ,  ,  ,  }, { ,  ,  1, -16, 30, -16, 1 ,  ,  ,  }, { ,  ,  ,  1, -16, 30, -16, 1,  ,  }, { ,  ,  ,  ,  1, -16, 30, -16, 1 ,  }, { ,  ,  ,  ,  , 1, -16, 30, -16, 1 }, { ,  ,  ,  ,  ,  ,  1, -16, 30, -16}, { ,  ,  ,  ,  ,  ,  , 1 , -16, 30}}))+(1/(6(Subscript[h, k])))(({{-2, 0, 1,  ,   ,  ,  ,  ,  ,  }, {0, 2, 0, -1 ,  ,  ,  ,  ,  ,  }, { 1, 0, -2, 0, 1 ,  ,  ,  ,  ,  }, { ,  -1, 0, 2, 0, -1 ,  ,  ,  ,  }, { ,  ,  1, 0, -2, 0, 1 ,  ,  ,  }, { ,  ,  ,  -1, 0, 2, 0, -1,  ,  }, { ,  ,  ,  ,  1, 0, -2, 0, 1 ,  }, { ,  ,  ,  ,  , -1, 0, 2, 0, -1 }, { ,  ,  ,  ,  ,  ,  1, 0, -2, 0}, { ,  ,  ,  ,  ,  ,  , -1 , 0, 2}}))=:A+B

So, the symbol reads as follows:

(Subscript[K, k_])[θ_]:=(1/(6(Subscript[h, k])))(({{30-16 (E^(-I θ))-16(E^(I θ))+(E^(-2I θ))+(E^(2I θ)), -2+(E^(-2I θ))+(E^(2I θ))}, {-2+(E^(-2I (θ+π)))+(E^(2I (θ+π))), 30-16 (E^(-I (θ+π)))-16(E^(I (θ+π)))+(E^(-2I (θ+π)))+(E^(2I (θ+π)))}}))

The next step is to determine the intergrid-taransfer matrix:

L=10;

(Subscript[ϕ, k_])[θ_,x_]:=(∑_(i=-L)^L)(Subscript[φ, k,i])[x]Exp[i θ I]

Here, we need 4 frequencies to be able to reconstruct a coarse-grid function on the fine grid:

Simplify[(Subscript[ϕ, k-1])[2θ,0] == A0 (Subscript[ϕ, k])[θ,0]+A1 (Subscript[ϕ, k])[θ+π/2,0]+A2 (Subscript[ϕ, k])[θ+2π/2,0]+A3 (Subscript[ϕ, k])[θ+3π/2,0],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

Simplify[(Subscript[ϕ, k-1])[2θ,(Subscript[h, k])/2] == A0 (Subscript[ϕ, k])[θ,(Subscript[h, k])/2]+A1 (Subscript[ϕ, k])[θ+π/2,(Subscript[h, k])/2]+A2 (Subscript[ϕ, k])[θ+2π/2,(Subscript[h, k])/2]+A3 (Subscript[ϕ, k])[θ+3π/2,(Subscript[h, k])/2],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

Simplify[(Subscript[ϕ, k-1])[2θ,(Subscript[h, k])] == A0 (Subscript[ϕ, k])[θ,(Subscript[h, k])]+A1 (Subscript[ϕ, k])[θ+π/2,(Subscript[h, k])]+A2 (Subscript[ϕ, k])[θ+2π/2,(Subscript[h, k])]+A3 (Subscript[ϕ, k])[θ+3π/2,(Subscript[h, k])],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

Simplify[(Subscript[ϕ, k-1])[2θ,3(Subscript[h, k])/2] == A0 (Subscript[ϕ, k])[θ,3(Subscript[h, k])/2]+A1 (Subscript[ϕ, k])[θ+π/2,3(Subscript[h, k])/2]+A2 (Subscript[ϕ, k])[θ+2π/2,3(Subscript[h, k])/2]+A3 (Subscript[ϕ, k])[θ+3π/2,3(Subscript[h, k])/2],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

Solve[{%,%%,%%%,%%%%},{A0,A1,A2,A3}]

{a0,a1,a2,a3}={A0,A1,A2,A3}/.(%[[1]])//Expand

Here, the symbol of the intergrid transfer is not only a 4-dimensional vector. Becuase we already have 2x2-symbols on the coarse grid, the intergrid transfer is a 2x4-matrix, mapping between

	span{(Subscript[(UnderBar[ϕ]), k-1])(2θ),(Subscript[(UnderBar[ϕ]), k-1])(2θ+π)} 	and	span{(Subscript[(UnderBar[ϕ]), k])(θ),(Subscript[(UnderBar[ϕ]), k])(θ+π/2),(Subscript[(UnderBar[ϕ]), k])(θ+π),(Subscript[(UnderBar[ϕ]), k])(θ+3π/2)} 
	
The symbol reads as follows:

{{a0,a1,a2,a3},{a3,a0,a1,a2}/.θ->θ+π/2}//FullSimplify//MatrixForm

P[θ_]:=(({{-((Cos[(θ/2)])^4) (-2+Cos[θ]), (((Sin[θ])^3)/4), (2+Cos[θ]) ((Sin[(θ/2)])^4), -(1/4) ((Sin[θ])^3)}, {-(1/4) ((Cos[θ])^3), (1/4) (((Cos[(θ/2)]-Sin[(θ/2)]))^4) (2+Sin[θ]), (((Cos[θ])^3)/4), -(1/4) (((Cos[(θ/2)]+Sin[(θ/2)]))^4) (-2+Sin[θ])}}))

As outlined in the paper, we have to set up also the symbols for the mass matrix and the stiffnes matrix for the four-dimensional basis

	((Subscript[(UnderBar[ϕ]), k])(θ),(Subscript[(UnderBar[ϕ]), k])(θ+π/2),(Subscript[(UnderBar[ϕ]), k])(θ+π),(Subscript[(UnderBar[ϕ]), k])(θ+3π/2)) 

Here, we define the symbol based on the symbol for the two-dimensional basis ((Subscript[(UnderBar[ϕ]), k])(θ),(Subscript[(UnderBar[ϕ]), k])(θ+π)) as follows:

(Subscript[M2, k_])[θ_]:=(({{(Subscript[M, k])[θ][[1,1]], 0, (Subscript[M, k])[θ][[1,2]], 0}, {0, (Subscript[M, k])[θ+π/2][[1,1]], 0, (Subscript[M, k])[θ+π/2][[1,2]]}, {(Subscript[M, k])[θ][[2,1]], 0, (Subscript[M, k])[θ][[2,2]], 0}, {0, (Subscript[M, k])[θ+π/2][[2,1]], 0, (Subscript[M, k])[θ+π/2][[2,2]]}}))

(Subscript[K2, k_])[θ_]:=(({{(Subscript[K, k])[θ][[1,1]], 0, (Subscript[K, k])[θ][[1,2]], 0}, {0, (Subscript[K, k])[θ+π/2][[1,1]], 0, (Subscript[K, k])[θ+π/2][[1,2]]}, {(Subscript[K, k])[θ][[2,1]], 0, (Subscript[K, k])[θ][[2,2]], 0}, {0, (Subscript[K, k])[θ+π/2][[2,1]], 0, (Subscript[K, k])[θ+π/2][[2,2]]}}))

Now, we check if the Galerkin identity is satisfied (up to scaling):

FullSimplify[P[θ].(Subscript[M2, k])[θ].Transpose[P[θ]]==(1/2) (Subscript[M, k-1])[2θ],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

FullSimplify[P[θ].(Subscript[K2, k])[θ].Transpose[P[θ]]==(1/2) (Subscript[K, k-1])[2θ],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]

The symbol of the whole operator (Subscript[G, k]) =(1/((Subscript[h, k])^2))( I - (Subscript[Π, k]) ) ((Subscript[K, k])^(-1))(Subscript[M, k])  reads as follows:

(Subscript[G, k_])[θ_]=FullSimplify[(IdentityMatrix[4]-FullSimplify[Transpose[P[θ]].FullSimplify[Inverse[(1/2) (Subscript[K, k-1])[2θ]]].P[θ].(Subscript[K2, k])[θ],Assumptions->{(Subscript[h, k])>0,(Subscript[h, k-1])==2(Subscript[h, k])}]).Inverse[(Subscript[K2, k])[θ]].(Subscript[M2, k])[θ]/((Subscript[h, k])^2)]

Eigenvalues[%]

In this case, we do not have to determine the supremum anymore, because the spectral radius takes the value (1/10) for all frequencies θ.