Reduced basis methods for a simultaneous space-time variational formulation of parametric parabolic PDEs

(1)

MSc Mathematics

Master Thesis

Reduced basis methods for a

simultaneous space-time variational

formulation of parametric parabolic

PDEs

Author: Supervisor:

Femke Madsen

prof. dr. R. P. Stevenson

Examination date:

(2)

Abstract

In this thesis, reduced basis methods for solving parametric PDEs will be discussed. First, the theory for problems that involve a coercive bilinear form, which arise from variational formulations of elliptic PDEs, is given. Next, we compare two ways of solv-ing parabolic PDEs with reduced basis methods: ussolv-ing a time-marchsolv-ing scheme and a space-time variational formulation. Numerical experiments will be presented and we will compare the outcomes. Finally, conclusions about the comparison will be discussed.

Title: Reduced basis methods for a simultaneous space-time variational formulation of parametric parabolic PDEs

Author: Femke Madsen, femke.madsen@student.uva.nl, 10754164 Supervisor: prof. dr. R. P. Stevenson

Second Examiner: dr. C. C. Stolk Examination date: May 13, 2020

Korteweg-de Vries Institute for Mathematics University of Amsterdam

Science Park 105-107, 1098 XG Amsterdam http://kdvi.uva.nl

(3)

Introduction

Numerous physical phenomena can be described by partial differential equations (PDEs), for example heat-transport, quantum mechanics and electrodynamics. Since they play such a central role in these physical models, and thus also in every day life, already a lot is known about solving PDEs. Sometimes analytical solutions can be obtained, but in most cases this is not possible. We then turn to numerical solutions, which can be obtained by a discretization of the PDE. Some examples of such discretization techniques for solving PDEs are Finite Element methods, Finite Volume methods or Finite Difference methods. All these techniques result into solving a very high dimensional sparse discrete system. For solving the high dimensional discrete system a lot of storage and computation power is needed. Also, the computation time for obtaining one solution can be very high.

One can imagine that if solutions for many instances of a PDE are required, the use of the above techniques give rise to problems. This can for example be the case in multi-query problems, where more simulations of a problem with a varying quantity are required. Think for example about parameter studies, some optimization problems or uncertainty quantification. Also, the above techniques are not suitable if computation of a solution is required really fast. This can be the case in real-time problems, where we for example have a problem with interaction with humans or when real-time processes are involved. For these classes of problems another discretization technique is developed: the Reduced Basis method.

The general setting for Reduced Basis methods is the following: we assume to have a problem that is dependent on a parameter µ, which lies in some parameter space P, and we are looking for a parameter-dependent solution u(µ) in some space V . The manifold of all solution u(µ) will be denoted by M. The crucial insight of the Reduced Basis method is that in many cases the solution u(µ) depends smoothly on µ, so that this solution manifold M can be approximated by some low dimensional space, which we will denote by VN. Thus, if we could find a space VN that approximates the solution

manifold M well, we could instead look for an approximate solution uN(µ) ∈ VN, which

now lies close to the real solution u(µ). Since the approximation space is supposed to be small, the solutions uN(µ) can be obtained very quickly for a lot of µ’s and thus this

method is suitable for multi-query and real-time problems.

The construction of such an approximation space, also called the Reduced Basis space or RB-space, is generally done by computing some so called snapshots, which are solu-tions u(µ(i)) for some carefully chosen µ(i) ∈ P, i = 1, . . . , N , after which VN is defined

as the span of those snapshots. Of course, the real solutions u(µ(i)) cannot be obtained analytically. Instead we take some approximate solution uH(µ(i)), which is obtained

(6)

close to the real solution u(µ(i)).

This construction of the approximation space VN, spanned by some very good

approx-imations {uH(µ(i))}, may take a very long time: we need to obtain N solutions using for

example the Finite Element method. The crucial insight here is that this construction of this approximation space needs to be done only once. Therefore, we can divide the computational procedure of the Reduced Basis method into a so called offline and online phase.

During the offline phase, the RB-space VN is constructed and some high dimensional

quantities are precomputed. The operation count here depends on H, defined as being the dimension of the finite element space, but only µ-independent quantity’s are precom-puted. This means that the offline phase should only be performed once. Then, during the online phase, µ-dependent approximate solutions uN(µ) can be obtained. Here the

operation count should not depend on H but on N : the online phase can be performed very fast. The Reduced Basis method thus only makes sense if we need the solution for a lot of µ’s, otherwise the high cost of the offline phase causes that it is more efficient simply to solve the problem for the individual µ’s of interest.

The Reduced Basis method described above is very well developed for problems which involves a continuous, coercive bilinear form. Such problems arise for example from a variational formulation of an elliptic PDE. In this case, rigorous a-posteriori error bounds can be given which are crucial for selecting the snapshots, as well as to certify the accuracy of the reduced basis solution. Furthermore, upper bounds for the effectivity factor of the a-posteriori error bounds can be given, see [13], and exponential decay for the error between the real and RB-solution can be shown, as in [2]. Numerical experiments as in [13] show that RB-methods also work very well from a computational point of view.

If we now move to the case where we do not have a coercive bilinear form, which for example is the case for parabolic PDEs, the theory is less well developed. In most references, see e.g. [13] and [9], a time-marching scheme is used to solve parabolic PDEs with Reduced Basis methods. This means that the time-variable is discretized and then for each time step tk a remaining elliptic problem has to be solved with say the Finite Element method. In this case the reduced basis method is applied to replace this finite element space by a reduced basis space of much smaller dimension. On the other hand, no complexity reduction takes place in the time direction, since for any time step still a, now low dimensional, system has to be solved.

Another way of solving the parabolic PDEs using RB-methods, is via a so called space-time variational formulation. In this formulation, the time-variable is simply seen as an extra space-dimension. Application of the RB-method on such a formulation was already given in [9].

In this thesis, we are going to look at another space-time formulation, recently intro-duced in [8], and apply RB-methods on this. This space-time formulation reformulates the problem of a parabolic PDE to a problem which involves a continuous and coercive bilinear form. On this reformulation, we can thus apply the whole well developed theory for problems with continuous and coercive bilinear forms. We will give numerical results

(7)

of the application of RB-methods on both the time-marching scheme and the space-time variational formulation. We will see that the RB-space in the space-time variational formulation is of higher dimension than the RB-space in the time-marching case. But because of the fact that in the online phase of the time-marching scheme we still need to find an RB-solution for every time-step, the online-phase of the space-time variational problem will be faster.

This thesis is divided into five chapters. In the first chapter, some basic notions on functional analysis and numerical analysis are given that are needed in this thesis. In the second chapter, the theory on Reduced Basis methods for problems that involve continuous and coercive bilinear forms is presented. Here we will follow mainly [13], from which we also reproduce some numerical experiments. In the third chapter, we will discuss both the time-marching scheme, as in [13], [9], and the space-time variational formulation, as in [8], for parabolic PDEs. We will also give the theory for the application of the Reduced Basis methods in both cases. In the fourth chapter, we will discuss our numerical experiments and compare the two methods. Finally, in the last chapter, we will draw some conclusions and give some ideas for further research.

(8)

1 Basic notions in Functional and

Numerical Analysis

This first chapter is meant to give a very short introduction in some Functional Analysis and Numerical Analysis topics that are used in this thesis. For most topics a reference for further reading is given. We advice readers with no background in these topics to read the references as well.

Readers with strong background in these topics can skip the first chapter or use the first chapter as a reference while reading the thesis. Each time a topic is used in the thesis, we will mention the corresponding sections from this chapter.

1.1 Bilinear Forms

In the theory of solving partial differential equations, bilinear forms play a central role as they arise in weak formulations of e.g. parabolic and elliptic PDEs. In this subsection we will give the definition and some properties of bilinear forms. We conclude with the Lax-Milgram theorem.

Definition 1.1. Let H be a Hilbert-space. A map a(·, ·) : H ×H → R is called a bilinear form if it is linear in both components, i.e., for v1, v2, v3∈ H and α ∈ R it holds that

a(v1+ v2, v3) = a(v1, v3) + a(v2, v3), a(v1, v2+ v3) = a(v1, v2) + a(v1, v3),

a(αv1, v2) = αa(v1, v2), a(v1, αv2) = αa(v1, v2).

If a bilinear form a is also dependent on a parameter µ ∈ P, where P is some parameter space, we call a a parametric bilinear form and denote this with a(·, ·; µ).

Definition 1.2. A parametric bilinear form a(·, ·; µ) over a Hilbert space H is continuous if there exist γ(µ) ∈ R with

γ(µ) := sup

u,v∈H\{0}

a(u, v; µ) kukkvk < ∞.

It is uniformly continuous with respect to µ if sup_µ∈Pγ(µ) ≤ ∞, i.e. there exist a ¯γ ≤ ∞ such that γ(µ) ≤ ¯γ for all µ ∈ P.

Definition 1.3. A parametric bilinear form a(·, ·; µ) over a Hilbert space H is called coercive if there exist α(µ) ∈ R with

α(µ) := inf

u∈H\{0}

a(u, u; µ) kuk2 > 0.

(9)

The coercivity is uniform with respect to µ if infµ∈Pα(µ) > 0, i.e. there exists an ¯α

such that α(µ) ≥ ¯α > 0 for all µ ∈ P.

The following theorem is used widely to show well-posedness of weak formulations of PDEs.

Theorem 1.4. (Lax-Milgram) Let H be a Hilbert space with inner product h·, ·i and assume that a : H × H → R is a continuous, coercive bilinear form on H with coercivity constant α. Then, given any functional φ ∈ H0, there exist a unique element u ∈ H such that

a(u, v) = φ(v), ∀v ∈ H. Moreover, the solution u can be bounded by

kuk ≤ kφkH0 α . Proof. For a proof, see for example [1], Theorem 2.3.

1.2 Sobolev spaces and Bochner spaces

Remember that for Ω ⊂ Rn, L2(Ω) is defined as the space of equivalence classes of the

a.e. square integrable functions over Ω, i.e.

L2(Ω) = f : Ω → R | Z Ω |f (x)|2_dx 1/2 < ∞ .

Similarly, Bochner spaces, which are Lp spaces over function spaces, can be introduced.

Definition 1.5. Let I ∈ R be an interval and H a Hilbert space. The Bochner space L2(I; H) is then defined as

L2(I; H) = f : I → H | Z I kf (t)k2dt 1/2 < ∞ . This is a Hilbert space with inner product

hf, gi_L₂_(I:H) = Z

I

hf, giHdt.

It can be shown that the dual space of the Bochner space L2(I; H) is equal to L2(I; H0).

In this thesis we will only need L2Bochner spaces, but of course Bochner spaces Lp(I; H)

can also be defined for arbitrary p ∈ [1, ∞].

Further, remember that the Solobev spaces Hk(Ω) ⊂ L2(Ω), k ∈ N are defined as all

the functions from Ω to R for which all the weak derivatives up to order k lie in L2(Ω),

(10)

Definition 1.6. Let Ω be a domain in Rn. The Solobev space Hk(Ω) is then defined as Hk(Ω) = {f : Ω → R | Dαf ∈ L2(Ω), |α| ≤ k}.

Here, the derivative Dα of f is taken in the weak sense. This means that g is the α-th weak derivative of f if (−1)|α| Z Ω g φ dx = Z Ω f (Dαφ) dx,

for all test functions φ ∈ C₀∞(Ω). We then write g = Dαf . The Sobolev space Hk(Ω) is a Hilbert space with inner product and norm

hf, gi_Hk_(Ω)= X |α|≤k hDαf, DαgiL2(Ω), kf kHk(Ω) = hf, f i 1/2 Hk_(Ω).

The space H₀k(Ω) is defined the closure of C₀∞(Ω) in Hk(Ω) and is thus given by all the functions in Hk_{(Ω) that are zero on the boundary ∂Ω. The norm on H}k

0(Ω) is equivalent

to the following semi-norm on Hk(Ω)

kf k_Hk 0 = X 1≤|α|≤k hDαf, Dαf iL2(Ω) 1/2 .

The dual space of Hk(Ω) is denoted as H−k(Ω).

Analogue to the definition of Bochner spaces, we can define the space Hk(I; H) as Hk(I; H) = {f : I → H | Dαf ∈ L2(I; H)}.

Here a function f : I → H has a weak-time derivative g = Dαf if (−1)|α| Z I g φ dt = Z I f (Dαφ) dt,

for all φ ∈ C₀∞(I). An integral over a function that takes values in H is defined as in [15], Section 1.2.

For more information on Sobolev spaces, see for example [10], Section 4.2. For more details on Bochner spaces, see [15].

1.3 Numerically solving elliptic PDEs

In this section we will very briefly discuss how to solve an elliptic PDE using Finite Elements approximations. In this thesis we make use of Finite Element approximations for elliptic PDEs a couple of times, and it is a very important ingredient for solving a parametrized PDE using the Reduced Basis method. We will here follow [21], but more details and precise formulations on Finite Elements can be found in [7].

(11)

In most cases, elliptic PDEs are solved by converting the problem to its weak formu-lation. For this, let Ω be a bounded open set and suppose we want to solve the second order elliptic PDE

− n X i,j=1 ∂ ∂xj (aij ∂u ∂xi ) + n X i=1 bi(x) ∂u ∂xi + c(x)u = f (x), x ∈ Ω, (1.1)

where the coefficients aij satisfy aij ∈ C1( ¯Ω) and satisfy the uniform ellipticity condition,

and bi, c, f ∈ C( ¯Ω). The uniform ellipticity condition is satisfied for aij if n X i,j=1 aijyiyj ≥ C n X i=1 y2_i, y = (y1, . . . , yn) ∈ Rn, x ∈ ¯Ω,

for some C > 0. Also assume that we have Dirichlet boundary value conditions u(x) = 0 x ∈ ∂Ω.

A solution u of the elliptic problem with u ∈ C2(Ω) ∩ C( ¯Ω) is called a classical solution. This classical solution exists provided that ai,j, bi, c, f and Ω are smooth enough.

In a lot of cases, for example when the function f on the right hand side is a dis-continuous function or when Ω has a re-entrant corner, it is not possible to find such a classical solution. It is then possible to look at the weak formulation of our elliptic PDE. For this, we first multiply both sides with a test function v ∈ C₀1(Ω) and then integrate over Ω. Then, by integration by parts and the fact that v = 0 on the boundary of Ω, it follows that the equation in (1.1) can be rewritten as

n X i,j=1 Z Ω aij ∂u(x) ∂xi ∂v(x) ∂xj dx + n X i=1 Z Ω bi(x) ∂u ∂xi v(x) dx+ Z Ω c(x)u(x)v(x) dx = Z Ω f (x)v(x) dx, v ∈ C₀1(Ω). (1.2)

Now, instead of the requirement that u ∈ C2(Ω) ∩ C( ¯Ω), we now require that u ∈ H₀1(Ω), which is a weaker requirement. Also, the requirements on ai,j, bi, c, f can be weakened.

Thus, this now gives us that the weak solution of the elliptic problem is given by the following:

Definition 1.7. Let aij, bi ∈ L∞(Ω), i, j = 1, . . . , n, c ∈ L∞(Ω) and f ∈ L2(Ω). A

function u ∈ H₀1(Ω) is called a weak solution of the elliptic problem if it satisfies

n X i,j=1 Z Ω aij ∂u ∂xi ∂v ∂xi dx + n X i=1 Z Ω bi(x) ∂u ∂xi v dx+ Z Ω c(x) u v dx = Z Ω f (x) v dx, (1.3)

(12)

If we define the bilinear form a : H₀1× H1 0 → R as a(u, v) = n X i,j=1 Z Ω aij ∂u ∂xi ∂v ∂xi dx + n X i=1 Z Ω bi(x) ∂u ∂xi vdx + Z Ω c(x)uvdx

and the linear form g : H₀1 → R as g(v) =

Z

Ω

f (x)v(x)dx,

then the problem in the definition of the weak solution can be rewritten as: find a solution u ∈ H₀1(Ω) such that

a(u, v) = g(v), ∀v ∈ H₀1(Ω). (1.4) If the bilinear form a is coercive, which is certainly the case if aij satisfies the uniform

ellipticity condition, bi≡ 0 and c ≥ 0 on Ω, this problem now has a unique solution due

to the Lax-Milgram Theorem 1.4. For the case when bi 6= 0, we refer to [21].

1.3.1 Galerkin approximation

Generally, it is not possible to find a solution of (1.4) analytically, since we are searching a solution in an infinite dimensional function space. We thus want to find an approximate solution to (1.4). To do this, we replace the function space in which we are looking for a solution, by a finite-dimensional subspace VH of the function space. This

finite-dimensional space VH can be constructed using the theory of Finite Elements, see [7]

Chapter 1.

In the case of the elliptic problem (1.3) (and in the numerical experiments in this thesis) this finite-dimensional space VH consists of piecewise polynomials of maximal

degree p with respect to a mesh on the domain Ω, that are zero on the boundary ∂Ω. This means that VH ⊂ H01(Ω). We are thus looking for a solution of the Galerkin system:

find uH ∈ VH such that

a(uH, vH) = g(vH), ∀vH ∈ VH. (1.5)

This approximation method is also called the finite element approximation method and results in solving a matrix-vector equation. The bilinear form in (1.5) is coercive and continuous if it was coercive and continuous in (1.3), since coercivity and continuity are inherited onto subspaces. The problem then satisfies the conditions of the Lax-Milgram Theorem 1.4. In this case (1.5) has a unique solution.

It is further immediate that for a solution u of (1.4) and a solution uH of (1.5) it holds

that

a(u − uH, vH) = 0 for all vH ∈ VH.

This property is also known as Galerkin Orthogonality. This can be used to prove Cea’s Lemma, which is a very usefull statement for error analysis:

(13)

Lemma 1.8. (Cea’s Lemma) Let u be a solution of (1.4) and uH be a solution of (1.5), then ku − u_Hk_H1 0(Ω) ≤ c1 c2 min vH∈VH ku − v_Hk_H1 0(Ω),

for some constants c1, c2 ∈ R.

Proof. By coervivity and continuity of the bilinear form a and the Galerkin Orthogo-nality, it follows that

ku − uHk2_H1 0(Ω) ≤ 1 c2 a(u − uH, u − uH) = 1 c2 (a(u − uH, u − vH) + a(u − uH, vH − uH)) ∀vH ∈ VH = 1 c2 a(u − uH, u − vH) ≤ c1 c2 ku − u_Hk_H1 0(Ω)ku − vHkH01(Ω).

Dividing both sides by ku − uHkH1

0(Ω) gives the result for all vH ∈ VH, so in particular

for the minimum.

1.4 Least Squares discretization

In this section we give a short introduction into Least Squares discretization. For further reading see for example [11].

Let H, V be two Hilbert spaces and let G ∈ L(H, V0). Further, assume that G is a homeomorphism from H to ran(G) ⊂ V0, i.e.

kGwkV0 _{h kwk}_H w ∈ H. (1.6)

Suppose now that, for given g ∈ V0 we want to find a u ∈ H such that Gu = g.

Since ran(G) ⊂ V0, the above problem does not always have a solution. We will therefore rewrite the problem to a least-squares formulation, which will then result in solving a system that involves a continuous and coercive bilinear form which has a unique solution. This means that, instead of solving Gu = g, we want to find a solution u ∈ H satisfying

u = argmin_w∈H1

2kGw − gk

2

V0. (1.7)

This is thus the minimizer of the least-squares functional. The minimizer u of (1.7) can be found by putting the derivative with respect to w equal to zero. This then results in solving the Euler-Lagrange equation for u, i.e. find u ∈ H such that

(14)

Note that, if we denote G∗ as the Hilbert-space adjoint of G, the above equation is exactly given by the normal equations in H, thus by

G∗Gu = G∗g, in H,

which is a very well known technique to solve least-squares problems.

Since assumption (1.6) holds, the above system (1.8) has a unique solution u by the Lax-Milgram Theorem 1.4. The solution u is therefore the unique minimizer of the least-squares functional (1.7). Obviously, if g ∈ ran G, then the solution u of (1.7) satisfies Gu = g.

Again, it is generally not possible to find a solution for (1.7) analytically. We are therefore going to find an approximate solution. For this, assume that k · kV0 is a

computable quantity (e.g. when V ' V0 is an L2-space). Given a closed linear subspace

Hδ of H, e.g. Hδ being a finite dimensional subspace, let uδ = argmin_w∈Hδ

1

2kGw − gk

2 V0

be the minimizer of the least-squares functional over Hδ. From the above, this is the unique solution in Hδ of the Galerkin system

hGuδ, GwiV0 = hg, Gwi_V0, w ∈ Hδ.

From Cea’s Lemma 1.8 it further follows that

ku − uδkH ≤ Cminw∈Hδku − wk_H,

for some contant C. This means that uδ is the best quasi-approximation with respect to the H-norm. If now g ∈ ran G, it follows from (1.6), that

ku − uδ_k

H ' kG(u − uδ)kV0 = kg − Guδk_V0.

The right-most expression is now computable since we assumed that k·kV0 is computable.

Hence it gives us a reliable and efficient a-posteriori estimator of the error u − uδ _{in the}

H-norm.

1.5 Crank-Nicolson time-stepping

In Section 1.3 we saw how to solve a elliptic PDE numerically. If we now move to time-dependent PDEs, in particular parabolic PDEs, we also have to deal with the time-parameter. We now have to discretize the time-derivative as well. In this section, we will show how to do this using a Crank-Nicolson time-stepping scheme. We follow here the lecture notes of [21], which can also be used as reference for further reading.

We will consider the following model problem in one space dimension, but all the theory can easily be generalised to different problems, also in more space dimensions.

(15)

Let us denote the space-time cylinder Q = Ω × (0, T ], where Ω = (0, 1), T > 0. We then want to find a solution u(x, t) of the following second order parabolic problem,

∂u ∂t − ∂2u ∂x2 + ∂u ∂x+ u = f (x, t), (x, t) ∈ Q u(0, t) = 0, u(1, t) = 0, t ∈ (0, T ) u(x, 0) = u0(x), x ∈ Ω. (1.9)

Again, we reformulate this problem to its weak form. In order to do so we multiply both sides with a test function v and integrate over space only. This results in the following problem: find u(t) ∈ H₀1(Ω) such that for any t ∈ (0, T ],

Z Ω ∂u ∂t v dx + Z Ω ∂u ∂x ∂v ∂x+ ∂u ∂x v + u v dx = Z Ω f (x, t) v dx. Z Ω u(x, 0) v dx = Z Ω u0 v dx, (1.10)

for all test functions v ∈ H₀1(Ω). The above formulation is also referred to as a semi-variational formulation.

To find an approximate solution of this semi-variational formulation, we discretize Q as follows. For K, H ∈ N, H ≥ 2, K ≥ 1, let h = 1/H be the mesh size in the space direction and ∆t = T /K be the mesh size in the time direction. Define now the uniform mesh ¯QH,K on ¯Q as

¯

QH,K = {(xj, tk) : xj = jh, 0 ≤ j ≤ H; tk= k∆t, 0 ≤ k ≤ K}.

Let VH ⊂ H01(Ω) denote the approximation space consisting of all continuous piecewise

linear functions defined on the mesh on Ω, that vanish at the boundary points.

The system (1.10) can now be approximated by the finite element method, using the approximation space VH for the space direction. The time derivative can be

approxi-mated in two ways: we can use the forward divided difference or the backward difference in t. These methods are respectively called the forward Euler scheme and backward Euler scheme.

Forward Euler Scheme: Find uk_H ∈ VH, 0 ≤ k ≤ K such that

uk+1_H − uk H ∆t , vH L2(Ω) + a(uk_H, vH) = hf (·, tk), vHiL2(Ω) ∀vH ∈ VH, hu0_H − u0, vHiL2(Ω) = 0 ∀vH ∈ VH. (1.11) Here uk

H ∈ VH is an approximation of u(·, tk) ∈ VH and a(·, ·) is defined as

a(u, v) = Z 1

0

(16)

Rewriting (1.11) gives us, for all 0 ≤ k ≤ K − 1, huk+1_H , vHiL2(Ω) = hu k H, vHiL2(Ω)− ∆ta(u k H, vH) + ∆thf (·, tk), vHiL2(Ω) ∀vH ∈ VH hu0 H, vHiL2(Ω) = (u0, vH)L2(Ω) ∀vH ∈ VH.

This system is not always stable, but it can be shown, see [21], that it is stable under some restrictions on the size of the time-steps.

Let now {ψi} be the one-dimensional piecewise linear element basis functions of the

approximation space VH. Given uk_H, to find uk+1_H , we need to solve a system with matrix

M , which has entries Mi,j = hψj, ψii. This matrix is symmetric positive definite and has

size (H − 1) × (H − 1). The same matrix M has to be inverted when determining u0_H. Backward Euler Scheme: Find uk_H ∈ V_H, 0 ≤ k ≤ K such that

uk+1_H − uk H ∆t , vH L2(Ω) + a(uk+1_H , vH) = hf (·, tk+1), vHi ∀vH ∈ VH, hu0_H − u₀, vHiL2(Ω) = 0 ∀vH ∈ VH. (1.12) Again, uk

H and a(·, ·) are defined as in the Forward Euler Scheme. Rewriting (1.12) gives

us, for all 0 ≤ k ≤ K − 1, huk+1_H , vHiL2(Ω)+ ∆ta(u k+1 H , vH) = hukH, vHiL2(Ω)+ ∆thf (·, t k+1_{), v} HiL2(Ω) ∀vH ∈ VH hu0_H, vHiL2(Ω)= hu0, vHiL2(Ω) ∀vH ∈ VH.

This backward Euler scheme is always stable (see again [21]), so no restrictions have to be made on the size of the time-steps.

Now, given uk_H, to find uk+1_H , we need to solve a system with matrix A, which has entries Ai,j = hψj, ψii + ∆t(hψ0j, ψi0i + hψj0, ψii + hψj, ψii). Here again {ψi} are the

one-dimensional piecewise linear element basis functions of the approximate space VH. Note

that when determining u0_H, still only the matrix M has to be inverted.

Generally, the two above schemes can be studied by looking at the following θ-scheme: find uk_H ∈ VH, 0 ≤ k ≤ K, such that

uk+1_H − uk H ∆t , vH L2(Ω) + a(θuk+1_H + (1 − θ)uk_H, vH) = hθf (·, tk+1) + (1 − θ)f (·, tk), vHiL2(Ω) ∀vH ∈ VH, hu0_H − u₀, vHiL2(Ω) = 0 ∀vH ∈ VH, (1.13)

where 0 ≤ θ ≤ 1. Note that for θ = 0, this gives the forward Euler scheme and for θ = 1 this gives the backward Euler scheme. For θ = 1₂, this scheme also called the Crank-Nicolson scheme.

It can be shown, see for example [21], that the scheme in (1.13) is unconditionally stable for θ ∈ [1/2, 1]. This means that it is stable without any restrictions on the size of the time-steps.

(17)

2 Reduced Basis methods

In this chapter the details of the Reduced Basis method in the case of problems with coercive and continuous bilinear forms, that for example arise from variational formula-tions of elliptic PDEs, will be given. We will discuss the error analysis, the Offline-Online decomposition and the construction of the RB-space. We further give some reproduc-tions of numerical experiments. In this chapter we will follow [13] closely and adapt most of the notation.

2.1 General problem

2.1.1 Setting

Let H be a real, separable Hilbert space with inner product h·, ·i, norm k·k and dual space H0. Let further V ,−→ H ' H0_,−_{→ V}0 _{be a so called Gelfand triple (see e.g. [3]). The inner}

product on V will be denoted as h·, ·iV and its norm by k · kV. We also assume to have

a parameter domain, which we will denote by P ⊂ Rp and assume to have a parametric bilinear form a(·, ·; µ) on V and parametric linear forms l(·; µ), f (·; µ), which lie in the dual space V0 for all µ ∈ P. We further assume that a, l, f are uniformly continuous with respect to µ, with continuity constants ¯γ, ¯γland ¯γf and that a is uniformly coercive

with respect to µ with coercivity constant ¯α (see also Section 1.1). The general problem we now want to solve can be formulated as follows.

Problem 2.1. For µ ∈ P, find a solution u(µ) ∈ V such that a(u(µ), v; µ) = f (v; µ), ∀v ∈ V.

Often we are interested in a functional evaluated at the solution rather then the solution itself. In that case we are additionally looking for an output s(µ) ∈ R of the form

s(µ) = l(u(µ)).

The following proposition states that Problem 2.1 is well-posed and is stable.

Proposition 2.2. The Problem 2.1 admits a unique solution u(µ) ∈ V . Further, the solution and output satisfy

ku(µ)k_V ≤ kf (·; µ)kV0 α(µ) ≤ ¯ γf ¯ α, |s(µ)| ≤ kl(·; µ)k_V0ku(µ)k_V ≤ ¯γlγ¯f ¯ α .

(18)

Proof. The uniqueness, existence and the first bound for u(µ) follow immediately from the Lax-Milgram Theorem 1.4. The second and parameter independent bound for u(µ) now follows from the assumptions of uniform continuity and coercivity. The bound for the output functional follows by definition of the output functional s and uniform continuity of l.

We will refer to a solution u(µ) of the general problem 2.1 as the real solution. Now that we have established solvability of the problem, we can introduce the solution man-ifold M as follows.

Definition 2.3. The solution manifold M of the Problem 2.1 is given by M := {u(µ) | u(µ) solves Problem 2.1 for µ ∈ P} ⊂ V.

Another assumption that needs to be made is so called parameter-separability of the (bi)linear forms a, l, f . This assumption is a very important property for an efficient implementation of the RB-method, which will be discussed in detail in Section 2.3. Definition 2.4. A parametric bilinear form a is called parameter-separable if there exist parameter-independent continuous bilinear forms aq(·, ·) : V × V → R and coefficient

functions θa

q : P → R for q = 1, . . . , Qa∈ N such that

a(u, v; µ) =

Qa

X

q=1

θ_qa(µ)aq(u, v), µ ∈ P, u, v, ∈ V.

A similar definition can be given for a linear form l with coefficient functions θl q and

parameter-independent continuous linear forms lq.

From now on a, l, f are assumed to be parameter-separable with (bi)linear forms aq, lq, fq, coefficient functions θqa, θql, θ

f

q and numbers of components Qa, Ql, Qf.

2.1.2 Properties of M

Reduced Basis methods rely on the fact that the solution manifold M often can be very well approximated by a low dimensional subspace. We are therefore interested in properties of M that may give some insights on its complexity or how well the solution manifold may be approximated by a low-dimensional subspace. Note first that Propo-sition 2.2 implies the boundedness of M, since the solution is bounded. Next, we will also see that Lipschitz-continuity of the solution holds.

Proposition 2.5. Let a, l, f be such that the coefficient functions θa q, θql, θ

f

q are Lipschitz

continuous with respect to µ. It then follows that the (bi)linear forms a, l, f , the solution u(µ) and the output s(µ) are also Lipschitz-continuous with respect to µ.

(19)

Proof. We proceed similar to [6]. The Lipschitz continuity of a, l, f follow immediately by linearity. We thus have that, for µ1, µ2 ∈ P and u, v ∈ V ,

|a(u, v; µ1) − a(u, v; µ2)| ≤ LakukVkvkV|µ1− µ2|

|f (v; µ1) − f (v; µ2)| ≤ LfkvkV|µ1− µ2|,

|l(u; µ₁) − l(u; µ2)| ≤ LlkukV|µ1− µ2|.

To show Lipschitz continuity of the solution u(µ) note first that for µ1, µ2∈ P,

a(u(µ1), v; µ1) = f (v; µ1), ∀v ∈ V,

a(u(µ2), v; µ2) = f (v; µ2), ∀v ∈ V.

Thus, by linearity of a it follows that

a(u(µ1) − u(µ2), v; µ1) = f (v; µ1) − f (v; µ2) + a(u(µ2), v; µ2) − a(u(µ2), v; µ1), (2.1)

for all v ∈ V . Considering v = u(µ1) − u(µ2) in (2.1) and by Lipschitz continuity of a

and f , we see that

|a(u(µ₁) − u(µ2), u(µ1) − u(µ2); µ1)| ≤ |f (u(µ1) − u(µ2); µ1) − f (u(µ1) − u(µ2); µ2)|

+ |a(u(µ2), u(µ1) − u(µ2); µ2) − a(u(µ2), u(µ1) − u(µ2); µ1)|

≤ L_fku(µ1) − u(µ2)kV|µ1− µ2|

+ Laku(µ2)kVku(µ1) − u(µ2)kV|µ1− µ2|

≤ (Lf + Laku(µ2kV)ku(µ1) − u(µ2)kV|µ1− µ2|.

By uniform coercivity of a and the fact that ku(µ2)kV ≤ γ¯_α_¯f, we get that

¯

αku(µ1) − u(µ2)k2V ≤ (Lf + Laku(µ2)kV)ku(µ1) − u(µ2)kV|µ1− µ2|

≤ (Lf + La

¯ γf

¯

α)ku(µ1) − u(µ2)kV|µ1− µ2|.

Dividing both sides by ku(µ1) − u(µ2)kV and ¯α shows Lipschitz continuity of u. By

linearity and Lipschitz continuity of l, it now immediately follows that the output s is also Lipschitz continuous.

As was noted in [17], also differentiability of the solution function with respect to µ can be shown.

Proposition 2.6. Let a, l, f be such that the coefficient functions θa_q, θ_ql, θqf are

differ-entiable with respect to µ. It then follows that the solution u(µ) is differdiffer-entiable with respect to µ. Furthermore, the partial derivative ∂µiu(µ) ∈ V satisfies the sensitivity

problem

a(∂µiu(µ), v; µ) = ˜fi(v; u(µ), µ)

where ˜fi(·; u(µ), µ) is defined as

˜ fi(·; u(µ), µ) = Qf X q=1 (∂µiθ f q(µ))fq(·) − Qa X q=1 (∂µiθ a q(µ))aq(u(µ), ·; µ).

(20)

It can also be shown that similar statements hold for higher order differentiability. From this we can conclude that smoothness of the coefficient functions of a, l, f is in-herited by the solution u(µ), and thus by the solution manifold M. The smoother the manifold, the better it can be expected to be approximated by a low-dimensional space. Thus, if our Problem 2.1 has (bi)linear forms with smooth coefficient functions, the real solution u(µ) will also depend smooth on µ. This now gives hope that the Reduced Basis method can result in very good approximations.

2.1.3 Thermal block example

We now give an example of the problem stated above. We will consider a specific case of the Thermal Block example from [13]. In this example, we assume to have some solid block consisting of four subblocks with each a different heat conductitivy. We want to solve a parametric elliptic PDE that models the heat transport through this block. The values of the conductivities are the parameters in the problem. The block is heated on one part of the boundary, while it is insulated on two other parts of the boundary and cooled to a chosen temperature on the boundary part that remains. In the end, the output in which we are interested is the average temperature on the boundary part that is heated.

For a precise formulation of the above example, let Ω = (0, 1)2 _{and divide Ω in 4}

subblocks of the same size. We denote these subblocks as Ωi, i = 1, . . . , 4, see Figure 2.1.

On the bottom boundary part the block is heated. It has unit outward normal denoted by n(x) and here there is a unit flux into the domain. We denote this boundary part as ΓN,1. The right and left boundary are the boundary parts that are insulated. Here

no flux will go into the domain and these boundaries are denoted by ΓN,0. The upper

boundary is a homogeneous Dirichlet boundary denoted by ΓD. Here the temperature

is cooled to a chosen temperature. For now, we assign 0 as temperature value on ΓD.

(0,0) ΓN,1 (1,0) (1,1) (0,1) Ω1 Ω2 Ω3 Ω4 ΓN,0 ΓN,0 ΓD

Figure 2.1: Thermal Block Example

The heat conductivityies in the block are denoted by the parameters (µ1, µ2, µ3, µ4) ∈

(21)

k : Ω × µ → R as k(x; µ) := 4 X q=1 µqχΩq(x).

The parametric elliptic PDE we now want to solve is the following: find u(x; µ) such that

−∇ · (k(x; µ)∇u(x; µ)) = 0, x ∈ Ω, u(x; µ) = 0, x ∈ ΓD,

(k(x; µ)∇u(x; µ)) · n(x) = i, x ∈ ΓN,i, i = 0, 1.

The weak formulation of this problem is then naturally given by: find u(·; µ) ∈ H_Γ1

D(Ω) := {u ∈ H1_{(Ω) : u|} ΓD = 0} such that p X q=1 Z Ωq µq∇u(x; µ) · ∇v(x)dx = Z ΓN,1 v(x)dx, (2.2) for all v ∈ H_Γ1

D(Ω). The output in which we now are interested is the average temperature

of the solution on the bottom boundary part, i.e., s(µ) :=

Z

ΓN,1

u(x; µ)dx.

Note that the weak formulation (2.2) of the Thermal Block problem fits exactly in the framework of the Problem 2.1 with H = L2(Ω) and V = H_Γ1_D(Ω). Also, defining the

bilinear form a as a(u, v; µ) := p X q=1 Z Ωq µq∇u(x; µ) · ∇v(x)dx

and the linear forms f, l as

f (v; µ) := Z γN,1 v(x)dx l(u; µ) := Z ΓN,1 u(x; µ)dx,

it is easy to see that a, f, l are parameter separable. Also, it is not hard to show that a is uniformly coercive and continuous and l, f are uniformly continuous.

If we require a lot of solutions for this Thermal block problem, which for example can be the case in an optimization problem for the output functional, it could be inconve-nient to obtain all these solutions with for example the Finite Element method. The Thermal block example is thus a good example where RB-methods can be of use: if we can find a low-dimensional space VN ⊂ V spanned by a couple solutions obtained by

(22)

hope that we can find a good approximate solution for a new µ in this low-dimensional space VN. This will result in solving a low-dimensional system, which is cheap to solve.

As seen in the end of the last section, it can be expected that if the coefficient functions are Lipschitz continuous/differentiable, the solution manifold inherits these properties. In the case of the Thermal Block example, we see that the coefficient functions are of the form µ 7→ µi, which are simply linear functions. These are Lipschitz continuous and

also smooth. This means that the solution manifold depends smoothly on µ and thus we expect that M can be very well approximated by a low-dimensional space.

Remark 2.7. In the next sections, we often want to use or say something about a real solution u(µ) of the Problem 2.1. From a computational point of view, it is not possible to find such a solution u(µ). Instead, we search for an approximate solution which lies very close to the real solution. We thus search for a solution uH(µ) in some

very high dimensional Finite Element space VH ⊂ V . This means that we instead solve

the following problem:

Problem 2.8. For µ ∈ P, find a solution uH(µ) ∈ VH such that

a(uH(µ), v; µ) = f (v; µ), v ∈ VH,

We will call a solution uH(µ) of Problem 2.8 the detailed solution. Note that this problem

also has a unique solution and satisfies the same bounds as in Prop. 2.2, since continuity and coercivity are inherited by the subspace VH.

2.2 Reduced Basis approach

We now want to find an approximate solution of the Problem 2.1 with the Reduced Basis method. For this, we assume to have a low dimensional space spanned by detailed solutions, the so called snapshots, given by

VN := span{uH(µ(1)), . . . , uH(µ(N ))} = span(ΦN) ⊂ VH.

For now, we assume that {uH(µ(i))}Ni=1 are linearly independent and thus form a basis

for VN. The RB-formulation can then be given by

Problem 2.9. For µ ∈ P, find a solution uN(µ) ∈ VN such that

a(uN(µ), v; µ) = f (v; µ), v ∈ VN.

The output sN(µ) ∈ R can then be given by

sN(µ) = l(uN(µ)).

Again, since continuity and coercivity are inherited by the subspace VN, well-posedness

and stability of the solution uN follow from Lax-Milgram with the same constants as

(23)

Proposition 2.10. The Problem 2.9 admits a unique solution uN(µ). Furthermore, the

solution and output sN(µ) satisfy

ku_N(µ)kV ≤ kf (·; µ)k_V0 H α(µ) ≤ ¯ γf ¯ α, |s_N(µ)| ≤ kl(·; µ)kV0 HkuN(µ)kV ≤ ¯ γl¯γf ¯ α .

From now on we will refer to the solution uN(µ) of Problem 2.9 as the reduced solution

or the RB-solution. Further, we will call the space VN the reduced space or RB-space.

The discrete formulation of Problem 2.9 can now be given, by linearity, as the following system of linear equations that needs to be solved.

Proposition 2.11. Let µ ∈ P and a reduced basis ΦN = {φ1, . . . , φN} be given. We can

then define the matrix AN(µ), right hand side vector fN(µ) and output vector lN(µ) as

AN(µ) := (a(φj, φi; µ))Ni,j=1∈ RN ×N

fN(µ) := (f (φi; µ))Ni=1∈ RN

lN(µ) := (l(φi; µ))Ni=1∈ RN.

Then, we can solve the following system for uN = (uN,i)Ni=1∈ RN:

AN(µ)uN(µ) = fN(µ),

and finally, the solution uN(µ) and output sN(µ) of Problem 2.9 are obtained by

uN(µ) = N

X

j=1

uN,jφj, sN(µ) = lTN(µ)uN(µ).

Remark 2.12. In addition to the analytical stability in Proposition 2.10, if we use an orthonormal reduced basis with respect to the inner product on V and assuming that the bilinear form a is symmetric, algebraic stability can also be guaranteed. It can then be shown that the condition number of AN is bounded by γ(µ)/α(µ), which is independent

of N , see [13].

2.2.1 RB-approximation errors

The above Remark 2.12 shows that, if an orthonormal basis is chosen, the condition number of AN does not increase with growing N . Furthermore, the RB-matrix AN is

small. If the Reduced Basis method additionally yields accurate approximations, which we will study next, the RB-method will be favorable over simply applying a Finite Element discretization for every µ.

(24)

Remark 2.13. It is important to note that in this section we will look into the error between the detailed solution and the RB-solution, so into

eH := uH(µ) − uN(µ).

In the end, we are of course interested in the error between the real solution and the RB-solution. Since the detailed solution is supposed to be very accurate, we will for now only consider the error between the detailed and RB-solution. This means that we do not take the error between the real solution and the detailed solution into account. We will come back to this in our numerical experiments in Chapter 4.

First of all, the RB-solution is always as good as the best-approximation from VN up

to a constant factor, i.e. for all µ ∈ P we have that ku_H(µ) − uN(µ)kV ≤

γ(µ) α(µ)v∈VinfN

ku_H(µ) − vkV.

This is a direct corollary of Cea’s Lemma 1.8. Note that if the solution manifold is expected to be well approximable by a low-dimensional subspace, the term on the right hand side of the inequality is supposed to be small for well chosen RB-space VN. This

means that also the error between the detailed solution and the reduced basis solution will be small.

For further error analysis of for example a-posteriori bounds, the following error-residual relation is very useful

Proposition 2.14. For µ ∈ P the residual rH(·; µ) ∈ VH0 is defined as

rH(v; µ) := f (v; µ) − a(uN(µ), v; µ), v ∈ VH.

The error, which is given by eH(µ) = uH(µ) − uN(µ) ∈ VH, then satisfies the following

problem

a(eH, v; µ) = rH(v; µ), v ∈ VH.

Proof. Writing out the definition of a(eH, v; µ) gives us immediately that

a(eH, v; µ) = a(uH − uN, v, µ) = a(uH, v; µ) − a(uN, v; µ)

= f (v) − a(uN, v; µ) = rH(v; µ),

for all v ∈ VH.

Another important property of the RB-solution is the property of reproduction of solutions. This property says that if a detailed solution lies in the space, the RB-scheme of Problem 2.9 will give this detailed solution as the reduced solution.

Proposition 2.15. If uH(µ) ∈ VN for some µ ∈ P, then uN(µ) = uH(µ). This means

(25)

Proof. Since uN(µ) ∈ VN, we have that if uH(µ) ∈ VN, also eH ∈ VN. Now, by coercivity it follows that α(µ)keHk2V ≤ a(eH, eH; µ) = a(uH(µ), eH; µ) − a(uN(µ), eH; µ) = f (eH; µ) − f (eH; µ) = 0, hence eH = 0.

This property can be used, among other things, to show uniform convergence of the RB-solution to the real solution. For this, first assume that the parameter-space P is compact and let ΦN be a snapshot basis, so φi= uH(µ(i)), i = 1, . . . , N . Denote further

the parameter set belonging to the snapshot basis by SN, i.e. SN = {(µ(i))}Ni=1. If it

now holds that a, l, f are Lipschitz continuous, it is easy to show, similarly to Prop. 2.5, that the RB-solution uN : P → VN is Lipschitz continuous with Lipschitz constant Lu,

which is independent of N . Then, it holds for all N, µ and µ∗ = arg min_µ0_∈S

Nkµ − µ 0_k Rp, that kuH(µ) − uN(µ)kV ≤ kuH(µ) − uH(µ∗)kV + kuH(µ∗) − uN(µ∗)kV + kuN(µ) − uN(µ∗)kV = Lukµ − µ∗kRp+ 0 + Lukµ − µ ∗_k Rp = 2Lukµ − µ∗kRp,

where the second to last equality holds because of Prop. 2.15. Suppose now that the sets SN get dense in P. Thus, if we define hN := sup dist(µ, SN), we assume that

lim

N →∞hN = 0.

This now means that, from the above,

kuH(µ) − uN(µ)kV ≤ 2hNLu,

and this gives us that

lim

N →∞_µ∈PsupkuH(µ) − uN(µ)kV = 0.

Thus uniform convergence of the RB-solution to the detailed solution holds. This con-vergence rate, which is linear in hN, is not always of real practical value: hN might decay

really slow with N and then the snapshot basis ΦN must become very large to guarantee

small errors. Later on, in Section 2.3.1, we see that if we choose our parameters µ(i) for the snapshot basis carefully, this can result in exponential convergence.

We next investigate an a-posteriori error bound and an effectivity index bound, which are important for the certification of the RB-methods. For this, we first of all assume that we can compute a lower bound for the coercivity constant, which we will denote by αLB. We should be able to compute this bound rapidly. Also, this constant is assumed

to be still relatively large, so ¯α ≤ αLB(µ). For the computational aspects of determining

(26)

Proposition 2.16. Let αLB be the computational lower bound as given above. For all

µ ∈ P, the error eH(µ) then satisfies

ku_H(µ) − uN(µ)kV ≤ ∆u(µ) := kr_H(·; µ)k_V0 H αLB(µ) , |sH(µ) − sN(µ)| ≤ ∆s(µ) := kl(·; µ)kV0 H∆u(µ).

Proof. By coercivity and the error-residual relation from Prop. 2.14 we obtain that αLB(µ)keHk2V ≤ α(µ)keHk2V ≤ a(eH, eH; µ) = rH(eH; µ) ≤ krH(·; µ)kV_H0 keHkV.

Dividing both sides by keHkV and αLB gives the bound for keHkV. The bound for the

output now follows immediately by continuity of the functional l. Since we assumed that VH is discrete space, the norm krHkV0

H is now a computable

quantity. In the analysis of Finite Element Methods, the residual norm krkV0 is not

available, since this norm is then computed over an infinite dimensional space. Since we can compute the residual norm (and thus the error bound in Prop. 2.16) after the computation of the RB-solution uN(µ), this error bound is called an a-posteriori error

bound for the RB-solution. Also, since this bound is a proven upper bound for the error, this error is a rigorous error bound. This now implies that the RB-method is a certified method: together with the RB-solution, we can simultaneously obtain a rigorous error bound, which immediately gives us the certification of the RB-method.

Next, we will study how tight this a-posteriori error bound actually is. For this, we first show that if the error is zero, the error bound also is zero.

Proposition 2.17. If the error between the reduced solution and the detailed solution is zero, i.e. if uH(µ) = uN(µ), then ∆u(µ) = ∆s(µ) = 0.

Proof. Since 0 = a(0, v; µ) = a(eH, v; µ) = rH(v; µ) for all v ∈ VH, we have that

krH(·; µ)kV_H0 = sup v∈VH

|r_H(v; µ)| kvkV

= 0.

Hence ∆u(µ) = 0 and also ∆s(µ) = 0.

Next, we show that the factor of overestimation for the a-posteriori error bound, also called the effectivity index, can be bounded by a constant independent of the parameter µ. Note that this also gives us a lower bound for the error.

Proposition 2.18. The effectivity index ηu(µ) is defined as and bounded by

ηu(µ) := ∆u(µ) ku_H(µ) − uN(µ)kV ≤ γ(µ) α_LB(µ) ≤ ¯ γ ¯ α,

where γ(µ) is the continuity constant and αLB is the computable lower bound for the

(27)

Proof. Let vr ∈ VH be the Riesz-representation of the residual rH(·; µ), i.e.

hv_r, viV = rH(v; µ), v ∈ VH, kvrkV = krH(·; µ)kV_H0.

Then, because of the error-residual relation from Prop. 2.14 and the continuity of a, it follows that

kvrk2V = hvr, vriV = rH(vr; µ) = a(eH, vr; µ) ≤ γ(µ)keHkV0

HkvrkV.

Dividing both sides with kvrkV and keHkV0

H we obtain that

kvrkV

keHkV0 h

≤ γ(µ).

Hence we conclude that ηu(µ) = ∆u(µ) ke_Hk_V0 H = kvrkV αLB(µ)keHkV_H0 ≤ γ(µ) αLB(µ) .

The parameter-independent bound now follows immediately from the uniform coercivity and continuity.

Example 2.19. Going back to the Thermal Block Example from Section 2.1.3, we can give some numerical experiments regarding the error and effectivity index bounds. The results are reproductions from the package RBmatlab [18], using the script rb_tutorial.m, and are the same as in [13].

For this, consider the Thermal Block example with µmin = 1/µmax= 0.1 and construct

a basis ΦN = {uH(µ(i))}5i=1 with µ(i) = (0.1 + 0.5(i − 1), 0.1, 0.1, 0.1)T, i = 1, . . . , 5. We

are then going to look for RB-approximations for µ = (µ1, 0.1, 0.1, 0.1), µ1 ∈ [0.1, 4], in

the RB-space with basis ΦN.

First, a plot of the a-priori error-estimator together with the true error is given in Figure 2.2. We see that the error estimator is indeed an upper bound for the true error. Also, in the sampling points µ(i), the error kuH(µ(i)) − uN(µ(i))kV is (numerically) zero.

This is exactly what we should see according to Prop. 2.15. Further, also the error-estimator is (numerically) zero in these point, which we also expected according to Prop. 2.17. Note that the true error is given as a dotted line, since evaluation of the real error is very expensive. The error-estimator, as we will see later, is not expensive to evaluate. Therefore the line of the estimator is finely resolved.

Secondly a plot of the effectivity index and the bound for the effectivity index is given, see Figure 2.3. Note that again the effectivity index is a dotted line, while the effectivity bound is a finely resolved line. We see that indeed, as in Prop. 2.18, the effectivity index bound is an upper bound for the effectivity index. This picture shows that the effectivity index is of order 10. This is considered as quite a good effectivity index of the error bound.

(28)

Figure 2.2: Error and error bounds

Figure 2.3: Effectivity and effectivity bound

2.3 Offline/Online decomposition

In the next section we will give the details on the computational aspects of the RB-method: the so called offline and online decomposition. For this, we first introduce some assumptions and notation.

Assumptions and notation. As noted in Remark 2.7 we take as detailed solution the approximate solution of the Problem 2.8, which will be denoted as uH(µ). We therefore

assume that VH is some very high dimensional Finite Element space with respect to a

uniform subdivision of some domain into triangles or rectangles. Further, let the basis of VH be given by {ψi}Hi=1, where ψi are the basis functions and H is some very large

number. The matrix A(µ), inner product matrix K and vectors f and l are then given by

A(µ) := (a(ψj, ψi; µ))Hi,j=1∈ RH×H, K := (hψi, ψjiV)Hi,j=1∈ RH×H,

f (µ) := (f (ψi; µ))Hi=1∈ RH, l(µ) := (l(ψi; µ))Hi=1∈ RH.

We can now obtain a detailed solution uH(µ) =PHi=1uH,i(µ)ψi by solving the system

(29)

where (uH)i= uH,i. The output solution is now obtained by s(µ) = l(µ)TuH(µ). As in

most discretizations of PDEs with basis functions of local support, we have that A(µ) is a sparse matrix.

From the above it follows that the complexity of obtaining a detailed solution is O(H2_{), while a reduced solution requires O(N}3_{) operations, since the resulting N × N}

-matrix from Prop. 2.11 for a reduced solution is a dense -matrix. This means that the RB-approach only makes sense if N H.

Of course, before we can compute a reduced solution, we have to construct a reduced basis space VN spanned by detailed solutions for different parameters {µi}. We will skip

the construction of the reduced basis for now, it will be discussed in Section 2.3.1. To make sure that we can compute a reduced solution efficiently, we want to decompose the computation of such a reduced solution into an Offline and Online phase.

• Offline Phase During the offline phase, which is the expensive phase, first the RB-basis is constructed by computing N detailed solutions. After, high dimensional quantities are precomputed which are µ independent. The operation count here depends on H.

• Online Phase During the online phase, which is the cheap phase, the precomputed data from the offline phase is used to obtain a µ dependent reduced system and from this system, the solution uN(µ) and sN(µ) can now be computed rapidly.

Here the operation count should depend on N and not on H.

To see how such an Offline/Online decomposition can be established, we first look into the complete computation of an RB-solution. The complete computation of such an RB-solution can be divided into the following four steps

1. N computations of the detailed Problem 2.8 for the construction of the Reduced Basis consisting of snapshots : O(N H2)

2. N2 evaluations of a(φj, φi; µ) for the construction of the matrix AN(µ): O(N2H)

3. N evaluations of f (φi; µ) for the construction of the vector fN(µ): O(N H)

4. Obtaining a solution of the N × N system in the reduced Problem 2.9: O(N3). We now want to divide the above steps into the offline and online phase. Clearly, step 1 belongs to the offline phase and step 4 to the online phase. It is not so clear where to put step 2 and 3, since in those steps we require both expensive as well as parameter-dependent computations. Here the parameter-separability (see Def. 2.4) of a, l, f comes in. From linearity it follows that this parameter-separability transfers to parameter-separability of the matrix and vectors AN, fN, lN. This now assures a proper

(30)

Proposition 2.20. (Offline Phase) After the construction of a reduced basis ΦN =

{φ1, . . . , φN} the parameter-independent component matrices and vectors can be

con-tructed in the following way:

AN,q:= (aq(φj, φi))Ni,j=1 ∈ RN ×N, q = 1, . . . , QA,

fN,q:= (fq(φi))Ni=1∈ RN, q = 1, . . . , Qf,

lN,q:= (lq(φi))Ni=1∈ RN, q = 1, . . . , Ql.

(Online Phase) For a given µ ∈ P, the coefficient functions θa_q(µ), θqf(µ) and θlq(µ)

can be evaluated for q in suitable ranges and the matrix and vectors can be assembled as follows AN(µ) = Qa X q=1 θa_q(µ)AN,q, fN(µ) = Qf X q=1 θf_q(µ)fN,q, lN(µ) = Ql X q=1 θl_q(µ)lN,q.

This results in the discrete system from Prop. 2.11 for which its solution uN(µ) and the

functional sN(µ) of the solution can be solved.

We see here that indeed the online phase does not depend on N , since its operation count is O(N2· Qa+ N · (Qf+ Ql) + N3).

Remark 2.21. The reduced component matrices AN,q and vectors fN,q, lN,q in the

Of-fline phase can be obtained in a very simple way. For this, write the reduced basis vector φj as φj =PHi=1φi,jψi, i.e. the expansion in the basis of our very high dimensional

dis-crete space VH. Let

ΦN = (φi,j)H,N_i,j=1∈ RH×N

be the coefficient matrix. If the component matrices Aq = (aq(ψj, ψi))Hi,j=1 and vectors

fq= (fq(ψi))H_i=1, lq = (lq(ψi))H_i=1of the Finite Element space for the detailed solution are

available, the reduced component matrices and vectors can be obtained by evaluating AN,q= ΦTNAqΦN, fN,q = ΦTNfq, lN,q = ΦTNlq.

As noted in Section 2.2.1, the certification of the RB-method is a very important and usefull topic. To ensure that the a-posteriori error estimators can also be calculated efficiently, we would like to Offline/Online decompose the residual norm, which is the important quanitity in the error estimators. This can be done since, by linearity, the parameter-separability also holds for the residual and its Riesz-representative. Remem-ber that the Riesz-representative vr ∈ V of a functional r ∈ V0 is defined by

r(v) = hv, vriV, ∀v ∈ V, krkV0 = kv_rk_V.

Proposition 2.22. Set Qr:= Qf + N Qa and define rH,q ∈ V_H0 , q = 1, . . . , Qr as

(31)

Let further uN(µ) =PNi=1uN,iφi be the RB-solution and define θrq(µ), q = 1, . . . , Qr as (θr₁, . . . , θ_Qr_r) := (θ₁f, . . . , θf_Q f, −θ a 1· uN,1, . . . , −θaQa· uN,1, . . . , −θ a 1· uN,N, . . . , −θaQa· uN,N).

Let vr, vr,q ∈ VH be the Riesz-representations of respectively rH and rH,q. Then we can

write the residual and its Riesz-representation as follows

rH(v, µ) = Qr X q=1 θr_q(µ)rH,q(v), vr(µ) = Qr X q=1 θr_q(µ)vr,q, µ ∈ P, v ∈ VH.

Proof. This follows directly from linearity of the bilinear form a and the Riesz-map. Since the residual is parameter-separable by the above, so it its norm. We can thus decompose the residual norm into an offline/online decomposition:

Proposition 2.23. (Offline Phase) After the offline-phase of Prop. 2.20, define the matrix Gr ∈ RQr×Qr as

Gr := (rq(vr,q0))Qr

q,q0₌₁ = (hvr,q, vr,q0i_V)Qr

q,q0₌₁.

(Online Phase) For given µ ∈ P, after the computation of the RB-solution uN(µ),

compute the coefficient vector θr(µ) := (θ1r(µ), . . . , θQrr(µ))

T_{. The residual norm can}

then be obtained by

kr(·; µ)k_V0 H =

q

θr(µ)TGrθr(µ).

Proof. Since the Riesz-map is an isometry and the Riesz-representative is parameter seperable, see Prop. 2.22, this gives us

kr(µ)k2_V0 H = kvr(µ)k 2 V = Qr X q=1 θ_qr(µ)vr,q, Qr X q0₌₁ θr_q0(µ)v_r,q0 V = θr(µ) T_G rθr(µ).

Note that in the Offline/Online decomposition of the residual norm, we need to cal-culate the Riesz-representatives vr,q of rq once in the Offline phase. As is shown in [13],

Lemma 2.31, this can be done by solving a sparse H ×H system, where the inner product matrix K (see assumptions and notations) has to be inverted.

Example 2.24. (Thermal Block) Going back to the Thermal Block Example 2.1.3, we are going to look into the offline/online decomposition. Assume that we have a reduced basis consisting of snapshots,

ΦN = {uH(µ(1)), . . . , uH(µ(N ))} =: {φ1, . . . , φN} ⊂ HΓ1D(Ω),

where the snapshots are calculated for some discretization (e.g. using FEM) of (2.2). We are then looking for a solution uN(µ) ∈ VN := span(ΦN) that satisfies

p X q=1 Z Ωq µq∇uN(x; µ) · ∇vN(x)dx = Z ΓN,1 vN(x)dx, (2.3)

(32)

for all vN ∈ VN.

In the Offline phase, we evaluate the component matrices AN,q as

(AN,q)i,j =

Z

Ωq

∇φj· ∇φidx, q = 1, . . . , 4,

and the component vectors fN and lN as

(fN)i= (lN)i=

Z

ΓN,1

φidx, q = 1, . . . , 4.

For a given µ ∈ P, a reduced solution uN(µ) can now be obtained by solving the

matrix-vector equation AN(µ)uN(µ) = fN, where

AN(µ) = 4

X

q=1

µqAN,q.

The output is then given by lT_NuN(µ).

2.3.1 Basis Generation

An important part of the RB-method is the generation of the basis of the RB-space VN. As mentioned before, this is a part of the offline phase, but since it is so important

we consider this separately. In this subsection, we will discuss the basis generation via the Greedy algorithm, but for example also an algorithm based on the POD, see also [5], can be used. In Chapter 3, the POD will be discussed as a part of the so called POD-Greedy algorithm for solving parabolic PDEs with the RB-method. The differences between the POD and the Greedy algorithm will also be discussed there, see Remark 3.8. One of the most simple reduced basis types, the one we also used in the previous sections, is the so called Lagrangian Reduced Basis.

Definition 2.25. Let the snapshots {uH(µ(i))} be linearly independent. We call

ΦN := {uH(µ(1)), . . . , uH(µ(N ))} ⊂ VH

a Lagrangian reduced basis.

The Lagrangian reduced basis can provide globally well approximations of the manifold M if the set of parameters SN := {µ(i)} for the snapshot basis is carefully chosen. We

are aiming for a globally accurate approximation of the manifold M. This means we are interested in finding a space VN with dim(VN) = N minimizing

d(M, VN) := sup µ∈P inf v∈VN ku_H(µ) − vkV =: sup µ∈P δ(uH(µ), VN). (2.4)

(33)

Note that we are actually interested in minimizing the maximal error between the de-tailed solution and the RB-solution, so, into minimizing

sup

µ∈M

kuH(µ) − uN(µ)kV. (2.5)

Cea’s Lemma 1.8 tells us that, up to a constant factor, (2.5) is majorized by (2.4). This means that minimizing the maximal projection error will also result in approximately minimizing the maximal error in the RB-solution. For now, we will consider minimizing (2.4), since this gives the more general idea of the problem.

We want to find a space VN of dimension N that minimizes the maximal projection

error δ(uH, VN). Of course, if the solution manifold M cannot be well approximated by a

low-dimensional linear subspace, we are also not going to find such a space. The so called Kolmogorov N -width gives us a quantity that tells us how well M can be approximated by some N -dimensional linear subspace. The Kolmogorov N -width is defined as

dN(M) := inf Y ⊂V, dim Y =N sup u∈M ku − P_YukV.

Thus, if the Kolmogorov N -width of M is small, we hope to find a space VN that is

(almost) as good as the optimal subspace.

Finding such a subspace VN is a very complex problem which cannot be solved

analyt-ically. We therefore introduce an algorithm that constructs an approximating subspace by adding new basis vectors. These new basis vectors are chosen in a way that aims for minimization of d(M, VN). Aiming for minimizing d(M, VN) is done by considering an

error indicator ∆(Y, µ) ∈ R+, which is an estimator for δ(uH(µ), VN) when VN = Y is the

current approximation space. The algorithm below is known as the Greedy-algorithm. Algorithm 1. (Greedy algorithm) Let Strain ⊂ P be a training set of parameters

and let tol > 0 be a given error tolerance. Set V0 := {0}, S0 := ∅, Φ0 := ∅ and define

iteratively while N := max µ∈Strain ∆(VN, µ) > tol µ(N +1) := arg max µ∈Strain ∆(VN, µ) SN +1:= SN ∪ {µ(N +1)} φN +1:= uH(µ(N +1)) ΦN +1:= ΦN∪ {φN +1} VN +1:= VN + span(φN +1).

To ensure termination of this algorithm we require that for any subspace Y ⊂ VH it

holds that

(34)

In this case, the algorithm terminates in at most |Strain| steps. The training set Strain

should be chosen such that it represents the parameter set P well: it should not miss any important parts of P. Note that since the error sequence {n} is generated on the

training set, it will happen that the error on the whole parameter set is not as small as on the training set. Later in this section we will comment on the choice of this training set.

The error-indicator can be chosen in different ways. Three options are the following; 1. The first option is the projection error,

∆(Y, µ) := inf

v∈Y kuH(µ) − vkV = kuH(µ) − PYuH(µ)kV,

which is exactly the quantity δ(uH(µ), VN) that we want to estimate. Using this

indicator, the greedy algorithm will simply construct a good approximation space, without taking the RB-model into account. This is thus exactly the general prob-lem we were looking at. This error-indicator is very expensive to evaluate: there are high dimensional operations required to evaluate the projection error. This means that the size of our training set Strain cannot be too big. Also, we have

to store all the snapshots uH(µ), which will limit the size of the training-set due

to memory contraints. When using this indicator in the greedy algorithm, the algorithm is called the strong greedy procedure.

2. Another indicator, which is similar to the projection error as indicator, is the error between the detailed solution and the RB-solution:

∆(Y, µ) := kuH(µ) − uN(µ)k.

For this indicator, we need to have an RB-model available. As in the projection error indicator, the error-indicator is expensive to evaluate and we have to store all the snapshots uH(µ), which again limits the size of Strain. The advantage is, of

course, that this is actually the quantity (2.5) we want to minimize. 3. The last option discussed here is the a-posteriori error:

∆(Y, µ) := ∆u(µ) =

krH(·; µ)kV0 H

αLB

.

For this estimator we both need an RB-model and an a-posteriori error estimator. The evaluation of this indicator is very cheap. This is due to the fact that the residual norm krH(µ)kV_H0 can be decomposed in an Offline/Online decomposition.

Thus, for evaluation of this indicator, in each loop of the Greedy algorithm, we need to precompute some expensive quantities once, and can then quickly evaluate the residual norms krH(µ)kV_H0 for all µ in the training set. Since this indicator is

cheap to evaluate, our training set Strain can be chosen a lot bigger than in the

previous cases. This means that the training set will be more representative for the whole parameter space. This version of the greedy algorithm is called the weak greedy algorithm.

(35)

Note that all three indicators satisfy the termination criterium 2.6. For the projection error, this is a trivial statement. The statement for the true RB-error follows from Prop. 2.15 and for the a-posteriori error it follows from Prop. 2.17.

The theoretical foundation for the greedy algorithm can be seen in the following propo-sition, which was formulated and proved in [2]. It states that if the Kolmogorov N -width decays polynomially or exponentially as a function of N , then so does N, and so does

d(M, VN) with VN produced in the Greedy algorithm.

Proposition 2.26. Let Strain = P be compact and let the error indicator ∆ be chosen

such that for some η ∈ (0, 1] and the sequence {µ(i)} generated by the Greedy Algorithm 1 it holds that ku(µ(N +1)) − PVNu(µ (N +1)_{)k ≥ η sup} u∈M ku − PVNuk. (2.7) Then

1. if for some α, M > 0 and all N ∈ N, it holds that dN(M) ≤ M N−α, and d0(M) ≤

M ,then

N ≤ CM N−α, n > 0

with suitable constant C > 0.

2. if for N ≥ 0 and M, a, α > 0, it holds that dN(M) ≤ M e−aN

α

, then N ≤ CM e−cN

β

, N ≥ 0 with β := α/(α + 1) and suitable constants c, C > 0,

Note that the projection error indicator obviously meets the requirement (2.7) with η = 1 and the error between the detailed solution and RB-solution meets the requirement as a consequence of Cea’s Lemma 1.8, with η = α_¯_γ¯. The a-posteriori error meets the requirement (2.7) due to again Cea’s Lemma 1.8, the effectivity from Prop. 2.18 and the error bound Prop. 2.16. It can be shown, see [13], Remark 2.42, that in this case η = α_¯_γ¯22.

Remark 2.27. A lot can be said about the choice of the training set Strain. There are

many ways of choosing this but since we will always choose the training set as a subset of P with equidistance or random points in our numerical experiments, we won’t go into the details here. For an overview of some choices for a training set, see [13], Remark 2.44. Note that it is always a good idea to check the error by calculating the error over a test training set consisting of random chosen parameters.

Example 2.28. We will end this chapter by looking one last time at the Thermal Block example 2.1.3. Again, we will show a reproduction of the numerical experiments shown in [13] and used for this the package RBmatlab [18], rb_tutorial.m.

In Figure 2.4 the decay of the error is plotted. For this, the parameter set P is chosen as [µmin, µmax] with µmin = 1/µmax = 0.5, the training set consists of 1000 random

(36)

Figure 2.4: Decay error indicator and test error for Greedy Algorithm

points, the tolerance tol is set to 10−6 and the error indicator for the greedy algorithm

is chosen as the a-posteriori error bound. The plot now depicts the maximal error estimator and the maximal error over a test training set of size 100, for the size of the reduced basis N ranging from 1 to 20. We see that plot confirms the statement in Prop. 2.26: the RB-solution converges exponentially.

Reduced basis methods for a simultaneous space-time variational formulation of parametric parabolic PDEs

MSc Mathematics

Master Thesis