An all-at-once symbolic Local Fourier Analysis for 2d distributed elliptic optimal control model problem

This Mathematica notebook accompanies the paper “A local Fourier convergence analysis of a multigrid mehtod using symbolic computation” (V. Pillwein and S. Takacs). A preprint version is available at 

https://www.dk-compmath.jku.at/publications

The notation here follows the notation of this paper. In this notebook the computations for the two-dimensional problem are summarized. For the computations of the one-dimensional case, please see 

http://www.risc.jku.at/people/vpillwei/sLFA/Analysis1D.nb

and for the rigorous, non-sharp analysis for the two-dimensional case, please see

http://www.risc.jku.at/people/vpillwei/sLFA/Analysis2D.nb

As mentioned in the paper, we obtain for the symbols of the stiffness matrix and the mass matrix:

(Subscript[K, k])= (1/3)(8-(E^(-I θ1))-(E^(I θ1))-(E^(-I θ2))-(E^(I θ2))-(E^(I (-θ1-θ2)))-(E^(I (-θ1+θ2)))-(E^(I (θ1-θ2)))-(E^(I (θ1+θ2))));

(Subscript[M, k])=((h^2)/36)(16+4(E^(-I θ1))+4(E^(I θ1))+4(E^(-I θ2))+4(E^(I θ2))+(E^(I (-θ1-θ2)))+(E^(I (-θ1+θ2)))+(E^(I (θ1-θ2)))+(E^(I (θ1+θ2))));

The preconditioning matrices are (for the case that the boundary is neglected) just scaled identity matrices. Therefore the symbol equals the scaling factor.

(Subscript[(OverHat[K]), k])=(8/3);

(Subscript[(OverHat[M]), k])=((16 (h^2))/36);

Id=IdentityMatrix[2];

The smoother: collective Jacobi smoother

(S=FullSimplify[Id-τ Inverse[(({{(Subscript[(OverHat[M]), k]), (Subscript[(OverHat[K]), k])}, {(Subscript[(OverHat[K]), k]), -(α^(-1))(Subscript[(OverHat[M]), k])}}))].(({{(Subscript[M, k]), (Subscript[K, k])}, {(Subscript[K, k]), -(α^(-1))(Subscript[M, k])}}))])//MatrixForm//Timing

(bigS = (({{S[[1,1]], S[[1,2]], 0, 0, 0, 0, 0, 0}, {S[[2,1]], S[[2,2]], 0, 0, 0, 0, 0, 0}, {0, 0, S[[1,1]]/.{θ1->θ1+π}, S[[1,2]]/.{θ1->θ1+π}, 0, 0, 0, 0}, {0, 0, S[[2,1]]/.{θ1->θ1+π}, S[[2,2]]/.{θ1->θ1+π}, 0, 0, 0, 0}, {0, 0, 0, 0, S[[1,1]]/.{θ2->θ2+π}, S[[1,2]]/.{θ2->θ2+π}, 0, 0}, {0, 0, 0, 0, S[[2,1]]/.{θ2->θ2+π}, S[[2,2]]/.{θ2->θ2+π}, 0, 0}, {0, 0, 0, 0, 0, 0, S[[1,1]]/.{θ1->θ1+π,θ2->θ2+π}, S[[1,2]]/.{θ1->θ1+π,θ2->θ2+π}}, {0, 0, 0, 0, 0, 0, S[[2,1]]/.{θ1->θ1+π,θ2->θ2+π}, S[[2,2]]/.{θ1->θ1+π,θ2->θ2+π}}}))) // MatrixForm

We apply the transformation (e^(i (Subscript[θ, j])))-> cos (Subscript[θ, j])+i sin (Subscript[θ, j]) =: (Subscript[c, j]) + i (Subscript[s, j]) for j=1,2.

(bigS=bigS/.{Cos[θ1]->c1,Cos[θ2]->c2})//MatrixForm

The coarse-grid correction

First introduce the intergrid transfer (projection) operator P=(Subscript[P, k]^(k-1)). Here we have a coupling between frequency ((Subscript[θ, 1]),(Subscript[θ, 2])) and ((Subscript[θ, 1])+π,(Subscript[θ, 2])), ((Subscript[θ, 1]),(Subscript[θ, 2])+π) and ((Subscript[θ, 1])+π,(Subscript[θ, 2])+π):

P=(({{1+(1/2)(E^(-I θ1))+(1/2)(E^(I θ1))+(1/2)(E^(-I θ2))+(1/2)(E^(I θ2))+(1/4)(E^(I (-θ1-θ2)))+(1/4)(E^(I (-θ1+θ2)))+(1/4)(E^(I (θ1-θ2)))+(1/4)(E^(I (θ1+θ2)))}, {1+(1/2)(E^(-I (θ1+π)))+(1/2)(E^(I (θ1+π)))+(1/2)(E^(-I θ2))+(1/2)(E^(I θ2))+(1/4)(E^(I (-(θ1+π)-θ2)))+(1/4)(E^(I (-(θ1+π)+θ2)))+(1/4)(E^(I ((θ1+π)-θ2)))+(1/4)(E^(I ((θ1+π)+θ2)))}, {1+(1/2)(E^(-I θ1))+(1/2)(E^(I θ1))+(1/2)(E^(-I (θ2+π)))+(1/2)(E^(I (θ2+π)))+(1/4)(E^(I (-θ1-(θ2+π))))+(1/4)(E^(I (-θ1+(θ2+π))))+(1/4)(E^(I (θ1-(θ2+π))))+(1/4)(E^(I (θ1+(θ2+π))))}, {1+(1/2)(E^(-I (θ1+π)))+(1/2)(E^(I (θ1+π)))+(1/2)(E^(-I (θ2+π)))+(1/2)(E^(I (θ2+π)))+(1/4)(E^(I (-(θ1+π)-(θ2+π))))+(1/4)(E^(I (-(θ1+π)+(θ2+π))))+(1/4)(E^(I ((θ1+π)-(θ2+π))))+(1/4)(E^(I ((θ1+π)+(θ2+π))))}}));

Now we may introduce the symbols for the involved system matrices on the fine grid. Here we have to take the block-structure of the system matrix 𝒜 and the preconditioner (OverHat[𝒜]) on the one hand and the fact that we want to consider both,  ((Subscript[θ, 1]),(Subscript[θ, 2])) and ((Subscript[θ, 1])+π,(Subscript[θ, 2])), ((Subscript[θ, 1]),(Subscript[θ, 2])+π) and ((Subscript[θ, 1])+π,(Subscript[θ, 2])+π), at the same time, on the other hand:

((Subscript[A, k])=(({{(Subscript[M, k]), (Subscript[K, k]), 0, 0, 0, 0, 0, 0}, {(Subscript[K, k]), -(α^(-1))  (Subscript[M, k]), 0, 0, 0, 0, 0, 0}, {0, 0, (Subscript[M, k])/.{θ1->θ1+π}, (Subscript[K, k])/.{θ1->θ1+π}, 0, 0, 0, 0}, {0, 0, (Subscript[K, k])/.{θ1->θ1+π}, -(α^(-1)) (Subscript[M, k])/.{θ1->θ1+π}, 0, 0, 0, 0}, {0, 0, 0, 0, (Subscript[M, k])/.{θ2->θ2+π}, (Subscript[K, k])/.{θ2->θ2+π}, 0, 0}, {0, 0, 0, 0, (Subscript[K, k])/.{θ2->θ2+π}, -(α^(-1)) (Subscript[M, k])/.{θ2->θ2+π}, 0, 0}, {0, 0, 0, 0, 0, 0, (Subscript[M, k])/.{θ1->θ1+π,θ2->θ2+π}, (Subscript[K, k])/.{θ1->θ1+π,θ2->θ2+π}}, {0, 0, 0, 0, 0, 0, (Subscript[K, k])/.{θ1->θ1+π,θ2->θ2+π}, -(α^(-1)) (Subscript[M, k])/.{θ1->θ1+π,θ2->θ2+π}}})))//MatrixForm

And for the intergrid transfer operator 𝒫=(Subscript[𝒫, k]^(k-1)) for two variables reads as:

(𝒫=(({{P[[1]], 0, P[[2]], 0, P[[3]], 0, P[[4]], 0}, {0, P[[1]], 0, P[[2]], 0, P[[3]], 0, P[[4]]}})))//MatrixForm

We apply the transformation (e^(i θ))-> cos θ+i sin θ =: c + i s

((Subscript[A, k])=(Re[(Subscript[A, k])]+I Im[(Subscript[A, k])]//FullSimplify)/.{Sin[θ1]->s1,Cos[θ1]->c1,Sin[θ2]->s2,Cos[θ2]->c2})//MatrixForm

Rewrite the exponential function in terms of sine and cosine. This yields an expression not depending on sine in which we replace Cos[θ1], Cos[θ2] by new variables c1 and c2:

(𝒫 =Simplify[𝒫/.Exp[x_]->Cos[I x]+I Sin[I x]]/.{Sin[θ1]->s1,Cos[θ1]->c1,Sin[θ2]->s2,Cos[θ2]->c2})//MatrixForm

Due to the fact that we have a conforming discretization, we have:

((Subscript[A, k-1])=Simplify[𝒫.(Subscript[A, k]).Transpose[𝒫]])//MatrixForm

Define the identity matrix:

Id=IdentityMatrix[8];

Derive the symbol for the coarse-grid correction operator:

CG=Id-Transpose[𝒫].Inverse[(Subscript[A, k-1])].𝒫.(Subscript[A, k]);

(CG=FullSimplify[CG,Assumptions->{h>0,α>0,0≤c1≤1,0≤c2≤1}];)//Timing

TG = bigS.CG.bigS//Factor;

Now try our favorite η-substitution:

TG=TG/.{α->h^4/η};

ℒ=DiagonalMatrix[{1,α^-1,1,α^-1,1,α^-1,1,α^-1}];

DiagonalMatrix[{1,(h^2) ,1,(h^2),1,(h^2) ,1,(h^2) }]//MatrixForm

Determine the normal matrix:
Instead of

	(ℒ^(-1/2))(TG^T)ℒ  TG (ℒ^(-1/2))

we use

	(ℵ^(-1)) (ℒ^(-1))(TG^T)ℒ  TG  ℵ,

where  ℵ = (({{1, 0, 0, 0, 0, 0, 0, 0}, {0, (h^2), 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, (h^2), 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, (h^2), 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 0, 0, (h^2)}})) as both have the same eigenvalues, but the latter matrix already scales more nicely in the matrix entries.

First we compute the left factor and already perform the η-substitution:

ℵ=DiagonalMatrix[{1,(h^2) ,1,(h^2),1,(h^2) ,1,(h^2) }];
LeftFac = Inverse[ℵ].Inverse[ℒ]/.{α->(h^4)/η};

NM = LeftFac.Transpose[TG].ℒ.TG.ℵ;

Now try our favorite η-substitution:

NM=NM/.{α->(h^4)/η};

Timing[NM = Together[NM];]

Variables[NM]

The substitution results in a matrix that does not depend on  h or α anymore.

The next step is to compute the eigenvalues for boundary values of (c1,c2,η):

A) Set (c1,c2) = (1,1) and compute the eigenvalues of NM for symbolic τ and η,

tmp = NM/.{c1->1,c2->1}//Factor;

evs = Eigenvalues[tmp]//Union//Factor

Note that the supremum over all η>0 of these eigenvalues can easily be computed by hand:

	(sup_(η>0))max{0,((((144(((2-3τ))^2)+η (((4-3τ))^2)))^2)/(256(((η+36))^2))),((((((η(τ-4))^2)+576 (((τ-1))^2)))^2)/(256(((η+36))^2)))}
		=(sup_(η>0))max(({0,((144(((2-3τ))^2)+η (((4-3τ))^2))/(16(η+36))),((((η(τ-4))^2)+576 (((τ-1))^2))/(16(η+36)))})^2)
		=max(({0,((144(((2-3τ))^2))/(16*36)),((576 (((τ-1))^2))/(16*36)),((((4-3τ))^2)/16),((((τ-4))^2)/16)})^2)
i.e, as:

bd11eta= Max[evs/.η->0,Limit[evs,η->Infinity]]//Factor

res11eta=FullSimplify[Sqrt[bd11eta],Assumptions->τ>0]

in other words,

res11eta = Max[(1/16) (((-4+τ))^2),(1/4) (((2-3 τ))^2)];

Plot[res11eta,{τ,0,1.5},PlotStyle->{Black,Thick,Dotted},PlotRange->{0,Sqrt[2]},AxesLabel->{"τ",Subscript["q","TG"]}]

B) As a second boundary value, we consider first the point (c2,η) = (1,0) and then compute the supremum over all c1:

tmp = NM/.{c2->1,η->0}//Factor;

evs = Eigenvalues[tmp]//Union//Factor

The eigenvalue 0 does not need to be considered, and we compute the maximum of the remaining 3 eigenvalues using CAD:

cadc110=Factor/@Resolve[0<τ&& ForAll[c1,0≤c1≤1,evs[[2]]≤B && evs[[3]]≤B && evs[[4]]≤B],{τ,B},Reals]

From this result we extract the following bound:

bdc110[τ_]:=Piecewise[{{(1/16) (((-2+3 τ))^2) (8-12 τ+9 (τ^2)),0≤τ≤Root[256-1536 #1+3168 ((#1)^2)-3024 ((#1)^3)+1215 ((#1)^4)&,1]},{(1/256) (((-4+3 τ))^4),Root[256-1536 #1+3168 ((#1)^2)-3024 ((#1)^3)+1215 ((#1)^4)&,1]<τ≤Root[256-1536 #1+3168 ((#1)^2)-3024 ((#1)^3)+1215 ((#1)^4)&,2]},{(1/16) (((-2+3 τ))^2) (8-12 τ+9 (τ^2)),Root[256-1536 #1+3168 ((#1)^2)-3024 ((#1)^3)+1215 ((#1)^4)&,2]<τ}}]

which yields after taking the square root rewritten in terms of a maximum the bound:

resc110 = Max[(1/16) (((4-3 τ))^2),(1/4) (Sqrt[8+3 τ (-4+3 τ)]) Abs[2-3 τ]];

Plot[{res11eta,resc110},{τ,0,1.5},PlotStyle->{{Black,Thick,Dotted},{Black,Thick,Dashed}},PlotRange->{0,Sqrt[2]},AxesLabel->{"τ",Subscript["q","TG"]}]

The remaining parts of the boundary of (([0, 1])^2)×(0,∞) do not yield any further information concerning the supremum, hence the overall maximum is given by Max[res11eta,resc110], i.e., by

bdTG[τ_]:=Max[(1/16) (((-4+τ))^2),(1/4) (Sqrt[8+3 τ (-4+3 τ)]) Abs[2-3 τ]];

pl=Plot[{res11eta,resc110,bdTG[τ]},{τ,0,1.5},PlotStyle->{{Black,Dotted},{Black,Dashed},{Gray,Thick}},PlotRange->{0,Sqrt[2]}]

Now we compare the obtained bound with the values computed using gridding in Matlab:

gridding={1.4008,1.1067,0.9025,0.8556,0.8100,0.7656,0.7225,0.6806,0.6400,0.6006,0.5625,0.7358,0.9744,1.2484,1.5602,1.9115}

Show[{pl,ListPlot[Transpose[{Table[i/10,{i,0,15}],gridding}],PlotStyle->{Black,Thick}]}]