• Keine Ergebnisse gefunden

control of the


Academic year: 2022

Aktie "control of the"

Mehr anzeigen ( Seite)





Robust preconditioning and error estimates for optimal

control of the

convection-diffusion-reaction equation with limited

observation in Isogeometric Analysis

K.-A. Mardal, J. Sogn, S. Takacs

RICAM-Report 2020-46





Abstract. In this paper we analyze an optimization problem with limited observation governed by a convection–diffusion–reaction equation. Motivated by a Schur complement approach, we arrive at continuous norms that enable analysis of well-posedness and subsequent derivation of error analysis and a preconditioner that is robust with respect to the parameters of the problem. We provide conditions for inf-sup stable discretizations and present one such discretization for box domains with constant convection. We also provide a priori error estimates for this discretization. The preconditioner requires a fourth order problem to be solved. For this reason, we use Isogeometric Analysis as a method of discretization. To efficiently realize the preconditioner, we consider geometric multigrid with a standard Gauss-Seidel smoother as well as a new macro Gauss-Seidel smoother. The latter smoother provides good results with respect to both the geometry mapping and the polynomial degree.

Key words. PDE-constrained optimization, optimal control, robust preconditioning, error estimates

AMS subject classifications. 49K20, 65F08, 65N22, 65N15

1. Introduction. In this paper, we consider an optimal control problem involv- ing a linear Convection–Diffusion–Reaction (CDR) problem, which reads as follows:

(1.1) Minimize J(u, q) := 1


2kqk2L2(Ω) for u∈U, q ∈L2(Ω) subject to

ưε∆u+β· ∇u+σu=fưq in Ω, u= 0 on ∂Ω.


Here and in what follows, Ω is a bounded open subset ofRd(d= 1,2,3) with Lipschitz boundary, f ∈ L2(Ω),ε, σ ∈ Rwith ε >0, σ ≥0, β ∈ L(Ω)d with∇ ·β = 0 and O ⊆ Ω is measurable in Rd. For certain choices of the parameters, like ε β, σ, convection–diffusion–reaction problems are singular perturbation problems exhibiting sharp gradients and a potential for loss of regularity. To overcome the problems associated to the loss of regularity, significant effort has been put in the development of methods with low regularity, such as discontinuous Galerkin methods [2, 9]. We take the opposite approach and investigate to what extent higher regularity may be used in the setting of optimal control problems.

Funding: The first author acknowledges support from the Research Council of Norway, grant 300305 and 301013. The second and the third author are supported by the Austrian Science Fund (FWF): P31048.

Department of Mathematics, University of Oslo, Oslo, Norway. Center for Biomedical Comput- ing, Simula Research Laboratory, Lysaker, Norway ([email protected]).

Johann Radon Institute for Computational Mathematics (RICAM), Austrian Academy of Sci- ences, Linz, Austria ([email protected]).

§Johann Radon Institute for Computational Mathematics (RICAM), Austrian Academy of Sci- ences, Linz, Austria ([email protected]).



There are two main problems with (1.1)–(1.2), namely: 1) potential sharp gradi- ents leading to non-physical oscillations in the numerical solution and 2) ill-posedness due to limited observations, this is, whenOis a subset of Ω. Motivated by the fact that higher regularity has been exploited in the cases with limited observations [17,25,4], we derive order optimal preconditioners via stability analysis in non-standard Sobolev spaces.

When solving the CDR problem, it is common to consider some stabilization method (e.g. the streamline upwind Petrov Galerkin (SUPG) method) or adaptive grids (e.g. Shishkin grids) to reduce the oscillatory behavior, cf. [10]. This is also the case in optimal control settings, see, e.g., [21, 3, 13,8]. We do not use any such stabilization techniques, but we remark that the trial and test functions involved in the state equation differ. That means, the state equation, if considered isolated, is discretized by a Petrov-Galerkin method, although the complete optimality system is discretized by a standard Galerkin method, this is, trial and test functions agree. In particular, in the continuous setting, the trial space and test spaces areH2(Ω)∩H01(Ω) andL2(Ω), respectively, with properly weighted norms.

By considering a Schur complement of the optimal control problem, we derive non- standard norms in which well-posedness is obtained, assuming extra regularity. From the well-posedness of the continuous system we subsequently analyse corresponding discrete systems to arrive at both error estimates and preconditioners that are robust with respect to the problem parameters α, ε, β and σ. In detail, we provide a con- dition for the discretization which ensures that the preconditioner is sparse and that the preconditioned system is stable. Further, we give an example of such a discretiza- tion based on Isogeometric Analysis (IgA) [14, 5]. For our approach, IgA provides useful discretization methods since the extra regularity leads H2(Ω)–conforming ap- proximation spaces. Using these discretization methods, a priori error estimates are derived, where we detail the dependencies of the problem parameters. We remark that the error estimates required extending some approximation error estimates for tensor-product B-splines, which is done in AppendixB.

Similar Schur complement preconditioners were used, on the linear algebra level, in [19,20] for optimal control problems of the CDR equation. [20] also considers mixed constraints. The preconditioners perform well for different values of the problem pa- rameters. The Schur complement preconditioners were replaced with approximations based on the factorization approach by [19]. However, this approach does not work well for problems with limited observation, i.e., whenO(Ω.

To solve the resulting linear system we use preconditioned Krylov subspace meth- ods. We consider two approaches to realize our preconditioner: sparse direct meth- ods and multigrid methods. For mid-sized problems, sparse direct solvers work well since each component of the preconditioner is symmetric and positive definite. For large-sized problems, we use a multigrid method to realize the fourth-order operator.

Combining the results from this paper and from [23, 24], it follows that the multi- grid method we consider is robust in the grid-size, however, not necessarily in any of the other problem parameters. Finding a multigrid method which is robust in the grid-size, the chosen spline degrees,α, ε, β, σ andO, remains an open problem.

The outline of the paper is as follows: in the next section we perform the analysis of the continuous problem. In Section3, we analyse the discrete problem and provide a condition for a stable discretization. IgA is then introduced in Section 4 along with the proposed discretization. In Section5, a priori error estimates are derived.

Section6 contains a discussion of the solution of an one-dimensional problem and in Section 7 we perform numerical experiments on two-dimensional problems and see


how the preconditioner behaves.

2. Analysis of the continuous problem. To obtain a standard (weak) vari- ational formulation of the state equation (1.2), one would choose the state variable u and test function ˜w to be in H01(Ω). Instead, we consider the strong variational formulation, this is, findu∈H2(Ω)∩H01(Ω) such that

(q,w)˜ L2(Ω)+ (ưε∆u+β· ∇u+σu,w)˜ L2(Ω)= (f,w)˜ L2(Ω) ∀w˜∈L2(Ω).

The Lagrangian functional associated to (1.1)–(1.2) is L(q, w, u) := 1



+ (q, w)L2(Ω)+ (ưε∆u+β· ∇u+σu, w)L2(Ω)ư(f, w)L2(Ω), where u ∈ H2(Ω)∩H01(Ω), q ∈ L2(Ω) and the Lagrangian multiplier w ∈ L2(Ω).

From the first order necessary optimality conditions


∂q(q, w, u) = 0, ∂L

∂w(q, w, u) = 0, ∂L

∂u(q, w, u) = 0, which are also sufficient here, we obtain the optimality system:

Problem 2.1. Find (q, w, u)∈L2(Ω)×L2(Ω)×H2(Ω)∩H01(Ω) such that α(q,q)˜L2(Ω)+ (w,q)˜L2(Ω)= 0 ∀q˜∈L2(Ω), (q,w)˜ L2(Ω)+ (ưε∆u+β· ∇u+σu,w)˜ L2(Ω)= (f,w)˜ L2(Ω) ∀w˜∈L2(Ω),

(w,ưε∆˜u+β· ∇˜u+σ˜u)L2(Ω)+ (u,u)˜ L2(O)= (ud,u)˜ L2(O) ∀˜u∈H2(Ω)∩H01(Ω).

Problem2.1can be written as

(2.1) A

 q w u


 0 M f M˜Oud

 with A:=

αM M 0

M 0 K

0 K0 MO


Here,M :L2(Ω)→(L2(Ω))0 represents theL2(Ω)-inner product, that is, we have hM q, wi= (q, w)L2(Ω),

whereh·,·idenotes the duality product. The notation ”0” is used to denote both dual spaces and dual operators. K:H2(Ω)∩H01(Ω)→(L2(Ω))0 is the state operator:

hKu, wi= (ưε∆u+β· ∇u+σu, w)L2(Ω).

Finally, MO : H2(Ω)∩H01(Ω) → (H2(Ω)∩H01(Ω))0 and ˜MO : L2(O) → (H2(Ω)∩ H01(Ω))0, and both represent theL2(O)-inner product on the subdomainO.

We observe that the block operatorA has a block tridiagonal form. Such tridi- agonal operators are studied in [25,4]. We use the Schur complement preconditioner proposed in [25]:

(2.2) S(A) :=

Sq 0 0 0 Sw 0

0 0 Su



where the components are

(2.3) Sq:=αM, Sw:= 1

αM, Su:=MO+αK0M−1K.

These Schur complements define weighted norms as follows:

kqk2Sq :=hSqq, qi=αkqk2L2(Ω), kwk2S

w :=hSww, wi= 1


kuk2Su :=hSuu, ui=kuk2L2(O)+αk −ε∆u+β· ∇u+σuk2L2(Ω). (2.4)

The last norm follows from K0M−1Ku, u

= sup


hKu, wi2

hM w, wi = sup


(−ε∆u+β· ∇u+σu, w)2L2(Ω)


= k −ε∆u+β· ∇u+σuk4L2(Ω)

k −ε∆u+β· ∇u+σuk2L2(Ω) =k −ε∆u+β· ∇u+σuk2L2(Ω). We show well-posedness with respect to the norms (2.4) by showing that the operator (2.5) A:L2(Ω)×L2(Ω)×H2(Ω)∩H01(Ω)→L2(Ω)0×L2(Ω)0×(H2(Ω)∩H01(Ω))0 is an isomorphism with respect to the norms (2.4). This is done by using the main result in [25], which for our problem reads as follows.

Theorem 2.2. Assume the Schur complements in (2.3) are well defined and pos- itive definite, that is,

(2.6) hSqq, qi ≥σqkqk2L2(Ω), hSww, wi ≥σwkwk2L2(Ω), hSuu, ui ≥σukuk2H2(Ω)

for some positive constants σq, σw and σu, which can depend on the given param- eters. Then, A in (2.5) is an isomorphism, moreover, the condition number of the preconditioned operatorS−1A is bounded:

(2.7) κ S−1A

≤ cos(π/7)

sin(π/14) ≈4.05.

The conditions in (2.6) ensure that the spaces L2(Ω), L2(Ω) and H2(Ω)∩H01(Ω), equipped with normsk·kSq,k·kSwandk·kSu, are complete. Before proving Condition (2.6), we provide a useful lemma which bounds theH2-norm. The proof of this lemma is presented in AppendixA.

Lemma 2.3. If the domainΩ has a Lipschitz boundary and

• the boundary is a polygon (polyhedron) or

• the domain is the image of a geometric mappingG:Ω := (0,b 1)d→Ω, where both k∇rGkL(bΩ) andk(∇rG)−1kL(bΩ) are bounded forr∈ {1,2,3}, then theH2-norm is bounded by the L2-norm of the Laplacian, i.e.,

(2.8) kukH2(Ω)≤Ck∆ukL2(Ω) ∀u∈H2(Ω)∩H01(Ω), for a constantC depending only onΩ.


Theorem 2.4. IfΩis a domain such that the conditions of Lemma2.3hold, then the assumptions of Theorem2.2are satisfied for Problem 2.1.

Proof. For simplicity, we prove this lemma only for σ= 0. An extension to the caseσ >0 is straight-forward.

The first two conditions are trivial sincehSqq, qi=αkqk2L2(Ω) and hSww, wi= α1kwk2L2(Ω). For the third condition, let

δ:= kβk


cP +kβk,

where kβk=kβkL(Ω) andcP is constant from the Poincar´e inequality, that is, we have


Note thatδ∈[0,1). Letu∈H2(Ω)∩H01(Ω) be arbitrary but fixed. We consider the two cases:

(2.9) kβkk∇ukL2(Ω)< δεk∆ukL2(Ω) and kβkk∇ukL2(Ω)≥δεk∆ukL2(Ω). First case. From the definition, we have

K0M−1Ku, u

= sup


(β· ∇u−ε∆u, w)2L2(Ω)



By settingw=−ε∆u, we get using the Cauchy–Schwarz inequality that K0M−1Ku, u

≥((β· ∇u,−ε∆u)L2(Ω)2k∆uk2L2(Ω))2 ε2k∆uk2L2(Ω)

≥(−kβkk∇ukL2(Ω)εk∆ukL2(Ω)2k∆uk2L2(Ω))2 ε2k∆uk2L2(Ω)

= (−kβkk∇ukL2(Ω)+εk∆ukL2(Ω))2. Using the first inequality in (2.9), we obtain

K0M−1Ku, u


= ε4

(ε+cPkβk)2k∆uk2L2(Ω). Second case. By settingw=u, we get

K0M−1Ku, u

≥((β· ∇u, u)L2(Ω)+εk∇uk2L2(Ω))2

kuk2L2(Ω) = ε2k∇uk4L2(Ω)


using integration by parts. Due to the homogeneous Dirichlet boundary conditions and∇ ·β = 0, the term (β· ∇u, w)L2(Ω) is skew symmetric and vanishes forw=u.

Finally, we use kukL2(Ω) ≤ cpk∇ukL2(Ω) and the second inequality in (2.9), which gives

K0M−1Ku, u


kuk2L2(Ω) ≥ ε2


≥ δ2ε4

kβk2c2Pk∆uk2L2(Ω)= ε4



To summarize, in both cases we get hSuu, ui=kuk2L2(O)

K0M−1Ku, u


K0M−1Ku, u

≥α ε4

(ε+cPkβk)2k∆uk2L2(Ω)≥α Cε4

(ε+cPkβk)2kuk2H2(Ω). The last inequality follows from Lemma2.3.

Theorem2.2and Theorem2.4show that Problem2.1is well-posed with respect to the norms in (2.4). The boundedness and coercivity constants are bounded indepen- dent from the regularization parameterαas well as the problem parametersε,βand σ. Consequently, the operator preconditioner (2.2) is a robust preconditioner for the optimality system, that is, the condition number is uniformly bounded independently of the above mentioned parameters. So far, we have only analyzed the problem on the continuous level. In the next section, we carry this analysis over to the discrete case and provide a computationally feasible preconditioner.

3. Analysis of the discrete problem. We consider conforming discretizations, that is, we choose the finite-dimensional spacesQh andUh such that they satisfy

Qh⊂L2(Ω) and Uh⊂H2(Ω)∩H01(Ω).

Applying Galerkins principle to (2.1) leads to the discrete variational problem for the functions (qh, wh, uh) ∈ Qh×Qh×Uh, which we immediately write in matrix- vector notation. We denote the vector representation of functions in these spaces by underlined versions of the corresponding symbols, i.e., forqh∈Qh the corresponding coefficient vector isq

h ∈RdimQh. With a slight abuse of notation, we use the same notation also for the right-hand-side vectors fh and ud,h, which are obtained by testing the corresponding linear functionals with the basis functions in Qh and Uh, respectively, see, e.g., [18, Section 6] for further details. Furthermore, operators with subscripthdenote matrix representations of the operators.

Using this notation, the discrete problem reads as follows.

Problem 3.1. Find (qh, wh, uh)∈RdimQh×RdimQh×RdimUh such that

(3.1) Ah

 qh

wh uh


 0 fh ud,h

 with Ah:=

αMh Mh 0

Mh 0 Kh

0 KhT MO,h


The exact Schur complement preconditioner (2.2) of the discretized system is

(3.2) S(Ah) =

αMh 0 0

0 α1Mh 0

0 0 MO,h+αKhTMh−1Kh


Under the mild conditionUh⊆Qh this preconditioner is symmetric positive definite.

This is a straight forward extension of [25, Lemma 4.4] by using the fact that (β·

∇uh, uh)L2(Ω)= 0. In this case, Theorem 2.2yields the following condition number bound:

κ (S(Ah))−1Ah

≤ cos(π/7)

sin(π/14) ≈4.05.


This preconditioner cannot be efficiently realized since the matrixMh−1 is dense. So, we use the following preconditioner motivated by the norms (2.4) on the discretization spaces instead:

(3.3) Sh:=

αMh 0 0

0 1αMh 0

0 0 MO,h+αBh


where Bh is the matrix representation of the linear operator B : H2(Ω)∩H01(Ω) → (H2(Ω)∩H01(Ω))0,

(3.4) hBu,ui˜ = (−ε∆u+β· ∇u+σu,−ε∆˜u+β· ∇˜u+σ˜u)L2(Ω)

onUh. On the continuous level, the operatorsB andK0M−1K coincide. In general, this does not carry over to the discrete case. The following lemma gives sufficient conditions that guarantee thatBh andKhTMh−1Kh coincide.

Lemma 3.2. If

(3.5) (−ε∆ +β· ∇+σ)Uh⊂Qh,

thenKhTMh−1Kh=Bh and thusS(Ah) =Sh.

Proof. Letuh∈Uhbe arbitrary but fixed with coefficient vector uh. The defini- tions yield

KhTMh−1Khuh, uh

= sup


hKhuh, whi2 hMhwh, whi

= sup


(−ε∆uh+β· ∇uh+σuh, wh)2L2(Ω)



Since (−ε∆ +β· ∇+σ)Uh ⊂Qh, the supremum is attained for wh=−ε∆uh+β·

∇uh+σuh, and we have KhTMh−1Khuh, uh

= sup


(−ε∆uh+β· ∇uh+σuh, wh)2L2(Ω)


=k −ε∆uh+β· ∇uh+σuhk2L2(Ω)=hBhuh, uhi. Therefore,KhTMh−1Kh=Bh and thusS(Ah) =Sh.

4. Isogeometric analysis. Due to requirementUh⊂H2(Ω)∩H01(Ω), we need a smooth discretization space. We achieve this by using IgA. We give a brief intro- duction to the approximation spaces in use. Let Sp,k,`(0,1) be the space of B-spline functions on the unit interval (0,1) which arek-times continuously differentiable and piecewise polynomials of degreepon a uniform grid with grid size 2−`. For the space of B-splines with maximum continuity, that is, withk=p−1, we only writeSp,`.

On the parameter domain Ω := (0,b 1)d, we use a tensor-product B-spline space, denoted by

Sp,k,`d :=






For ease of notation, we assume to have the same spline degreep, the same continuity kand the same number of uniform refinement steps`, for each spatial dimension. We assume that the domain Ω can be parametrized by a geometry mappingG:Ωb →Ω = G(bΩ) with the property

(4.1) k∇rGkL(bΩ)≤c1 and k(∇rG)−1kL(bΩ)≤c2, for r= 1,2,3, for some constants c1 and c2. The discretization space Sp,k,` on the domain Ω is defined using the pull-back principle as

Sp,k,`(Ω) :=

f◦G−1:f ∈Sp,k,`d .

For more information on IgA, see the survey article [5] and the references therein. We use spline spaces with maximum smoothness as the discrete state space and reduce the smoothness forQh accordingly. More precisely, we use

(4.2) Qh:=Sp,p−3,`(Ω) and Uh:=Sp,`(Ω)∩H01(Ω) with p≥2.

The following lemma shows that, if we consider the special case of box domains (if trivially parametrized) and constant convection, the condition of Lemma3.2holds.

Theorem 4.1. If Ω = (0,1)d, Qh =Sp,p−3,`d and Uh =Sp,`d ∩H01(Ω) and if the convectionβ is constant, then

S(Ah) =Sh and κ Sh−1Ah

≤ cos(π/7)

sin(π/14) ≈4.05.

Proof. For sake of simplicity, we restrict the proof to the two-dimensional case.


f0 ∈Sp−1,k−1,`(0,1) ∀f ∈Sp,k,`(0,1) together with

Sp−1,k−1,`(0,1)⊂Sp,k−1,`(0,1) and Sp,k,`(0,1)⊂Sp,k−2,`(0,1), see [5]. Letu∈Sp,k,`2 =Sp,k,`(0,1)⊗Sp,k,`(0,1). Then,


∂x21 ∈Sp−2,k−2,`(0,1)⊗Sp,k,`(0,1)⊂Sp,k−2,`2 ,


∂x22 ∈Sp,k,`(0,1)⊗Sp−2,k−2,`(0,1)⊂Sp,k−2,`2 , and, by combining these results, also

∆u= ∂2u

∂x21 +∂2u

∂x22 ∈Sp,k−2,`2 . Sinceβ= (β1, β2) is constant we also have

β· ∇u=β1∂u




∈Sp,k−2,`2 . To summarize, for everyu∈Sp,k,`2 , we have

−ε∆u+β· ∇u+σu∈Sp,k−2,`2 ,

hence Condition (3.5) in Lemma 3.2holds and we haveS(Ah) =Sh. The condition number bound follows from Theorem2.2.


If the domain is not a box domain or if the convection is not constant, then Theorem 4.1 cannot be applied. Nevertheless, the numerical results presented in Section7indicate that the spaces proposed in (4.2) work well even for more complex domains and variable convection.

5. Error estimates. In this section, we derive discretization error estimates for Problem 3.1. Let A(x,x) denote the bilinear form in Problem˜ 2.1, wherex,x˜ ∈ X denotes the triplet (q, w, u)∈L2(Ω)×L2(Ω)×H2(Ω)∩H01(Ω), with the norm

kxk2S =kqk2Sq+kwk2Sw+kuk2Su.

The discrete triplet (qh, wh, uh)∈ Qh×Qh×Uh is denoted byxh ∈Xh. Galerkin orthogonality reads as follows:

(5.1) A(x−xh,yh) = 0 ∀yh∈Xh.

Since Problem2.1is well-posed (Theorem2.2), we have boundedness (5.2) A(x,y)≤ckxkSkykS ∀x,y∈X

and inf-sup stability

(5.3) sup



kykS ≥ckxkS ∀x∈X,

where c/c =κ(S−1A). The boundedness also holds for a conforming discretization spaceXh ⊂Xwith the same constant c. If the condition (3.5) in Lemma3.2holds, then also the inf-sup holds with the same constantc. Using this and the ideas of [1], we derive the following discretization error estimate.

Lemma 5.1. If xh∈Xh is the solution to (3.1) for a discretization space satisfy- ing condition (3.5) and ifx∈Xis the solution of (2.1), then we have the estimate (5.4) kx−xhkS ≤(1 +κ(S−1A)) inf



Proof. Letyhandzhbe arbitrary functions inXh. Due to Galerkin orthogonality (5.1), we have

A(xh−yh,zh) =A(xh−x+x−yh,zh) =A(x−yh,zh).

Combining this with the boundedness condition (5.2) gives

A(xh−yh,zh)≤ckx−yhkSkzhkS ∀yh,zh∈Xh. Using the discrete inf-sup gives

ckxh−yhkS ≤ sup


A(xh−yh,zh) kzhkS

≤ckx−yhkS ∀yh∈Xh.

Finally, we use the triangle inequality andc/c=κ(S−1A) to obtain the desired result kx−xhkS ≤ kx−yhkS+kxh−yhkS ≤(1 +κ(S−1A))kx−yhkS ∀yh∈Xh.


Let Ω = Ω := (0,b 1)d and let the convection β be constant. Since we assume a trivial parametrization, we have the discretization spaces Qh = Sp,pư3,`d and Uh = Sp,`d ∩H01(bΩ).

To derive the error estimates, we assume that the solution of (2.1) satisfies the regularity assumption (q, w, u)∈H1(bΩ)×H1(bΩ)×H3(bΩ)∩H01(bΩ).

Now, we estimate the approximation error term in (5.4) from above. First, we observe that



kxưyhkS ≤ inf

qh∈Sp,pư3,`d αkqưqhk2L2(bΩ)+ inf



αkwưwhk2L2(bΩ)+ inf



The two first terms can be bounded by using the following approximation error esti- mate

(5.5) inf

qh∈Sdp,pư3,`kqưqhkL2(bΩ)≤ h 4√

3k∇qkL2(bΩ), see [22, Corollary 1]. The estimate for the last term


uh∈Sp,`d ∩H01(bΩ)


is more involved. Before handling this term, we need a convenient notation and some auxiliary approximation error estimates.

Notation 5.2. In what follows, c is a generic positive constant independent ofα, σ, β, ε, h and p, but may depend on the spatial dimension d and the observation domainO.

In Appendix B, we extend some of the results of [22, 23, 24]. The main result is summarized in the following theorem.

Theorem 5.3. Let Πp : H3(bΩ)∩H01(Ω)b →Sdp,`∩H01(bΩ) be the H2-orthogonal projector, wherep≥3 and`≥1. Then,

k∇2(IưΠp)ukL2(bΩ)≤ch k∇3ukL2(bΩ), (5.6)

k∇(IưΠp)ukL2(bΩ)≤ch2k∇3ukL2(bΩ), (5.7)

k(IưΠp)ukL2(bΩ)≤ch3k∇3ukL2(bΩ) ∀u∈H3(bΩ)∩H01(bΩ).


Remark 5.4. In the statement of Theorem5.3, the parameter domainΩ is consid-b ered. This result can be extended to physical domains if the corresponding geometry functionGis sufficiently smooth, cf. [24].

With Theorem5.3, we can derive an error estimate for our problem.

Theorem 5.5. Let Ω :=Ω := (0,b 1)d with d∈N, let β ∈Rd be constant and let (qh, wh, uh)∈Sp,pư3,`d ×Sdp,pư3,`×Sp,`d ∩H01(bΩ)withp≥3and`≥1, be the solution to (3.1). If(q, w, u)∈H1(bΩ)×H1(bΩ)×H3(bΩ)∩H01(bΩ)is the solution of (2.1), then we have the following estimates:

kqưqhkSq+kwưwhkSw+kuưuhkSu≤ ch

αk∇qkL2(bΩ)+ 1

√αk∇wkL2(bΩ)+√ αmax


σ+ 1





Proof. Let ˜uh:=Πpu. Then, we have for allu∈H3(bΩ)∩H01(Ω)b kuưu˜hk2Su =kuưu˜hk2L2(O)+αk(ưε∆ +β· ∇+σ)(uưu˜h)k2L2(bΩ). For the first term, we simply extend the norm fromOtoΩ and obtainb

kuưu˜hk2L2(O)≤ kuưu˜hk2

L2(bΩ)≤c h6k∇3uk2


For the second term, we use the triangle inequality and Theorem5.3, and we get αk(ưε∆ +β· ∇+σ)(uư˜uh)k2L2(bΩ)



≤c α ε2h2+kβk2h42h6

k∇3uk2L2(bΩ). Combining these results gives

(5.9) inf


uh∈Sp,`d ∩H01(bΩ)

kuưu˜hkSu ≤c√

α hmax

ε,kβkh,(σ+ 1/√

α)h2 k∇3ukL2(bΩ).

The theorem follows from combining Lemma 5.1 with Equation (5.5) and Equa- tion (5.9).

Remark 5.6. Theorem5.5is somewhat restrictive since it requires the domain to be a box domain and the convection to be constant. These requirements are needed for the proof of the discrete inf-sup stability. The numerical results presented in Section7 indicate that the restrictions are not needed.

6. Numerical experiments: Accuracy of the solution. In this section we compare the solution of the forward problem to the (state) solution of the optimal control problem to investigate the fact that the optimal control problem naturally in- troduces a non-standard Petrov-Galerkin method for the state equation. We consider a well-known one-dimensional problem [7]:

ư∂xu(x)ưε ∂xxu(x) = 0 in (0,1), u(0) = 0, u(1) = 1, whose exact solution is

u(x) =eưx/εư1 eư1/εư1.

This problem is used as state equation in our optimal control problem (Problem1.1), where we choose β =ư1, σ= 0, f = 0, ud to be the exact solution of the forward problem, and the boundary conditions on uto be as for the forward problem. The analytical solution of the optimal control problem is

q(x) = 0, w(x) = 0, u(x) = eưx/εư1 eư1/εư1. We use the discretization spaces

Qh=Sp,pư3,`(0,1) and Uh={uh∈Sp,`(0,1) : u(0) = 0, u(1) = 1}

for the optimal control problem, which satisfy the condition (3.5). We compare the numerical solution for the state with the numerical solution of the forward problem,


where we use Uh as trial and test space. No stabilization techniques are used. The diffusion is set to ε = 0.01 and α= 0.001 for the optimal control problem. Three observation domains are considered: Full observation, that is, O = Ω = (0,1), and partial observation on O = (0,14) and on O = (34,1). The numerical solutions are displayed in Figures 1 to4. The plots indicate that the forward solution is unstable for coarse discretizations. The non-physical oscillations start in the boundary layer and propagate into the remainder of the computational domain. These kinds of insta- bilities are often remedied by upwind and/or Petrov-Galerkin schemes [7]. We remark though that our Petrov-Galerkin like approach for the state equation does not resem- ble any of the common Petrov-Galerkin schemes for this equation, as far as we know.

The state solution (of the optimal control problem) does not have these instabilities.

In fact, the state operatorKhis discretized with a PetrovGalerkin method as the trial space isUh and the test space isQh.

0.0 0.2 0.4 0.6 0.8 1.0

x 0.0

0.2 0.4 0.6 0.8 1.0 1.2


Optimal control solution Forward problem solution Exact solution

0.0 0.2 0.4 0.6 0.8 1.0

x 0.0

0.2 0.4 0.6 0.8 1.0


Optimal control solution Forward problem solution Exact solution

Fig. 1.Full observation onO= (0,1)withp= 2(both),`= 4(left),`= 6(right).

0.0 0.2 0.4 0.6 0.8 1.0

x 0.0

0.2 0.4 0.6 0.8 1.0 1.2


Optimal control solution Forward problem solution Exact solution

0.0 0.2 0.4 0.6 0.8 1.0

x 0.0

0.2 0.4 0.6 0.8 1.0


Optimal control solution Forward problem solution Exact solution

Fig. 2.Partial observation onO= (0,14)withp= 2(both),`= 4(left),`= 6(right).

In Figure2, we consider the optimal control problem with observation on (0,14).

This is only a quarter of the whole domain, but it is located at the boundary layer.

The solutions for the state variable are almost identical to those obtained for the case of full observation.

Next we look at the solution where the observation domain is (34,1). Here the solution is almost constant (u(x)≈1). From Figures3and4, we see that the approx- imation is not good. In the left plot of Figure3, we see that the boundary layer is not captured. However, the error does not propagate into the observation domain. Forh- refinement, see Figure3(right) and Figure4(left), we observe that the approximation


0.0 0.2 0.4 0.6 0.8 1.0 x

0.0 0.2 0.4 0.6 0.8 1.0 1.2


Optimal control solution Forward problem solution Exact solution

0.0 0.2 0.4 0.6 0.8 1.0

x 0.0

0.2 0.4 0.6 0.8 1.0


Optimal control solution Forward problem solution Exact solution

Fig. 3.Partial observation onO= (34,1)withp= 2(both),`= 4(left),`= 6(right).

0.0 0.2 0.4 0.6 0.8 1.0

x 0.0

0.2 0.4 0.6 0.8 1.0


Optimal control solution Forward problem solution Exact solution

0.0 0.2 0.4 0.6 0.8 1.0

x 0.0

0.2 0.4 0.6 0.8 1.0


Optimal control solution Forward problem solution Exact solution

Fig. 4.Partial observation onO= (34,1)withp= 2,`= 8(left) andp= 4,`= 6(right).

improves slowly. Forp-refinement, see Figure4 (right), the approximation improves significantly. Since we use splines, increasing the spline degree by one means that the number of degrees of freedom is only increased by one, while eachh-refinement step doubles the number of degrees of freedom.

Remark 6.1. The effect of increasing the spline degree compared toh-refinement as shown in Figure4 is somewhat surprising. We are not completely sure why larger spline degrees are so effective. Unfortunately, the error estimate in Theorem5.5does not provide any explanation for this behavior. Further analysis is needed to explain this properly.

7. Numerical experiments for exact and inexact preconditioners. In this section, we analyze the convergence of Krylov space solvers when the proposed preconditioner is used. In the first subsection, we consider an exact realization of the preconditioner. A multigrid approximation is then considered in the second subsec- tion.

We have done the numerical experiments for two model domains, in both cases ford= 2. The first domain is a box-domain, more precisely, Ω is the unit square, see Figure 5 (left), which is parameterized with the identity function. For this domain, the conditions of Theorem 4.1 are satisfied. In Model problem 7.1, we have full observation and in Model problem7.2, the observation is restricted to the subdomain represented by the smaller area in Figure5 (left). In all model problems, the desired state ud is a step function with value ud = 1 inside a circle and with value ud = 0 outside the circle. The support ofud is shown as the dashed lines in Figure5.


Fig. 5.Computational domainsΩ, partial observation domainsO(dark blue), and support of desired state (inside of dashed circles).

Model problem7.1 (Unit square and constant convection with full observa- tion). LetΩ :=O:= (0,1)2 be the computational domain, which is also the observa- tion domain. The convection is β= (−2,1)and there is no reaction σ= 0 or source term f = 0. The desired state is

ud(x, y) =

(1 if (x−38)2+ (y−58)2161 0 otherwise.

The diffusionεand regularization parameterαwill vary.

Model problem7.2 (Unit square and constant convection with limited obser- vation). Let Ω = (0,1)2 be the computational domain andO= (14,34)2 be the obser- vation domain. The remainder of this problem is the same as for Model problem7.1.

Furthermore, we consider a non-trivial geometry Ω, which is a approximation of a quarter annulus by means of a B-spline parameterization, see Figure5(right). Again, Model problem7.3is a problem with full observation and the observation domain in Model problem7.4 is the smaller area in Figure 5 (right). We observe that for this domain, the conditions of Theorem4.1are not satisfied.

Model problem7.3 (Quarter annulus and varying convection with full obser- vation). Let Ω =O=G((0,1)2) withG: (0,1)2→R2 and

(7.1) G(bx) =

(1 +bx1)(1−xb2)(1 + 2(√

2−1)bx2) (1 +bx1)bx2(2√



be the computational domain, which is also the observation domain. The convection isβ = (y,1 +x2)and there is no reaction σ= 0or source term f = 0. The desired state is

ud(x, y) =

(1, if (x−x0)2+ (y−y0)2161 0, otherwise,

where(x0, y0) =G(38,58). The diffusionεand regularization parameterαwill vary.

Model problem7.4 (Quarter annulus and varying convection with limited ob- servation). Let Ω =G((0,1)2) be the computational domain andO=G((14,34)2) be the observation domain, whereGis as in (7.1). The remainder of this problem is the same as for Model problem7.3.


For all model problems, we consider a discretization of the optimality system using the spaces given in (4.2) as outlined in Section 4. The resulting linear system of equations


is solved using the MINRES method, preconditioned with the proposed Schur comple- ment preconditioner (3.3). We use a random initial guessxh,0. The stopping criterion is

krkk ≤10−8kr0k,

whererk :=bh− Ahxh,k denotes the residual andk · kis the Euclidean norm.

7.1. Results with exact preconditioner. In this section, we present the re- sults for the Schur complement preconditioner (3.3) when realized using a sparse Cholesky decomposition.

Table1shows the iteration numbers needed to reach the stopping criteria for full and partial observation (Model problems 7.1 and 7.2). In these tables, αand ε are varied, while p= 2 and` = 6 are fixed. In Table2, we set ε = 10−3 and vary the refinement level`andα. From the tables, we observe that for the partial observation problem, we need a few more iterations for small values ofα. This is probably because MO,h is singular in case of partial observation. The iteration numbers are relatively small for all considered values ofα,εand`. This is predicted by the theory as Model problems7.1and7.2satisfy the conditions of Theorem4.1.

Table 1

Iteration numbers: Model problem7.1(left) and7.2(right),p= 2,`= 6.

ε\α 100 10−3 10−6 10−9

100 12 26 60 72

10−3 15 47 26 11

10−6 15 46 26 11

10−9 15 46 26 11

ε\α 100 10−3 10−6 10−9

100 12 20 57 78

10−3 15 41 54 19

10−6 14 41 53 19

10−9 14 41 53 19

Table 2

Iteration numbers: Model problem7.1(left) and7.2(right),p= 2,ε= 10−3.

`\α 100 10−3 10−6 10−9

4 15 46 14 7

5 15 47 19 8

6 15 47 26 11

7 15 46 38 11

`\α 100 10−3 10−6 10−9

4 15 46 35 15

5 15 43 44 17

6 15 41 54 19

7 15 39 55 22

Next, we consider Model problems7.3and7.4, which are the problems where the computational domain is a quarter annulus. The iteration numbers are shown in Ta- bles3and4. Note that the conditions of Theorem4.1are not satisfied. Nevertheless, the iteration numbers are comparable with those of Tables1and2.

Remark 7.5. Forε= 1, the iteration numbers in Table1(and Table3) are grow- ing as α becomes smaller. This may appear strange since we have proven that the condition number (for Table 1) is less then 4.05. Describing convergence estimates for Krylov subspace methods in term of only the condition number can be misleading and/or insufficient, cf. [16]. The different iteration numbers for various values of ε and α can be explained by the distribution of the eigenvalues. For small iteration numbers the eigenvalues are more clustered. Forε= 1, the iteration numbers starts decreasing whenα <10−9.



In terms of whether the courts of the home nations of the UK would adopt a similar approach, it is, perhaps, important to note that that we should not become overly fixated on the

„ The role of small airway dysfunction in asthma control and exacerbations: a longitudinal, observational analysis using data from the ATLANTIS study.. Kraft M

A scale analysis indicated that indices that represent terrain convergence (such as ln(a) or ln(a/tan(b))) explain variability at all scales from 10 m up to the catchment scale

Our convergence rate analysis of this class of accelerated Newton-Landweber methods contributes to the analysis of Newton-type regularization methods in two ways: first, we show

In this paper, we consider a Bayesian approach to parameter estimation, using Markov chain Monte Carlo (MCMC) methods, which is capable of dealing both with continuous time HMMs as

The analysis is presented within the proof of forthcoming Proposition 1, whose key conclusion is that under the scenario in consideration, sets A λ must have small

A first key observation in our analysis, which happens to be useful for the imaging later on, is that the Minnaert resonance and the sequence of body wave resonances are

A key point of the analysis is to show that the Foldy-Lax field appearing in the meso-scale approximation, derived in [14], is a discrete form of a (continuous) system of

Time for assembling and solving the systems that generate u h and y h as well as the time spent on e/w evaluation of error, majorant, and residual error estimator w.r.t..

A priori error analysis T-matrix based multi-particle scattering (adaptive) control of (accuracy of) forward solvers in offline stage Compression of T-matrix: optimal

In Section 4, we introduce the algebraic analysis ap- proach to linear systems theory and, using homological algebra techniques, we ex- plain that this approach is an

In our approach, we construct the preconditioner based the complete LDU factorization for the the linearized 3 × 3 coupled system and the ap- proximated Schur complement itself is

1) A priori error estimates for space-time finite element discretization of parabolic optimal control problems (in cooperation with D. We developed a priori error analysis for

An a posteriori error analysis of adaptive finite element methods for distributed elliptic control problems with control constraints.. ESAIM

The paper presents a new approach to EU governance by stressing the interdependence of governance and integration. It suggests that EU gov- ernance is not just shaped by the

An analysis of the financial account by sectors shows once more that, with a volume of EUR 6,690 million (mainly from the sale of Austrian securities abroad), the government has

Visualization and Visual Analysis of Multi-faceted Scientific Data:..

Model 2 includes the same dummy variables for secondary formal debt instruments but replaces the bank loan dummy with a dummy variable for broad bank debt (bank loan, overdraft,

Analysis of the Constitutional Treaty of the European Union shows that there is a serious discrepancy between the power of Member States as computed by the two main power

Travel distances and times decrease with the size of the municipalities, but even for smaller municipalities with less than 2,000 inhabitants the mean distance seems to be

An analysis of the situation in Austria has to take account of the fact that the different tax treatment of pub- lic- and private-sector wages leads to a distortion of

By contrast, the article at hand explores goods and service exports, starting with an analysis of changes in market shares that is meant to highlight how the patterns of

&#34;Numerical control&#34; (2) means the automatic control of a process performed by a device that makes use of numeric data usually introduced as the operation is in progress