symbolic Local Fourier Analysis for 1d 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 one-dimensional problem are summarized. For the computations of the two-dimensional case, please see 

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

and

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

As mentioned in the paper, the construction of stiffness and mass matrix for the one dimensional case can be done similarly to the case of the two-dimensional case. We obtain for the symbols of the stiffness matrix and the mass matrix:

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

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

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])=(2/h);

(Subscript[(OverHat[M]), k])=((2 h)/3);

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

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

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, θ and θ+π, at the same time, on the other hand:

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

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

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

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

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

((Subscript[A, k])=(Subscript[A, k])/.{(E^(I θ))->c+I s, (E^(-I θ))->c-I s,(E^(I (θ+π)))->-(c+I s), (E^(-I (θ+π)))->-(c-I s)})//MatrixForm

(𝒫=Simplify[𝒫/.{(E^(I θ))->c+I s, (E^(-I θ))->c-I s,(E^(I (θ+π)))->-(c+I s), (E^(-I (θ+π)))->-(c-I s)}])//MatrixForm

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

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

Define the identity matrix:

Id=DiagonalMatrix[{1,1,1,1}];

Derive the symbol for the smoother:

(S=Id-τ Inverse[(Subscript[(OverHat[A]), k])].(Subscript[A, k])//FullSimplify)//MatrixForm//Timing

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}])//MatrixForm//Timing

Derive the symbol for the two-grid-operator for 1 pre- and 1 post-smoothing step:

(TG=FullSimplify[S.CG.S,Assumptions->{h>0,0≤c≤1,α>0}];)//Timing

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

(NM=Transpose[TG].ℒ.TG.Inverse[ℒ]//Simplify;)//Timing

(ev=Eigenvalues[NM])//Timing

We have two eigenvalues (0 and (Subscript[λ, max])) with multiplicity 2 each:

ev[[3]]-ev[[4]]//FullSimplify//Timing

((σ =ev[[3]])//FullSimplify)//Timing

Do the replacement (h^4)/α -> η:

(σ =(σ/.{h->(((η α))^(1/4))})//FullSimplify)//Timing

D[σ,α]//FullSimplify

(σ =(σ/.{α->1})//FullSimplify)//Timing

We see that in this formula only even powers of c exist. Therefore we replace (c^2)-> γ

σ=σ/.{(c^(n_)):>(γ^(n/2))/;EvenQ[n]}

Now we want to compute the supremum of σ with respect to all η>0  and all 0≤γ≤1. Doing it directly is not possible in acceptable time:

(*CylindricalDecomposition[ForAll[{η,γ},η>0&&0≤γ≤1,σ≤λ],{τ,λ}]*)

We sample for η->0 and η->∞:

σ0=Limit[σ,η->0]

(res0 =Resolve[ τ>0 && ForAll[γ,0≤γ≤1,σ0≤λ],{τ,λ},Reals])//Timing

σinf=Limit[σ,η->∞]

(resinf = Resolve[τ>0 && ForAll[γ,0≤γ≤1,σinf≤λ],{τ,λ},Reals])//Timing

Compute the maximum of both:

res=CylindricalDecomposition[res0&&resinf,{τ,λ}]

We now guess that the maximum of (Subscript[q, 0]) and (Subscript[q, ∞]) is q. 

For every fixed τ we can find the smallest λ satisfying the logical expression above. We obtain the following piecewise defined function:

q0[τ_]:=2-12 τ+28 (τ^2)-32 (τ^3)+16 (τ^4);
q1[τ_]:=(1/16) (16-32 τ+24 (τ^2)-8 (τ^3)+(τ^4));
{τ1,τ2}={Root[16-160 #1+424 ((#1)^2)-504 ((#1)^3)+255 ((#1)^4)&,1],Root[16-160 #1+424 ((#1)^2)-504 ((#1)^3)+255 ((#1)^4)&,2]};

q[τ_]:=Piecewise[{{Sqrt[q0[τ]],τ≤τ1},{Sqrt[q1[τ]],τ1<τ≤τ2},{Sqrt[q0[τ]],τ2<τ}}]

Still we have to proof that the guess is really an upper bound. (If we can show this, we are by construction sure that it is an sharp bound.) We may again use CAD to do the proof:

CylindricalDecomposition[ForAll[{τ,γ,η},0≤γ≤1&&η>0&&0<τ≤τ1,σ≤q0[τ]],{}]//Timing

CylindricalDecomposition[ForAll[{τ,γ,η},0≤γ≤1&&η>0&&τ1<τ≤τ2,σ≤q1[τ]],{}]//Timing

(*CylindricalDecomposition[ForAll[{τ,γ,η},0≤γ≤1&&η>0&&τ2<τ,σ≤q0[τ]],{}]//Timing*)

This is too slow.

Now consider smaller parts:

CylindricalDecomposition[ForAll[{τ,γ,η},0≤γ≤1&&η>0&&τ2<τ≤ 7/10,σ≤q0[τ]],{}]//Timing

CylindricalDecomposition[ForAll[{τ,γ,η},0≤γ≤1&&η>0&&7/10≤ τ≤1,σ≤q0[τ]],{}]//Timing

Thus we have proven the identity for τ∈[0,1].

We may plot the solution:

Plot[{q[τ],1},{τ,0,1},PlotRange->{0,Sqrt[2]},PlotStyle->{Thick,{Thick,Dashed}},AxesLabel->{"τ",Subscript["q","TG"]}]

Alternatively, we may split the proof of the last part following the idea exploited in the 2D case:

num = Numerator[Together[σ]];
den = Denominator[Together[σ]];

Resolve[ForAll[{η,γ,τ},η>0 && 0≤γ≤1 && τ>0,den>0],{},Reals]

expr = Collect[q0[τ]*den-num,η,Factor]

Hence this is a polynomial of degree 3 in η. Since η>0 it is sufficient to verify that all coefficients are nonnegative to prove the last part of the bound:

cl = CoefficientList[expr,η];

Resolve[ForAll[{γ,τ},0≤γ≤1 && τ>τ2,cl[[1]]≥0],{},Reals]

Resolve[ForAll[{γ,τ},0≤γ≤1 && τ>τ2,cl[[2]]≥0],{},Reals]

Resolve[ForAll[{γ,τ},0≤γ≤1 && τ>τ2,cl[[3]]≥0],{},Reals]

Resolve[ForAll[{γ,τ},0≤γ≤1 && τ>τ2,cl[[4]]≥0],{},Reals]

This finishes the proof.