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

For the heuristic but sharp analysis for the two-dimensional case, please see

http://www.risc.jku.at/people/vpillwei/sLFA/Analysis2Dallatonce.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))))//FullSimplify

(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))))//FullSimplify

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);

Compute the smoothing rate

Id=IdentityMatrix[2];

(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

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.

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

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

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

(ev=Simplify[Eigenvalues[NM]])//MatrixForm//Timing

This is a double eigenvalue:

ev[[1]]-ev[[2]]

σ=ev[[1]];

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

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

Now we want to compute the supremum of σ with respect to all η>0  and all high frequencies (-1 ≤ c1 ≤ 0 and  -1 ≤ c2 ≤ 1) 

Doing it directly is possible if we multiply the inequality with the denominator of σ:

{num,den}= {Numerator[σ],Denominator[σ]};

(res=Resolve[ForAll[{η,c1,c2},η>0&&-1≤c1≤0&&-1≤c2≤1,num≤λ*den],{τ,λ},Reals])//Timing

FullSimplify[Sqrt[Piecewise[Table[{res[[i,2,2]],res[[i,1]]},{i,1,Length[res]}]]]]

qSM[τ_]:=(Piecewise[{{1-(τ/4), 0<τ≤(8/7)}, {(1/2) Abs[2-3 τ], True}}])

Plot[qSM[τ],{τ,0,(4/3)},PlotRange->{0,1},PlotStyle->{Thick},AxesLabel->{"τ",Subscript["q","SM"]}]

As explained in the paper, moreover we have to show that the iteration is non-expansive for low frequencies and all τ.

(Resolve[ForAll[{η,c1,c2},η>0&&-1≤c1≤1&&-1≤c2≤1,num≤λ*den],{τ,λ},Reals])//Timing

We see that here λ ≥ 1 is valid for all τ ∈ [0, (8/9)]. This can be done much faster as follows:

Resolve[ForAll[{η,c1,c2},η>0&&-1≤c1≤1&&-1≤c2≤1,num≤den]]//Timing

Smoothing-rate-based Two-grid convergence analysis

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

Since we consider 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

Ansatz for smoother as introduced in (5.1) and (5.2) in the paper:

ℐ=DiagonalMatrix[{1,1,q,q,q,q,q,q}];

Define the matrix for the two-grid step:

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

The next step is the computation of the norm. Therefore, we use the normal matrix:

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

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

As always, we may replace ((h^4)/α) by η:

(𝒩=Simplify[Together[𝒩/.{h^n_:>(α η)^(n/4)}],Assumptions->α>0];)//Timing

The problem depends on the following variables:

Variables[𝒩]

Here, the eigenvalues cannot be computed directly in a reasonable time. Therefore, we consider all denominators in the matrix...

(𝒩=Together[𝒩];)//Timing

denominators= Denominator[𝒩]//Flatten//Factor//Union

... and find a common denominator:

den = 256 (36+36 (c1^2)+9 (c1^4)+36 (c2^2)-126 (c1^2) (c2^2)-72 (c1^4) (c2^2)+9 (c2^4)-72 (c1^2) (c2^4)+144 (c1^4) (c2^4)+η+4 (c1^2) η+4 (c1^4) η+4 (c2^2) η+16 (c1^2) (c2^2) η+16 (c1^4) (c2^2) η+4 (c2^4) η+16 (c1^2) (c2^4) η+16 (c1^4) (c2^4) η);

Note that the denominators are independent of q. Moreover the denominator den is positive:

Resolve[ForAll[{c1,c2,η},0≤c1≤c1≤1 && η>0,den>0],{},Reals]//Timing

We need the maximal eigenvalue of the matrix 𝒩: based on interpolation in q we guess that 

	0,(q^4), epol + (Sqrt[dpol]), epol - (Sqrt[dpol]) 

with epol and dpol defined below and are each eigenvaluesof multiplicity two.

The functions epol and dpol are provided as their own files:

epol = << epol.m;
dpol = << dpol.m;

Variables[epol]
Variables[dpol]

We verify this guess by considering the factorization of the characteristic polynomial obtained by guessing and the determinant definition of the characteristic polynomial and factor the difference of the two:

interpol = ((λ^2)(((λ-(q^4)))^2)((((λ^2)- 2 epol λ + (epol^2)-dpol ))^2));
charpol = Det[λ Id-𝒩];

Factor[interpol-charpol]//Timing

Based on considering extreme cases for c1 and c2, we determine a bound for the maximal eigenvalue 

	λ(c1,c2,η,q) = epol + (Sqrt[dpol]):

FullSimplify[(epol+Sqrt[dpol])/.{c1->0,c2->0},Assumptions->q>0]

FullSimplify[(epol+Sqrt[dpol])/.{c1->1,c2->0},Assumptions->q>0]

FullSimplify[(epol+Sqrt[dpol])/.{c1->0,c2->1},Assumptions->q>0]

FullSimplify[(epol+Sqrt[dpol])/.{c1->1,c2->1},Assumptions->q>0]

Note that λ(1, 1, 0, q) is not defined. Therefore, we use the following limit:

Simplify[(epol+Sqrt[dpol])/.{η->0,c1->1}]
Limit[%,c2->1]
FullSimplify[%,Assumptions->q>0]

We guess that

	(sup_(c1, c2, η))   λ(c1,c2,η, q) =(((max { (((q^2)+3)/4),q(Sqrt[(q^2)+1])}))^2).

We define these two functions

qGUESS1[q_]:=(((q^2)+3)/4);
qGUESS2[q_]:=q(Sqrt[(q^2)+1]);
qGUESS[q_]:=Max[qGUESS1[q],qGUESS2[q]];

Compute the intersection point:

Simplify[Solve[qGUESS1[q]==qGUESS2[q]]]

(Q=(Sqrt[(1/15) (-5+4 (Sqrt[10]))]))//N

The "guessed" function for the smoothing rate looks as follows:

Plot[{qGUESS1[q],qGUESS2[q],1,qGUESS[q]},{q,0,1},PlotStyle->{{Thick,Dashing[0.02],Red},{Thick,Dashing[0.02],Orange},{Thick,Dashing[0.01],Black},{Thick,Blue}},PlotRange->All]

It remains to prove this upper bound: this needs to be split into several smaller steps because of the overall computational complexity

In either case (left or right bound), we need to show that:

	e(q) + (Sqrt[d (q)]) ≤ qGUESS(((q))^2)
	
or, equivalently, 

	 (Sqrt[d (q)]) ≤ ((qGUESS(q))^2)-e(q).

Note that 

 	∀ x,y:   y ≥ 0 ∧ (x^2) ≤ (y^2) ⟹ x ≤ y

Hence we prove the following:

	 qGUESS(((q))^2)-e(q)≥ 0
	 (((((qGUESS(q))^2)-e(q)))^2) - d(q) ≥ 0

The proof of these two items is split into smaller steps according to the general idea described next.

Step 1:
On the uppermost level we run into cases, where we need to prove for all η>0 that
	
	(Subscript[p, 1])(η) = (Subscript[γ, 0]) + (Subscript[γ, 1])  η ≤ 0

For this we show that both coefficients (Subscript[γ, 0]),(Subscript[γ, 1]) are non-positive, which is a sufficient condition.

Step 2:
For showing (non)positivity, (non)negativity of the coefficients (Subscript[γ, i])=(Subscript[γ, i])(c1,c2,q2) we iterate the idea above. These coefficients are polynomials in q2 of degree 2 or 3.

First we derive general conditions on the coefficients  (Subscript[a, k])of this polynomial to be non-negative for q2 in the corresponding range, i.e., with the generic polynomial

	 (Subscript[p, d])(x)=(∑_(k=0)^d)(Subscript[a, k])(x^k)

we determine the quantifier-free formula ψ((Subscript[a, 0]),(Subscript[a, 1]),...,(Subscript[a, d])) equivalent to ∀ 0 < x ≤ Q2: (Subscript[p, d])(x)>0 (or smaller,...).
Then we verify that for our specific coefficients (Subscript[a, k])(c1,c2) these conditions are satisfied (i.e., resolving quantifiers in the next higher level). In the course of this final step, we make  use of the easy to check logical identities

	A => (B ∧ C) ⧦ (A => B) ∧ (A => C)

and

	A => (B ∨ C) ⧦ (A => B) ∨ (A => C)

and stop the computations as soon as we reach True.

define the generic polynomial of degree d=2,3:

p[2,x_]:=a0+a1 x+a2 x^2;
p[3,x_]:=a0+a1 x + a2 x^2 + a3 x^3;

Step I: Proof of the left bound (0<q ≤ Q)

Step I-A:

Step I-A1: (linear expansion in η)

tmp = Together[(((qGUESS1[q])^2)-epol)*den];

cl = CoefficientList[tmp,η]//Factor;
Length[cl]

We find out that here only even powers of q occur. Therefore, we replace (q^2) by q2.

cl=cl/.{q^n_:>q2^(n/2)};

Step I-A2a: (The constant coefficient is a quadratic polynomial in q2 = (q^2))

cl1 = CoefficientList[cl[[1]],q2]//Factor;
Length[cl1]

To simplify computations we derive some simple conditions on the coefficients:

Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[1]]≥0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[3]]<0],{},Reals]//Timing

and determine general conditions on the coefficients using the knowledge we obtained above:

Timing[normalform = Resolve[a0≥0 &&a2≤0&&ForAll[q2,0<q2≤(Q^2),p[2,q2]≥0],{a0,a1,a2},Reals];]

normalform

The normalform is  of the form: A && (B || C)
The truth of A has been verified before

inpC = normalform[[2,2]]/.Thread[{a0,a1,a2}->cl1];

Timing[cadC = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpC],{c1,c2}]]

Hence we have

	True ∧ (B ∨ True) = True.

Step I-A2b: (The linear coefficient is a quadratic polynomial in q2 = (q^2))

cl1 = CoefficientList[cl[[2]],q2]//Factor;
Length[cl1]

We derive simple conditions on the coefficients:

Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[1]]≥0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[2]]>0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[3]]<0],{},Reals]//Timing

and determine general conditions on the coefficients using the knowledge we obtained above:

(normalform = CylindricalDecomposition[a0≥0 && a1>0 &&a2≤0&&Resolve[ForAll[q2,0<q2≤(Q^2),p[2,q2]≥0],{a0,a1,a2},Reals],{a0,a1,a2}];)//Timing

normalform

this is so simple that we do not split further and conclude:

inp = normalform/.Thread[{a0,a1,a2}->cl1];
Timing[cad = CylindricalDecomposition[Implies[0≤c1≤c2<1,inp],{c1,c2}]]

Step I-B:

tmp = Factor[((((((qGUESS1[q])^2)-epol))^2)-dpol)(den^2)];

We divide through (1/16) den ((q^2)-1)<0 and obtain 

	R = R[c1,c2,q2,η]= ((((qGUESS[q])^2)(qGUESS[q] - epol)^2 - dpol) den (16/((q^2)-1))

and need to show that  R≤0. 

For this we proceed again by considering R as linear polynomial in η and we show that the coefficients are non-positive:

	R[c1,c2,q2,η] = r0[c1,c2,q2] + r1[c1,c2,q2] η

R = Factor[tmp/(1/16 den ((q^2)-1))];

Step I-B1: (linear expansion in η)

cl = CoefficientList[R, η];
Length[cl]

We find out that here only even powers of q occur. Therefore, we replace ((q )^2) by q2.

cl=Together[cl]/.{q^n_:>q2^(n/2)};

Step I-B2: (The constant coefficient is a cubic polynomial in q2 = (q^2))

cl1 = CoefficientList[cl[[1]],q2]//Factor;
Length[cl1]

We derive simple conditions on the coefficients:

Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[1]]≤0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[2]]≤0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[3]]≥0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[4]]≥0],{},Reals]//Timing

and determine general conditions on the coefficients using the knowledge we obtained above:

Timing[normalform = Resolve[a0≤0 && a1≤0 && a2≥0 && a3≥0 &&ForAll[q2,0<q2≤(Q^2),p[3,q2]≤0],{a0,a2,a1,a3},Reals];]

normalform

The result is of the form:

	normalform = A &&( B1 || B2 || B3)

we replace the generic coefficients by our particular ones

inpA = normalform[[1]]/.Thread[{a0,a1,a2,a3}->cl1];
inpB1=normalform[[2,1]]/.Thread[{a0,a1,a2,a3}->cl1];
inpB2=normalform[[2,2]]/.Thread[{a0,a1,a2,a3}->cl1];
inpB3=normalform[[2,3]]/.Thread[{a0,a1,a2,a3}->cl1];

Timing[cadA = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpA],{c1,c2}]]

Timing[cadB1 = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpB1],{c1,c2}];]

Timing[cadB2 = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpB2],{c1,c2}];]

Timing[cadB1orB2 = CylindricalDecomposition[cadB1 || cadB2,{c1,c2}]]

Hence:

	normalform = A &&( B1 || B2 || B3) = True && (True || B3) = True

Step I-B2: (The linear coefficient is a cubic polynomial in q2 = (q^2))

cl1 = CoefficientList[cl[[2]],q2]//Factor;
Length[cl1]

We derive simple conditions on the coefficients:

Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[1]]≤0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[2]]≤0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[3]]≥0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[4]]≥0],{},Reals]//Timing

and determine general conditions on the coefficients using the knowledge we obtained above:

Timing[normalform = Resolve[a0≤0 && a1≤0 && a2≥0 && a3≥0 &&ForAll[q2,0<q2≤(Q^2),p[3,q2]≤0],{a0,a2,a1,a3},Reals];]

normalform

The coefficient has the same structure as the previous coefficient (but with different coefficients):

	normalform = A ∧( B1 ∨ B2 ∨ B3)

inpA = normalform[[1]]/.Thread[{a0,a1,a2,a3}->cl1];
inpB1=normalform[[2,1]]/.Thread[{a0,a1,a2,a3}->cl1];
inpB2=normalform[[2,2]]/.Thread[{a0,a1,a2,a3}->cl1];
inpB3=normalform[[2,3]]/.Thread[{a0,a1,a2,a3}->cl1];

Timing[cadA = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpA],{c1,c2}]]

Timing[cadB1 = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpB1],{c1,c2}]]

Timing[cadB2 = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpB2],{c1,c2}]]

Timing[cadB1orB2 = CylindricalDecomposition[cadB1 || cadB2,{c1,c2}]]

Hence:

	normalform = A ∧( B1 ∨ B2 ∨ B3) = True ∧ (True ∨ B3) = True

This concludes the proof for the case 0 < q ≤ Q.

Step II: Proof of the right bound (q > Q)

Step II-A:

Step II-A1: (linear expansion in η)

tmp = Together[(qGUESS2[q]^2-epol)*den];

cl = CoefficientList[tmp, η] // Factor;
Length[cl]

We find out that here only even powers of q occur. Therefore, we replace (q^2) by q2.

cl=cl/.{q^n_:>q2^(n/2)};

Step II-A2a: (The constant coefficient is a quadratic polynomial in q2 = (q^2))

cl1 = CoefficientList[cl[[1]],q2]//Factor;
Length[cl1]

We determine some obvious conditions on the coefficients:

Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[1]]<0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[2]]>0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[3]]>0],{},Reals]//Timing

and determine general conditions on the coefficients using the knowledge we obtained above:

(normalform = Resolve[a0<0 && a1>0 && a2>0 &&ForAll[q2,(Q^2)<q2<1,p[2,q2]≥0],{a0,a1,a2},Reals])//Timing

We have the following structure:

	normalform = A ∧ (B1 ∨B2 ∨ B3)

inpA = normalform[[1]]/.Thread[{a0,a1,a2}->cl1];
inpB1 = normalform[[2,1]]/.Thread[{a0,a1,a2}->cl1];
inpB2 = normalform[[2,2]]/.Thread[{a0,a1,a2}->cl1];
inpB3 = normalform[[2,3]]/.Thread[{a0,a1,a2}->cl1];

Timing[cadA = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpA],{c1,c2}]]

Timing[cadB3 = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpB3],{c1,c2}]]

Hence

	normalform = A ∧ (B1 ∨B2 ∨ B3)= True ∧ (B1 ∨B2 ∨ True)=True

Step II-A2b: (The linear coefficient is a quadratic polynomial in q2 = (q^2))

cl1 = CoefficientList[cl[[2]],q2]//Factor;
Length[cl1]

We determine some obvious conditions on the coefficients:

Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[1]]<0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[2]]>0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[3]]>0],{},Reals]//Timing

and determine general conditions on the coefficients using the knowledge we obtained above:

(normalform = Resolve[a0<0 && a1>0 && a2>0 &&ForAll[q2,(Q^2)<q2<1,p[2,q2]≥0],{a0,a1,a2},Reals])//Timing

We have the following structure:

	normalform = A ∧ (B1 ∨B2 ∨ B3)

inpA = normalform[[1]]/.Thread[{a0,a1,a2}->cl1];
inpB1 = normalform[[2,1]]/.Thread[{a0,a1,a2}->cl1];
inpB2 = normalform[[2,2]]/.Thread[{a0,a1,a2}->cl1];
inpB3 = normalform[[2,3]]/.Thread[{a0,a1,a2}->cl1];

Timing[cadA = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpA],{c1,c2}]]

Timing[cadB3 = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpB3],{c1,c2}]]

Hence

	normalform = A ∧ (B1 ∨B2 ∨ B3)= True ∧ (B1 ∨B2 ∨ True)=True

Step II-B:

tmp = Factor[((((((qGUESS2[q])^2)-epol))^2)-dpol)(den^2)];

Since 
		tmp = - den * (q^2) * S[c1,c2,q2,η] 

we can factor out (- den * (q^2))  and show for the remaining part 

		 S[c1,c2,q2,η]≤0 

proceeding by first considering coefficients of η>0 and then splitting the coefficient where necessary as above

S = Factor[tmp/(-den*(q^2))];

Step II-B1: (linear expansion in η)

cl = CoefficientList[S,η]//Factor;
Length[cl]

We find out that here only even powers of q occur. Therefore, we replace (q^2) by q2.

cl=cl/.{q^n_:>q2^(n/2)};

Step II-B2a: (The constant coefficient is a quadratic polynomial in q2 = (q^2))

cl1 = CoefficientList[cl[[1]],q2]//Factor;
Length[cl1]

We determine some obvious conditions on the coefficients:

Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[1]]>0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[2]]<0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[3]]<0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[4]]≥0],{},Reals]//Timing

and determine general conditions on the coefficients using the knowledge we obtained above:

(normalform=Resolve[a0>0 && a1<0&&a2<0&&a3≥0&&ForAll[x,(Q^2)<x<1,p[3,x]≤0],{a1,a3,a2,a0},Reals])//Timing

We have the following structure:

	normalform = A ∧ (B1 ∨B2 ∨ B3∨ B4)

inpA = normalform[[1]]/.Thread[{a0,a1,a2,a3}->cl1];
inpB1 =normalform[[2,1]]/.Thread[{a0,a1,a2,a3}->cl1];

Timing[cadA = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpA],{c1,c2}]]

Timing[cadB1 = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpB1],{c1,c2}]]

Hence

	normalform = A ∧ (B1 ∨B2 ∨ B3∨ B4)= True ∧ (True ∨B2 ∨ B3∨ B4)=True

Step II-B2b: (The linear coefficient is a cubic polynomial in q2 = (q^2))

cl1 = CoefficientList[cl[[2]],q2]//Factor;
Length[cl1]

We determine some obvious conditions on the coefficients:

Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[1]]>0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[2]]<0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[3]]<0],{},Reals]//Timing
Resolve[ForAll[{c1,c2},0≤c1≤c2<1,cl1[[4]]≥0],{},Reals]//Timing

and determine general conditions on the coefficients using the knowledge we obtained above:

(normalform=Resolve[a0>0 && a1<0&&a2<0&&a3≥0&&ForAll[x,(Q^2)<x<1,p[3,x]≤0],{a1,a3,a2,a0},Reals])//Timing

We have the following structure:

	normalform = A ∧ (B1 ∨B2 ∨ B3∨ B4)

inpA = normalform[[1]]/.Thread[{a0,a1,a2,a3}->cl1];
inpB1 =normalform[[2,1]]/.Thread[{a0,a1,a2,a3}->cl1];

Timing[cadA = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpA],{c1,c2}]]

Timing[cadB1 = CylindricalDecomposition[Implies[0≤c1≤c2<1,inpB1],{c1,c2}]]

Hence

	normalform = A ∧ (B1 ∨B2 ∨ B3∨ B4)= True ∧ (True ∨B2 ∨ B3∨ B4)=True

Interpretation

This finishes the of proof, i.e., we know that qGUESS is correct.

The plot of the Function qGUESS, discribing the relation between the two-grid rate and the smoothing rate, looks as follows.

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

We can plug the smoothing rate qSM into the function qGUESS, to obtain the two-grid rate  as a function of the damping parameter τ. Here, the blue graph represents 1+1 smoothing steps and the purple graph represents 2+2 smoothung steps.

Plot[{qGUESS[qSM[τ]],qGUESS[((qSM[τ])^2)],1},{τ,0,1.4},PlotStyle->{{Thick},{Thick},{Thick,Dashed}},PlotRange->{0,Sqrt[2]},AxesLabel->{"τ",Subscript["q","TG"]}]

Note that q1≤1 could only be shown for τ ∈[0,8/9]. For τ>8/9, we obtain q1>1. Therefore, the graph above is not correct in this case.

Nonetheless, we can also use the proposed approach there. Note that ℐ(q1,q2)=q1 ℐ(1,q2/q1) holds. Hence, we can use the formula derived above also for τ>8/9 and obtain the following result.

qB[τ_]:=Piecewise[{{(1/16) (16-72 τ+81 (τ^2)),τ≤0},{1,0<τ≤(8/9)},{(1/16) (16-72 τ+81 (τ^2)),τ>(8/9)}}]

qTotal[τ_,ν_]:=qGUESS[(((((qSM[τ])/(qB[τ]))))^ν)] ((qB[τ])^(2ν))

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

The optimal value for ν = νpre + νpost = 1 + 1 is:

qTotal[(8/9),1]

which is approximately

N[qTotal[(8/9),1]]

The optimal value for ν = νpre + νpost = 2 + 2 is:

qTotal[(8/9),2]

which is approximately

N[qTotal[(8/9),2]]

list=Transpose[{Table[i/10,{i,1,15}],{1.1136,0.8990,0.8551,0.8061,0.7581,0.7112,0.6789,0.6306,0.5893,0.5567,0.7715,1.0202,1.3055,1.6293,1.9945}}]

Show[{Plot[{qTotal[τ,1],qTotal[τ,2],1},{τ,0,1.4},PlotStyle->{{Thick},{Thick},{Thick,Dashed}},PlotRange->{0,Sqrt[2]},AxesLabel->{"τ",Subscript["q","TG"]}],ListPlot[list,PlotRange->{0,Sqrt[2]},AxesOrigin->{0,0}]}]