# Finite Element Method

An Introduction to the Finite Element Method (FEM) for Di?erential Equations

Mohammad Asadzadeh January 20, 2010

Contents

0 Introduction 0.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.2 Trinities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 5 6

1 Polynomial approximation in 1d 15 1.1 Overture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.2 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2 Polynomial Interpolation in 1d 2.1 Preliminaries . . . . . . . . . . . . . . 2.2 Lagrange interpolation . . . . . . . . . 2.3 Numerical integration, Quadrature rule 2.3.1 Gauss quadrature rule . . . . . 2.4 Exercises . . . . . . . . . . . . . . . . . 41 41 48 52 58 62

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

3 Linear System of Equations 65 3.1 Direct methods . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.2 Iterative method . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4 Two-points BVPs 4.1 A Dirichlet problem . . . . . . . . 4.2 Minimization problem . . . . . . . 4.3 A mixed Boundary Value Problem 4.4 The ?nite element method. (FEM) 4.5 Error estimates in the energy norm 4.6 Exercises . . . . . . . . . . . . . . . 3 83 83 86 87 90 91 97

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

4 5 Scalar Initial Value Problem 5.1 Fundamental solution and stability . . . . 5.2 Galerkin ?nite element methods (FEM) for 5.3 An a posteriori error estimate for cG(1) . 5.4 A dG(0) a posteriori error estimate . . . . 5.5 A priori error analysis . . . . . . . . . . . 5.6 The parabolic case (a(t) ≥ 0) . . . . . . . 5.6.1 Short summary of error estimates . 5.7 Exercises . . . . . . . . . . . . . . . . . . . 6 The 6.1 6.2 6.3 6.4 heat equation in 1d Stability estimates . . . . FEM for the heat equation Error analysis . . . . . . . Exercises . . . . . . . . . .

CONTENTS 103 . 103 . 105 . 108 . 114 . 116 . 120 . 124 . 127 129 . 130 . 135 . 141 . 146

. . . IVP . . . . . . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

7 The wave equation in 1d 149 7.1 FEM for the wave equation . . . . . . . . . . . . . . . . . . . 150 7.2 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 8 Piecewise polynomials in several dimensions 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Piecewise linear approximation in 2 D . . . . . . . . . . . 8.2.1 Basis functions for the piecewise linears in 2 D . . 8.2.2 Error estimates for piecewise linear interpolation . 8.2.3 The L2 projection . . . . . . . . . . . . . . . . . . 8.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 157 . 157 . 160 . 160 . 164 . 165 . 166

. . . . . .

. . . . . .

Riesz and Lax-Milgram Theorems 169 9.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 9.2 Riesz and Lax-Milgram Theorems . . . . . . . . . . . . . . . . 174 9.3 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 Poisson Equation 181 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 Error Estimates for FEM . . . . . . . . . . . . . . . . . . . . . 182 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191

10 The 10.1 10.2 10.3

CONTENTS

5

11 The heat equation in RN 195 11.1 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196 11.2 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 12 The wave equation in RN 203 12.1 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204 13 Convection - di?usion problems 13.1 A convection-di?usion model problem . . . . . . . 13.1.1 Finite Element Method . . . . . . . . . . . 13.1.2 The Streamline - di?usion method (SDM) 13.2 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207 207 210 212 214

6

CONTENTS

Chapter 0 Introduction

This note presents an introduction to the Galerkin ?nite element method (FEM), as a general tool for numerical solution of partial di?erential equations (PDEs). Iteration procedures and interpolation techniques are also employed to derive basic a priori and a posteriori error estimates, necessary for, e.g. solution properties, such as stability and convergence. Galerkin’s method for solving a general di?erential equation (both PDEs and ODEs) is based on seeking an approximate solution, which is 1. easy to di?erentiate and integrate 2. spanned by a set of nearly orthogonal basis functions in a ?nite-dimensional space.

0.1

Preliminaries

? A di?erential equation is a relation between an unknown function u and its derivatives u(k) , 1 ≤ k ≤ N , where k and N are integers. ? If the function u(x) depends on only one variable (x ∈ R), then the equation is called an ordinary di?erential equation, (ODE). ? The order of the di?erential equation is determined by the order of the highest derivative of the function u that appears in the equation. ? If the function u(x, t) depends on more than one variable, and the derivative with respect to more than one variable is present in the equation, then the 7

8

CHAPTER 0. INTRODUCTION

di?erential equation is called a partial di?erential equation, (PDE), e.g.: ut (x, t) ? uxx (x, t) = 0 is a homogeneous PDE of second order, whereas uyy (x, y ) + uxx (x, y ) = f (x, y ), is a non-homogeneous PDE of second order. ? A solution to a di?erential equation is a function; e.g. u(x), u(t, x) or u(x, y ). ? In general the solution u cannot be expressed in terms of elementary functions and numerical methods are the only way to solve the di?erential equation by constructing approximate solutions. Then the main question in here is: how close is the approximate solution to the exact solution? and how and in which environment does one measure this closeness? In which extent the approximate solution preserves the physical quality of the exact solution? These are some of the questions that we want to deal with in this notes. ? The linear di?erential equation of order n in time has the general form: where D = d/dt denotes the ordinary time derivative, and D k = k ≤ n. The corresponding linear di?erential operator is denoted by L(t, D) = dn d n?1 d + a ( t ) + . . . + a1 (t) + a0 (t). n ? 1 n n ? 1 dt dt dt

dk , dtk

L(t, D)u = u(n) (t) + an?1 (t)u(n?1) (t) + . . . + a1 (t)u′ (t) + a0 (t)u(t) = b(t), 1 ≤

0.2

Trinities

To continue we introduce the so called trinities classifying the main ingredients involved in the process of modeling and solving problems in di?erential equations, see Karl E .Gustafson [14] for details. The usual three operators involved in di?erential equations are Laplace operator Di?usion operator D’Alembert operator ?n = ?2 ?2 ?2 + + . . . + , ?x2 ?x2 ?x2 1 2 n (0.2.1) (0.2.2) (0.2.3)

d ? ?n , dt d2 = 2 ? ?n , dt

0.2. TRINITIES

9

where we have the space variable x := (x1 , x2 , x3 , . . . xn ) ∈ Rn , the time variable t ∈ R+ and ? 2 /?x2 i denotes the second partial derivative with respect to xi . We also recall a ?rst order operator: the gradient operator ?n which is a vector valued operator and is de?ned as follows: ?n = ? ? ? , ,..., . ?x1 ?x2 ?xn

Classifying the general second order PDEs in two dimensions The usual three classes of second order partial di?erential equations are elliptic, parabolic and hyperbolic ones. Second order PDEs with constant coe?cients in 2-D: Auxx (x, y )+2Buxy (x, y )+Cuyy (x, y )+Dux(x, y )+Euy (x, y )+F u(x, y )+G = 0. Here we introduce the discriminant d = AC ? B 2 : a quantity that speci?es the role of the coe?cients in determining the equation type. Discriminant Type of equation Solution behavior stationary energy-minimizing smoothing and spreading ?ow a disturbance-preserving wave

d = AC ? B 2 > 0 Elliptic

d = AC ? B 2 = 0 Parabolic

d = AC ? B 2 < 0 Hyperbolic

Example 1. Here are the class of the most common equations: Elliptic Potential equation d2 u d2 u + =0 dx2 dy 2 uxx (x, y ) + uyy (x, y ) = 0 A = C = 1, B = 0 d = AC ? B 2 = 1 > 0 Parabolic Heat equation du ? ?u = 0 dt ut (t, x) ? uxx (t, x) = 0 A = B = 0, C = ?1 d = AC ? B 2 = 0 Hyperbolic Wave Equation d2 u ? ?u = 0 dt2 utt (t, x) ? uxx (t, x) = 0 A = 1, B = 0, C = ?1 d = AC ? B 2 = ?1 < 0.

10

CHAPTER 0. INTRODUCTION

Second order di?erential equations with variable coe?cients in 2-D In the variable coe?cients case, one can only have a local classi?cation. Example 2. Consider the Tricomi operator of gas dynamics: Lu(x, y ) = yuxx + uyy . Here the coe?cient y is not a constant and we have A = y, B = 0, and C = 1. Hence d = AC ? B 2 = y and consequently, e.g. the domain of ellipticity is y > 0, and so on (see the Fig. below) y

elliptic parabolic

x

hyperbolic

Figure 1: Tricomi: an example of a variable coe?cient classi?cation.

?Summing up and generalizing to n space variables we have Mathematical Quantity ?n d ? ?n?1 dt d2 = dt 2 ? ?n?1 Laplacian Heat d’Alembertian Surname Physical Named Classi?cation Type

Potential operator Elliptic Di?usion operator Wave operator Parabolic Hyperbolic

0.2. TRINITIES The usual three types problems in di?erential equations

11

1. Initial value problems (IVP)

The simplest di?erential equation is u′ (x) = f (x) for a < x ≤ b, but also (u(x) + c)′ = f (x) for any constant c. To determine a unique solution a speci?cation of the initial value u(a) = u0 is generally required. For example for f (x) = 2x, 0 < x ≤ 1, we have u′ (x) = 2x and the general solution is u(x) = x2 + c. With an initial value of u(0) = 0 we get u(0) = 02 + c = c = 0. Hence the unique solution to this initial value problem is u(x) = x2 . Likewise for a time dependent di?erential equation of the second order (two time derivatives) the initial values for t = 0, i.e., u(x, 0) and ut (x, 0) are generally required. For a PDE such as the heat equation the initial value can be a function of the space variable. Example 3. The wave equation, on real line, associated with the given initial data: ? ? u ? u = 0, ?∞ < x < ∞ t > 0, tt xx ? u(x, 0) = f (x), u (x, 0) = g (x), ?∞ < x < ∞, t = 0.

t

2. Boundary value problems (BVP)

a. Consider the boundary value problems in R: Example 4. The stationary heat equation: ?[a(x)u′ (x)]′ = f (x), for 0 < x < 1.

To de?ne a solution u(x) uniquely, the di?erential equation is complemented by boundary conditions imposed at the boundaries x = 0 and x = 1: for example u(0) = u0 and u(1) = u1 , where u0 and u1 are given real numbers. b. Boundary value problems (BVP) in Rn : Example 5. The Laplace equation in Rn , x = (x1 , x2 , . . . , xn ): ? 2 2 2 ? ? u = ? u + ? u + . . . + ? u = 0, x ∈ ? ? Rn , n 2 2 ?x2 ?x ?x 1 2 n ? u(x, t) = f (x), x ∈ ? ? (boundary of ?).

12

CHAPTER 0. INTRODUCTION

Remark 1. In general, in order to obtain a unique solution for a (partial) di?erential equation, one should supply as many data as the sum of highest order (partial) derivatives involved in the equation. Thus in example 1, to determine a unique solution for the potential equation uxx + uyy we need to give 2 boundary conditions in the x-direction and another 2 in the y -direction, whereas to determine a unique solution for the wave equation utt ? uxx = 0, it is necessary to supply 2 initial and 2 boundary conditions.

3. Eigenvalue problems (EVP)

Let A be a given matrix. The relation Av = λv, v = 0 is a linear equation system where λ is an eigenvalue and v is an eigenvector. Example 6. In the case of a di?erential equation; e.g. the equation of a steady state vibrating string ?u′′ (x) = λu(x), u(0) = u(π ) = 0,

where λ is an eigenvalue and u(x) is an eigenfunction. u(0) = 0 and u(π ) = 0 are boundary values. The di?erential equation for a time dependent vibrating string with small displacement, and ?xed at the end points, is given by ? ? utt (x, t) ? uxx (x, t) = 0 0<x<π t≥0 ? u(0, t) = u(π, t) = 0, t ≥ 0, u(x, 0) = f (x) u (x, 0) = g (x).

t

Using separation of variables, see also Folland [11], this equation can be split into two eigenvalue problems: Insert u(x, t) = X (x)T (t) = 0 (a nontrivial solution) into the di?erential equation to get utt (x, t) ? uxx (x, t) = X (t)T ′′ (t) ? X ′′ (x)T (t) = 0. Dividing (0.2.4) by X (x)T (t) = 0 separates the variables, viz T ′′ (t) X ′′ (x) = =λ=C T (t) X (x) (independent of x and t ). (0.2.5) (0.2.4)

Consequently we get 2 ordinary di?erential equations (2 eigenvalue problems): X ′′ (x) = λX (x), and T ′′ (t) = λT (t). (0.2.6)

0.2. TRINITIES The usual three types of boundary conditions

13

1. Dirichlet boundary condition (the solution is known at the boundary of the domain), u(x, t) = f (x), for x = (x1 , x2 , . . . , xn ) ∈ ? ?, t > 0.

2. Neumann boundary condition (the derivative of the solution at the direction of outward normal is given) ?u = n · grad(u) = n · ?u = f (x), x =∈ ? ? ?n n = n(x) is the outward unit normal to ? ? at x ∈ ? ?, and grad(u) = ?u = ?u ?u ?u , ,..., . ?x1 ?x2 ?xn

3. Robin boundary condition (a combination of 1 and 2), ?u + k · u(x, t) = f (x), ?n k > 0, x = (x1 , x2 , . . . , xn ) ∈ ? ?.

2 n2 1 + n2 = 1

Example 7. For u = u(x, y ) we have n = (n1 , n2 ), with |n| = and n · ?u = n1 ux + n2 uy . y n = (n1 , n2 ) n2 P ? n1 x

Figure 2: The domain ? and its outward normal n at a point P ∈ ? ?.

14

CHAPTER 0. INTRODUCTION

Example 8. Let u(x, y ) = x2 + y 2. We determine the normal derivative of u in direction v = (1, 1). The gradient of u is the vector valued function ?u = 2x·e1 +2y ·e2 , where e1 = (1, 0) and e2 = (0, 1) are the unit orthonormal basis in R2 : e1 · e1 = e2 · e2 = 1 and e1 · e2 = e2 · e1 = 0. Note that v = e1 + e2 = (1, 1) is not a unit vector. The normalized v is obtained viz vn = v/|v |, i.e. (1, 1) e1 + e2 (1, 1) = √ . vn = =√ |e1 + e2 | 12 + 12 2 Thus with ?u(x, y ) = 2x · e1 + 2y · e2 , we get vn · ?u(x, y ) = which gives (1, 1) 2 x + 2y (1, 1) . vn · ?u(x, y ) = √ · [2x(1, 0) + 2y (0, 1)] = √ · (2x, 2y ) = √ 2 2 2 Thus √ 4 vn · ?u(1, 1) = √ = 2 2. 2 e1 + e2 (2x · e1 + 2y · e2 ). |e1 + e2 |

The usual three questions I.

In theory

1. Existence, at least one solution u 2. Uniqueness, either one solution or no solutions at all 3. Stability, continuous dependence of solutions to the data Remark. A property that concerns behavior, such as growth or decay, of perturbations of a solution as time increases is generally called a stability property.

II.

In applications

1. Construction, of the physical solution. 2. Regularity, how substitutable is the found solution. 3. Approximation, when an exact construction of the solution is impossible.

0.2. TRINITIES Three general approaches to analyzing di?erential equations

15

1. Transformation to a simpler problem: The separation of variables technique to reduce the (PDEs) to simpler eigenvalue problems (ODEs). Also called Fourier method, or solution by eigenfunction expansion (Fourier analysis). 2. The multiplier method: The multiplier method is a strategy for extracting information by multiplying a di?erential equation by a suitable function and then integrating. This usually is referred as variational formulation, or energy method (subject of our study). 3. Green’s Function, fundamental singularities, or solution by integral equations (an advanced PDE course).

16

CHAPTER 0. INTRODUCTION

Chapter 1 Polynomial approximation in 1d

Our objective is to present the ?nite element method as an approximation technique for solution of di?erential equations using piecewise polynomials. This chapter is devoted to some necessary mathematical environments and tools as well as a motivation for the unifying idea of using ?nite elements: A numerical strategy arising from the need of changing a continuous problem into a discrete one. The continuous problem will have in?nitely many unknowns (if one asks for u(x) at every x), and it cannot be solved exactly on a computer. Therefore it has to be approximated by a discrete problem with ?nite number of unknowns. The more unknowns we keep, the better will be the accuracy of the approximation and greater the expences.

1.1

Overture

Below we shall introduce a few standard examples of classical equations and some regularity requirements. Ordinary di?erential equations (ODEs) An initial value problem, (IVP), in population dynamics can be written as u ˙ (t) = λu(t), 0<t<1 u(0) = u0 , (1.1.1)

where u ˙ (t) = du and λ is a positive constant. For u0 > 0 this problem has the dt increasing analytic solution u(t) = u0 eλ·t , which would blow up as t → ∞. 17

18

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D

˙ (t) = F (u(t), t), where u(t) ∈ Rn is a time dependent Generally, we have u n ˙ = ? u(t)/?t ∈ Rn being its componentwise derivative vector in R , with u ˙ (t) = with respect to t ∈ R+ . Thus u(t) = [u1 (t), u2 (t), . . . , un (t)]T , u ′ ′ ′ T [u1 (t), u2 (t), . . . , un (t)] and F : Rn × R+ → Rn . Partial di?erential equations (PDEs) in bounded domains Let ? be a bounded, convex, subset of the Eucledean space Rn . Below is an example of a general boundary value problem in ? ? Rn with the ? Dirichlet boundary condition, ? ? ??u(x) + αb · ?u(x) = f, x ∈ ? ? Rn , (1.1.2) ? u(x) = 0, x ∈ ? ?. where α ∈ R, b = (b1 , b2 , . . . , bn ) ∈ Rn and u : Rn → R is a real-valued ?u ?u ?u ?2u ?2u ?2u function with ?u := ?x , , . . . , , ? u = 2 + ?x2 + . . . + ?x2 , ?x ?x ?x n 1 2 n 1 2 and ?u ?u ?u b · ?u = b1 + b2 + . . . + bn . ?x1 ?x2 ?xn The following Heat equation is an example of a boundary value problem with ? Neumann boundary condition ?u = ?u, ?t x ∈ ? ? Rk , ?u = 0, x ∈ ? ?, ?n

(1.1.3)

where n = (n1 , n2 , . . . , nk ) is the outward unit normal to the boundary ? ? at the point x ∈ ? ?, and ?u = n · ?u. (1.1.4) ?n Regularity requirements for classical solutions 1) u ∈ C 1 : every component of u has a continuous ?rst order derivative. 2) u ∈ C 1 : all ?rst order partial derivatives of u are continuous. 3) u ∈ C 2 : all partial derivatives of u of order 2 are continuous. 4) u ∈ C 1 R+ ; C 2 (?) :

?u ?t ?2u , ?xi ?xj

and

i, j = 1, 2, . . . , n are continuous.

1.1. OVERTURE

19

Remark 2. Note that, we tacitly understand that: u in 1) is a vector-valued function of a single variable as in the example of, general, dynamical system (1.1.1), whereas u in 2)-4) is a scalar (real-valued) function of several variables. ? Numerical solutions of (IVP) Example 9. A ?nite di?erence method. We descritize the IVP (1.1.1) with explicit (forward) Euler method based on a partition of the interval [0, T ] into N subintervals: t0 = 0 t1 t2 t3 tN = T

and an approximation of the derivative by a di?erence quotient at each subin+1 )?u(tk ) terval [tk , tk+1 ] as u ˙ (t) ≈ u(tk . Then an approximation of (1.1.1) is tk+1 ?tk given by u(tk+1) ? u(tk ) = λ · u(tk ), tk+1 ? tk and thus, letting ?tk = tk+1 ? tk , u(tk+1) = 1 + λ?tk u(tk ). (1.1.6) k = 0, 1 , . . . N ? 1 , with u(0) = u0 , (1.1.5)

Starting with k = 0 and the data u(0) = u0 , the solution u(tk ) would, iteratively, be produced at the subsequent points: t1 , t2 , . . . , tN = T . For a uniform partition, where all subintervals have the same length ?t, (1.1.6) would be u(tk+1 ) = 1 + λ?t u(tk ), k = 0, 1 , . . . , N ? 1 . (1.1.7)

There are corresponding ?nite di?erence methods for PDE’s. Our goal, however, is to study the Galerkin’s ?nite element method. To this approach we need to introduce some basic tools: Finite dimensional linear space of functions de?ned on an interval Below we give a list of some examples for ?nite dimensional linear spaces. Some of these examples are studied in details in the interpolation chapter.

20

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D I. P (q) (a, b) := { The space of polynomials in x of degree ≤ q, a ≤ x ≤ b}.

2 3 q A possible basis for P (q) (a, b) would be {xj }q j =0 = {1, x, x , x , . . . , x }. These are, in general, non-orthogonal polynomials and may be orthogonalized by Gram-Schmidt procedure. The dimension of P q is therefore q + 1.

II. An example of orthogonal bases functions, on (0, 1) or (?1, 1) are the Legendre polynomials: dk k 1 dn 2 k Pk (x) = (?1) [x (1 ? x) ] or Pn (x) = n (x ? 1)n , k n dx 2 n! dx

k

respectively. The ?rst four Legendre orthogonal polynomials on (?1, 1) are as follows: 3 1 5 3 P0 (x) = 1, P1 (x) = x, P2 (x) = x2 ? , P3 (x) = x3 ? x. 2 2 2 2 III. Periodic orthogonal bases on [0, T ] are usually represented by trigonometric polynomials given by

N

T

N

:= f (x) f (x) =

n=0

an cos

2π 2π nx + bn sin nx T T

IV. A general form of bases functions on an interval are introduced in the (q ) interpolation chapter: these are Lagrange bases {λi }q (a, b) i=0 ∈ P associated to a set of (q + 1) distinct points ξ0 < ξ1 < . . . < ξq in (a, b) determined by the requirement that ? q ? 1, i = j, x ? ξj . or λi (x) = λ i (ξ j ) = ? 0, ξi ? ξj i=j

j =1,(j =i)

A polynomial p ∈ P (q) (a, b), that has the value pi = p(ξi ) at the nodes x = ξi for i = 0, 1, . . . , q , expressed in terms of the corresponding Lagrange basis is then given by p(x) = p0 λ0 (x) + p1 λ1 (x) + . . . + pq λq (x). (1.1.8)

Note that for each node point x = ξi we have associated a base function λi (x), i = 0, 1, . . . , q . Thus we have q + 1 basis functions.

1.1. OVERTURE

21

Remark 3. Our goal is to approximate general functions, rather than polynomials, by piecewise polynomials of Lagrange type. Then, for a given function f the Lagrange coe?cients in (1.1.8) will be replaced by pi = f (ξi), 1 ≤ i ≤ q , and f (x) will be approximated by its Lagrange interpolant, viz

q

f (x) ≈

f (ξi )λi (x) := πq f (x).

i=0

(1.1.9)

We shall illustrate this in the next examples. Example 10. The linear Lagrange basis functions, q = 1, are given by (see Fig. 1.1.) λ0 (x) = ξ1 ? x ξ1 ? ξ0 and λ1 (x) = x ? ξ0 . ξ1 ? ξ0 (1.1.10)

1 λ0 (x) λ1 (x)

a

ξ0

ξ1

b

x

Figure 1.1: Linear Lagrange basis functions for q = 1. Example 11. A typical application of Lagrange bases is in ?nding a polynomial interpolant πq f ∈ P q (a, b) of a continuous function f (x) on an interval [a, b]. The procedure is as follows: Choose distinct interpolation nodes a = ξ0 < ξ1 < . . . < ξq = b and de?ne πq f (ξi) = f (ξi). Then πq f ∈ P (q) (a, b), de?nied as the sum in (1.1.9), interpolates f (x) at the nodes {ξi}, i = 0, . . . , q and using Lagrange’s formula (1.1.8), with pi = f (ξi), i = 0, 1, . . . , q , and x ∈ [a, b] yields πq f (x) = f (ξ0 )λ0 (x) + f (ξ1)λ1 (x) + . . . + f (ξq )λq (x).

22

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D For linear interpolant, i.e. q = 1, we only need 2 nodes and 2 basis functions. Choosing ξ0 = a and ξ1 = b, in (1.1.10) we get the linear interpolant π1 f (x) = f (a)λ0 (x) + f (b)λ1 (x), where λ0 (x) = i.e., b?x b?a and λ1 (x) = x?a , b?a

π1 f (x) = f (a) y

x?a b?x + f (b) b?a b?a

π1 f (x) f (x) x

a

b

Figure 1.2: The linear interpolant π1 f (x) on a single interval.

V. We shall frequently use the space of continuous piecewise polynomials on a partition of an interval into some subintervals. For example Th : 0 = x0 < x1 < . . . < xM < xM +1 = 1, with hj = xj ? xj ?1 and j = 1, . . . , M + 1, is a partition of [0, 1] into M + 1 subintervals. Let Vh denote the space of all continuous piecewise polynomial func(q ) tions of degree ≤ q on Th . Obviously, Vh ? P (q) (0, 1). Let also V h = { v : v ∈ Vh ,

? (q ) (q ) (q )

v (0) = v (1) = 0}.

Our motivation in introducing these function spaces is due to the fact that these are function spaces, adequate in the numerical study of the

1.1. OVERTURE y

23

x0

x1

h2

x2 xj ?1

hj

xj

xM

hM +1

xM +1 = 1

? (1)

x

Figure 1.3: Fig shows an example of V h .

boundary value problems using ?nite element methods for approximating solutions with piecewise polynomials. The standard basis for piecewise linears: Vh := Vh is given by the so called linear hat-functions ?j (x) with the property that ?j (x) is a piecewise linear function with ?j (xi ) = δij : ? x ? x j ?1 ? ? ? xj ?1 ≤ x ≤ xj ? ? 1, ? hj i = j, xj +1 ?x δij = i.e. ?j (x) = xj ≤ x ≤ xj +1 hj +1 ? ? 0, ? i = j, ? ? 0 x∈ / [xj ?1 , xj +1]. y 1 ?j (x)

(1)

x0

xj ?2 xj ?1

hj

xj

hj +1

xj +1

xM

xM +1

x

Figure 1.4: Fig shows a general piecewise linear basis function ?j (x).

24

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D

Vector spaces To establish a framework we shall use some basic mathematical concepts: De?nition 1. A set of functions or vectors V is called a linear space, or a vector space, if for all u, v ∈ V and all α ∈ R (real numbers), we have (i) u + αv ∈ V (ii) u + v = v + u (iii) ? u ∈ V, ? (?u) ∈ V (1.1.11) such that u + (?u) = 0.

Obseve that (iii) and (i), with α = 1 and v = (?u) implies that 0 (zero vector) is an element of every vector space. De?nition 2 (Scalar product). A scalar product is a real valued operator on V × V , viz (u, v ) : V × V → R such that for all u, v, w ∈ V and all α ∈ R, (i) (ii) u, v = v, u , (symmetry) u + αv, w = u, w + α v, w , (bi-linearity). (1.1.12)

De?nition 3. A vector space W is called a scalar product space if W is associated with a scalar product ·, · , de?ned on W × W . The function spaces C ([0, T ]), C k ([0, T ]), , P q , T q are examples of scalar product spaces associated with usual scalar product de?ned by

T

u, v =

0

u(x)v (x)dx,

(1.1.13)

De?nition 4 (Orthogonality). Two, real-valued, functions u(x) and v (x) are called orthogonal if u, v = 0. This orthogonality is also denoted by u ⊥ v . De?nition 5 (Norm). If u ∈ V then the norm of u, or the length of u, associated with the above scalar product is de?ned by u = u, u = u, u

1 2

T

=

0

|u(x)|2 dx

1 2

.

(1.1.14)

This norm is known as the L2 -norm of u(x). There are other norms that we will introduce later on.

1.1. OVERTURE

25

We also recall one of the most useful tools that we shall frequently use through out this note: The Cauchy-Schwarz inequality, | u, v | ≤ u A simple proof of (1.1.15) is given by using u ? av, u ? av ≥ 0, with a = u, v / v 2. v . (1.1.15)

Then by the de?nition of the L2 -norm and the symmetry property of the scalar product we get 0 ≤ u ? av, u ? av = u 2 u, v = u 2? v 2. v 4 Seting a = u, v / v

2 2

? 2a u, v + a2 v

2

and rearranging the terms we get u, v 2 ≤ u 2, v 2

which yields the desired result. ? Galerkin method for (IVP) For a solution u of the initial value problem (1.1.1) we use test functions v , in a certain vector space V , for which the integrals below are well-de?ned,

T T

u ˙ (t)v (t) dt = λ

0 0

u(t)v (t) dt,

?v ∈ V,

(1.1.16)

or equivalently

T 0

u ˙ (t) ? λ u(t) v (t)dt = 0, u ˙ (t) ? λ u(t) ⊥ v (t),

?v (t) ∈ V, ?v (t) ∈ V.

(1.1.17)

i.e. (1.1.18)

De?nition 6. If w is an approximation of u in the problem (1.1.16), then R w (t) := w ˙ (t) ? λw (t) is called the residual error of w (t)

26

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D

In general for an approximate solution w we have w ˙ (t) ? λw (t) = 0, otherwise w and u would satisfy the same equation and by the uniqueness we would get the exact solution (w = u). Our requirement is instead that w satis?es the equation (1.1.1) in average, or in other words we require that w satis?es (1.1.18), i.e, R w (t) ⊥ v (t), ?v (t) ∈ V. (1.1.19)

In our case the real solution belongs to C ((0, T )), or better to

s T

H s (0, T ) := {f :

? k f /?tk

k =0 0

2

dt < ∞}.

H s is a subspace of C ((0, T )) consisting of all function in L2 (0, T ) having also all their derivatives of order ≤ s in L2 (0, T ). We look for a solution U (t) in a ?nite dimensional subspace e.g. V (q) . More speci?cally, we want to look at an approximate solution U (t), called a trial function for (1.1.1) in the space of piecewise polynomials of degree ≤ q : V (q) = {U : U (t) = ξ0 + ξ1 t + ξ2 t2 + . . . + ξq tq }. (1.1.20)

Hence, to determine U (x) we need to determine the coe?cients ξ0 , ξ1 , . . . ξq . We refer to V (q) as the space of trial functions. Note that u(0) = u0 is given and therefore we may take U (0) = ξ0 = u0 . It remains to ?nd the real numbers ξ1 , . . . , ξq . These are coe?cients of q linearly independent monomials t, t2 , . . . , tq . To this approach we de?ne the test functions space: V

? (q )

= {v ∈ V (q) : v (0) = 0},

(1.1.21)

in other words v can be written as v (t) = ξ1 t + ξ2 t2 + . . . + ξq tq . Note that V

? (q )

= span[t, t2 , . . . , tq ].

(1.1.22)

For an approximate solution U we require that its residual R(U ) satisfy the orthogonality condition (1.1.19): R U (t) ⊥ v (t), ?v (t) ∈ V

? (q )

.

1.1. OVERTURE

27

Thus the Galerkin method for (1.1.1) is formulated as follows: Given u(0) = u0 , ?nd the approximate solution U (t) ∈ V (q) , for (1.1.1) such that (for simplicity we put T ≡ 1)

1 0 1

R U (t) v (t)dt =

0

˙ (t) ? λ U (t))v (t)dt = 0, ?v (t) ∈ V (U

? (q )

. (1.1.23)

Formally this can be obtained writing a wrong!!! equation by replacing u by U ∈ V (q) in (1.1.1), ? ? U ˙ (t) = λ U (t), 0<t<1 (1.1.24) ? U (0) = u ,

0

then, multiplying (1.1.24) by a function v (t) ∈ V space and integrating over [0, 1].

q j =1 ? (q )

? (q )

from the test function

q j j =1 ξj t ,

Now since U ∈ V (q) we can write U (t) = u0 +

˙ (t) = then U

jξj tk?1 . Further we have V is spanned by vi (t) = ti , i = 1, 2, . . . , q . ˙ and v = vj , j = 1, 2, . . . , q in (1.1.22) Inserting these representations for U, U we get

1 0 q q

kξj t

j =1

j ?1

? λu0 ? λ

j =1

ξj tj · ti dt = 0,

i = 1, 2, . . . , q,

(1.1.25)

which can be rewritten as

1 0 q 1

(jξj t

j =1

i+j ?1

? λ ξj t

i+j

) dt = λu0

0

ti dt.

(1.1.26)

Performing the integration (ξj :s are constants independent of t) we get ti+j +1 ti+j ?λ ξj j · i+j i+j+1 j =1 or equivalently

q q t=1 t=0

ti+1 = λ · u0 i+1

t=1 t=0

,

(1.1.27)

j =1

j λ λ ξj = ? · u0 i+j i+j+1 i+1

i = 1, 2, . . . , q,

(1.1.28)

28

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D

which is a linear system of equations with q equations and q unknowns (ξ1 , ξ2 , . . . , ξq ); given in the matrix form as AΞ = b, with A = (aij ), Ξ = (ξj )q j =1 , and b = (bi )q i=1 . (1.1.29)

But the matrix A although invertible, is ill-conditioned, mostly because {ti }q i=1 does not form an orthogonal basis. We observe that for large i and j λ j ? , are very the two last rows (columns) of A given by aij = i+j i+j+1 close to each others resulting to extreme small values for the determinant of A. If we insist to use polynomial basis up to certain order, then instead of monomials, the use of Legendre orthogonal polynomials would yield a diagonal (sparse) coe?cient matrix and make the problem well conditioned. This however, is a rather tedious task. Galerkin’s method and orthogonal projection Let u = (u1 , u2 , u3 ) ∈ R3 and assume that for some reasons we only have u1 and u2 available. Letting x = (x1 , x2 , x3 ) ∈ R3 , the objective, then is to ?nd U ∈ {x : x3 = 0}, such that (u ? U ) is as small as possible. For orthogonal projection we have z · n = 0, for all z ∈ {x : x · n = 0, x3 = 0}, where n is the normal vector to the plane {(x1 , x2 , 0)}. Obviously in this case U = (u1 , u2 , 0) and we have (u ? U ) ⊥ U . Note that, if m < n, and um is the projection of u = (u1 , u2, . . . , un?1, un ) on Rm , then um = (u1 , u2 , . . . , um , um+1 = 0, . . . , un = 0), and the Euclidean 2 2 distance: |u ? um | = u2 m+1 + um+2 + . . . + un → 0 as m → n. This meams the obvious fact that the accuracy of the orthogonal projection will improve by raising the dimension of the projection space. ? Galerkin method for (BVP) We consider the Galerkin method for the following stationary (u ˙ = du/dt = 0) heat equation in one dimension: ? ? ? d a(x) · d u(x) = f (x), 0 < x < 1; dx dx (1.1.30) ? u(0) = u(1) = 0. Let a(x) = 1, then we have ?u′′ (x) = f (x), 0 < x < 1; u(0) = u(1) = 0. (1.1.31)

1.1. OVERTURE x3 u = (u1 , u2, u3 ) n=u?U x2 U = (u1 , u2 , 0)

29

x1 Figure 1.5: Example of a projection on R2 .

Let now Th : 0 = x0 < x1 < . . . < xM < xM +1 = 1 be a partition of the interval (0, 1) into the subintervals Ij = (xj ?1 , xj ), with lengths |Ij | = hj = xj ? xj ?1 , j = 1, 2, . . . , M . We de?ne the ?nite dimensional space Vh0 = {v ∈ C (0, 1) : v is piecewise linear function on Th , and v (0) = v (1) = 0}, with the bases functions {?j }M j =1 . Due to the fact that u is known at the boundary points 0 and 1; it is not necessary to supply basis functions corresponding to the values at x0 = 0 and xM +1 = 1. Remark 4. If the Dirichlet boundary condition is given at only one of the boundary points; say x0 = 0 and the other one satis?es, e.g. a Neumann condition as ?u′′ (x) = f (x), 0 < x < 1; u(0) = b0 , u′ (1) = b1 , (1.1.32)

then the corresponding basis function ?0 will be unnecessary (no matter whether b0 = 0 or b0 = 0), whereas one needs to provide the half-base function ?M at xM +1 = 1 (dashed in the Fig below).

Now the Galerkin method for problem (1.1.31) (which is just the description of the orthogonality condition of the residual R(U ) = ?U ′′ ? f to the

30 y 1

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D

?1

?j

?M

?M +1

x0

x1

x2

xj ?1

xj hj

xj +1 hj +1

xM ?1

xM xM +1

x

Figure 1.6: General piecewise linear basis functions

test function space Vh0 ; i.e., R(U ) ⊥ Vh0 ) is formulated as follows: Find the approximate solution U (x) ∈ Vh0 such that

1 0

(?U ′′ (x) ? f (x))v (x)dx = 0,

?

?v (x) ∈ Vh0

(1.1.33)

Observe that if U (x) ∈ V h , then U ′′ (x) is either equal to zero or is not a well-de?ned function, in the latter case, the equation (1.1.33) does not make sense, whereas for U ′′ (x) = 0 and the data U (0) = U (1) = 0 we get the trivial approximation U (x) ≡ 0. This however, is relevant only if f ≡ 0, but then even u(x) ≡ 0 and we have a trivial case. If, however, we perform partial integration then

1 1

?

U (x)v (x)dx =

0 ? 0

′′

U ′ (x)v ′ (x)dx ? [U ′ (x)v (x)]1 0

(1.1.34)

and since v (x) ∈ V h ; v (0) = v (1) = 0, we get

1 1

?

U (x)v (x)dx =

0 ? 0

′′

U ′ (x)v ′ (x) dx

(1.1.35)

Now for U (x), v (x) ∈ V h , U ′ (x) and v ′ (x) are well-de?ned (except at the nodes) and the equation (1.1.33) has a meaning.

1.1. OVERTURE

31

Hence, The Galerkin ?nite element method (FEM) for the problem (1.1.30) is now reduced to: Find U (x) ∈ Vh0 such that

1 1

U ′ (x)v ′ (x) dx =

0 0

f (x)v (x)dx,

?v (x) ∈ Vh0 .

(1.1.36)

We shall determine ξj = U (xj ) the approximate values of u(x) at the node points xj . To this end using basis functions ?j (x), we may write

M M

U (x) =

j =1

ξj ?j (x)

which implies that

U (x) =

j =1

′

ξj ?′ (x). (1.1.37)

Thus we can write (1.1.36) as

M 1 1

ξj

j =1 0

?′j (x) v ′ (x)dx

=

0

f (x)v (x)dx,

?v (x) ∈ Vh0 .

(1.1.38)

Since every v (x) ∈ Vh0 is a linear combination of the basis functions ?j (x), it su?ces to try with v (x) = ?i (x), for i = 1, 2, . . . , M : That is to ?nd ξj (constants), 1 ≤ j ≤ M such that

M 0 1 1

?′i (x)

j =1

·

?′j (x)dx

ξj =

0

f (x)?i (x)dx,

i = 1, 2, . . . , M. (1.1.39)

This equation can be written in the equivalent matrix form as Aξ = b. Here A is called the sti?ness matrix and b the load vector:

1

(1.1.40)

A= ? b1 ?

{aij }M i,j =1 ,

aij =

0

?′j (x) · ?′i (x)dx, ? ξ1 ?

(1.1.41)

? ? ? ? ? b2 ? ? , with ? b=? ? ? ... ? ? ? bM

1

bj =

0

? ? ? ? ? ξ2 ? ? . (1.1.42) ? f (x)?i (x)dx, and ξ = ? ? ? ... ? ? ? ξM

32

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D

To compute the entries aij of the sti?ness matrix A, ?rst we need to determine ?′j (x). Note that ? ? ? ? ? ? ? ? ?

x?xi?1 hi xi+1 ?x hi+1

?i (x) =

xi?1 ≤ x ≤ xi xi ≤ x ≤ xi+1 else

=?

?′i (x) =

0

? ? ? ? ? ? ? ? ?

1 hi

? hi1 xi ≤ x ≤ xi+1 +1 0 else

xi?1 ≤ x ≤ xi

Sti?ness matrix A: If |i ? j | > 1, then ?i and ?j have disjoint non-overlapping supports, evideltly, we hence

1

aij =

0

?′i (x) · ?′j (x)dx = 0.

y 1 ?j ?1 ?j +1

xj ?2

xj ?1

xj

xj +1

xj +2

x

Figure 1.7: ?j ?1 and ?j +1 . As for i = j : we have that

hi xi hi+1

aii =

xi?1

1 hi

2

xi+1

dx+

xi

?

1 hi+1

2

dx =

1 1 xi ? xi?1 xi+1 ? xi + = + . 2 2 hi hi+1 hi hi+1

It remains to compute aij for the case j = i +1: A straightforward calculation (see the ?g below) yields

xi+1

ai,i+1 =

xi

?

1 hi+1

·

1 hi+1

dx = ?

xi+1 ? xi 1 =? . 2 hi+1 hi+1

(1.1.43)

1.1. OVERTURE y 1 ?j ?j +1

33

xj ?1

xj

xj +1

xj +2

x

Figure 1.8: ?j and ?j +1.

. Obviousely ai+1,i = ai,i+1 = ? hi1 +1 To summarize, we have

? ? ? a = 0, ? ? ij

1 aii = h + hi1 , i +1 ? ? ? 1 ? a i?1,i = ai,i?1 = ? hi ,

if |i ? j | > 1, i = 1, 2, . . . , M, i = 2, 3, . . . , M.

(1.1.44)

Thus by symmetry we ?nally have the sti?ness matrix for the stationary heat conduction as: ? ? ? ? ? ? ? ? ? ? ? ? (1.1.45)

? ? ? ? ? A=? ? ? ? ?

1 h1

+

1 h2

? h12 0 ... 0

1 h2

1 ?h 2

0 ? h13 ... ... 0

... 0 ... ... ? h1 M

1 hM

0 0 0 ? h1 M +

1 hM +1

+ ... 0 ...

1 h3

34

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D we get that ?1 2 ?1 ... 0 ?1 2 ... ... 0 ?1 ?1 0 0 ? ? ? ? ? ? ? ? ? ? ? ? ? ? (1.1.46)

With a uniform mesh, i.e. hi = h ? 2 ? ? ? ?1 ? ? 0 1 ? Aunif = · ? h ? ? ... ? ? ? ... ? 0

1 xi

... ... 0 ... 0 ?1 2

... ... ... ... 0 2 ?1

... ...

As for the components of the load vector b we have bi =

0

f (x)?i (x) dx =

x i ?1

f (x)

x ? xi?1 dx + hi

xi+1

f (x)

xi

xi+1 ? x dx. hi+1

? A ?nite di?erence approach To illustrate a ?nite di?erence approach we reconsider the stationary heat equation (1.1.31): ?u′′ (x) = f (x), 0 ≤ x ≤ 1; (1.1.47)

and motivate for its boundary conditions. The equation (1.1.47) is linear for the unknown function u, with inhomogeneous term f . There is some arbitrariness left in the problem, because any combination C + Dx could be added to any solution. The sum would constitute another solution, since the second derivative of C + Dx contributes nothing. Therefore the uncertainity left by these two arbitrary constants C and D will be removed by adding a boundary condition at each end point of the interval u(0) = 0, u(1) = 0. (1.1.48)

The result is a two-point boundary-value problem, describing not a transient but a steady-state phenomenon–the temperature distribution in a rode, for example with ends ?xed at 0 and with a heat source f (x). As our goal is to solve a discrete problem, we cannot accept more than a ?nite amount of information about f , say it values at equally spaced points x = h, x = 2h, . . . , x = nh. And what we compute will be approximate

?

1.1. OVERTURE

35

values u1 , u2 , . . . , un for the true solution u at these same points. At the ends x = 0 and x = 1 = (n + 1)h, we are already given the correct boundary values u0 = 0, un+1 = 0. The ?rst question is, How do we replace the derivative d2 u/dx2 ? Since every derivative is a limit of di?erence quotients, it can be approximated by a stopping at a ?nite stepsize, and not permitting h (or ?x) to approach zero. For du/dx there are several alternatives: du u(x + h) ? u(x) ≈ dx h or u(x) ? u(x ? h) h or u(x + h) ? u(x ? h) . 2h

The last, because it is symmetric about x, is the most accurate. For the second derivative there is just one combination that uses the values at x and x ± h: u(x + h) ? 2u(x) + u(x ? h) d2 u ≈ . (1.1.49) 2 dx h2 It also has the merit of being symmetric about x. (1.1.49) is obtained using d2 u u′ (x) ? u′(x ? h) ≈ . d2 x h Replacing the approximations u′ (x) ≈ in (1.1.49) we get

u(x+h)?u(x) h

(1.1.50)

u(x)?u(x?h) h

and u′ (x ? h) ≈

(u(x + h) ? u(x))/h ? (u(x) ? u(x ? h))/h d2 u ≈ 2 dx h u(x + h) ? 2u(x) + u(x ? h) = . h2

(1.1.51)

To repeat the right side approaches the true value of d2 u/dx2 as h → 0, but have to stop at a positive h. At a typical meshpoint x = jh, the di?erential equation ?d2 u/dx2 = f (x) is now replaced by this discrete analogue (1.1.51); after multiplying by h2 , ?uj +1 + 2uj ? uj ?1 = h2 f (jh). (1.1.52)

There are n equations of exactly this form, for every value j = 1, 2, . . . , n. The ?rst and last equations include the expressions u0 and un+1 , which are not unknowns–Their values are the boundary conditions, and they are shifted to the right side of the equation and contribute to the inhomogeneous terms

36

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D

(or at least, they might, if they were not known to be equal zero). It is easy to understand (1.1.52) as a steady-state equation, in which the ?ows (uj ? uj +1) coming from the right and (uj ? uj ?1 ) coming from the left are balanced by the loss h2 f (jh) at the center. The structure of the n equations (1.1.52) can be better visualized in matrix form Au = b viz ? ? ? ?? ? f (h) u1 2 ?1 0 . . . . . . 0 ? ? ? ?? ? ? ? ? ?? ? ? f (2h) ? ? ?1 2 ?1 0 . . . . . . ? ? u2 ? ? ? ? ?? ? ? ? ? ?? ? ? f (3h) ? ? 0 ?1 2 ?1 0 . . . ? ? u3 ? 2 ?, ?=h ? ?? ? (1.1.53) ? ? ? ?? ? ? ? ? ... ... ... ... ... 0 ?? · ? · ? ? ? ?? ? ? ? ? ?? ? ? ? ? . . . . . . 0 ?1 2 ?1 ? ? · ? · ? ? ? ?? ? f (nh) un 0 . . . . . . 0 ?1 2 which, once again, gives the structure of our uniform stifness matrix Aunif given in (1.1.46). So we conclude that, for this problem, the ?nite element and ?nite di?erence approximations are two equivalent approaches.

Remark 5. Unlike the matrix A for monomial approximation of IVP in (1.1.28), A has more desirable structure, e.g. A is a sparse, tridiagonal and symmetric matrix. This is due to the fact that the basis functions {?j }M j =1 are nearly orthogonal. The most important property of A is that it is positive de?nite. De?nition 7. A matrix A is called positive de?nite if

M

?η ∈ R , η = 0, η Aη > 0,

M

T

i.e.

i,j =1

ηi aij ηj > 0.

(1.1.54)

We shall use the positive de?niteness of A to argue that (1.1.40) is uniquely solvable. To this approach we prove the following well-known result: Proposition 1. If a square matrix A is positive de?nite then A is invertible and hence Aξ = b has a unique solution.

1.1. OVERTURE

37

Proof. Suppose Ax = 0 then xT Ax = 0, and since A is positive de?nite, then x ≡ 0. Thus A has full range and we conclude that A is invertible. Since A is invertible Aξ = b has a unique solution: ξ = A?1 b. Note however, that it is a bad idea to invert a matrix to solve the linear equation system. Finally we illustrate an example of the positive-de?niteness argument for Aunif . ? ? x Example 12. Assume M = 2 and let U (x, y ) = ? ? , then y U T Aunif U = (x, y ) ?

2

?

2 ?1

?1 2

?? ??

2

x y

?

= 2x ? xy ? xy + 2y = x + y + x ? 2xy + y = x2 + y 2 + (x ? y )2 ≥ 0.

2

? = (x, y ) ?

2 2

?

2x ? y ?x + 2y

2

? ?

(1.1.55)

Thus Aunif is positive de?nite. Since U T AU = 0 only if x = y = 0 i.e. U = 0. Remark 6. For a boundary value problem with, e.g. inhomogeneous Dirichlet boundary data, actually a direct approach would be with the test and trial function spaces having di?erent dimensions; test functions are zero at the boundary: v ∈ V h , trial function: u ∈ Vh (not necessarily zero at the boundary). This would yield a linear equation system AΞ = b with a rectangular matrix A instead of a quadratic one. Then to continue we use the least square method and instead solve AT AΞ = AT b. The solution Ξ is approximate and AΞ = b. Thus, the corresponding orthogonality condition of the residual to the test function space is now r := (AΞ ? b) ⊥ Cj , where Cj are columns in A. Summary: Roughly speaking, a systematic procedure for approximate solution for a di?erential equations would involve the following steps: 1. We need to approximate functions by polynomials agreeing with the functional values at certain points (nodes). This is the matter of Interpolation techniques which we shall introduce in the next chapter.

?

38

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D 2. The function f (x) is unknown and the elements of the vector b as well as the integrals that represent the elements of the coe?cient matrix are of involve character, for example when approximating by higher order polynomials and/or solving equations with variable coe?cients. Therefore we need to approximate di?erent integrals over subintervals of a partition. This may be done using Gauss quadrature rules. In simple case one may use usual or composite midpoint-, trapezoidal-, or Simpson’s-rules. In more involved cases one may employ Gauss quadrature rules. We shall brie?y introduce the idea of Gauss quadrature rule in the next chapter. 3. Finally we end up with linear systems of equations (LSE) of type (1.1.40). To solve LSE e?ciently we may use exact Gauss - elimination or the iteration procedures as Gauss-Seidel, Gauss-Jacobi or Over-relaxation methods. We discuss these concepts in the chapter of the numerical linear algebra.

1.2. EXERCISES

39

1.2

Exercises

Problem 1. Use the method of least squares to solve the following systems of linear equations. ? ? ? ? ? ? x1 + x2 = 3 ?x1 + x2 = 16 ? ? ? ? b. a. ?2x1 + 3x2 = 1 2x1 + x2 = ?9 ? ? ? ? ? ? ? 2x ? x ? x2 = 2 ? 2x2 = ?12 1 1 ? ? ? x1 ? ? ? ? ? ? ?2x1 ? ? c. x1 ? ? ? ? ? ?x1 ? ? ? ? ? 2x 1 + 2x2 = + x2 3 = ?4 = ?1 = 5 ? ? ? x1 + x2 + x3 ? ? ? ? ? ?x + x + x3 1 2 d. ? ? ?x2 + x3 ? ? ? ? ? x + x3 1 7 = 4 = 0 = 1 = 2

? 3x2 = ?2 + + x2 x2

Problem 2. Determine the line y = b + ct that ?ts the following pairs of data (t, y ) best. t y 1 1 2 3 5 2 4 5

? ? ? x1 ? ? ? ? ? x 1 e. ? ? x1 ? ? ? ? ? x1

+ x2 + x3 =

+ x2 ? x3 = ?1 ? x2 + x3 = ? x2 ? x3 = 1 3

a.

7 10

b.

t y

1 5

2

3

4

5 17

6 10 12

40 t y 1 2 2 3

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D 3 4 1 1 5 -2

c.

Problem 3. Determine the parameters a and b such that the parabolic curve y = ax2 + bx + c ?ts the following values of x and y best in the least squares sense. a. x -2 -1 0 y 2 1 1 1 2 2 3 b. x -1 0 1 y 2 2 1 2 0

Problem 5. Set up and solve the linear least squares system Ax ≈ b for ?tting the model function f (t, x) = x1 t + x2 et to the three data points (1, 2) (2, 3) and (3, 5). Problem 6. True or false: At the solution to a linear least squares problem Ax ≈ b, the residual vector r = b ? Ax is orthogonal to the column space of A. Problem 7. We want to ?nd a solution approximation U (x) to ?u′′ (x) = 1, 0 < x < 1, u(0) = u(1) = 0, using the ansatz U (x) = A sin πx + B sin 2πx.

Let r ? b ? Ax be the corresponding residual. Which of the following three vectors is a possible value for r ? ? ? ? ? ? ? ?1 ?1 1 ? ? ? ? ? ? ? ? ? ? ? ? ? 1? ??1? ?1? ? ? ? ? ? c. b. a. ? ? ? ? ? ? ? ? 1? ? 1? ?1? ? ? ? ? ? ? ?1 1 1

Problem 4. Let x be the solution of the ? 1 ? ? ?1 A=? ? ?1 ? 1

least squares problem Ax ≈ b, where ? 0 ? ? 1? ?. ? 2? ? 3

1.2. EXERCISES a. Calculate the exact solution u(x). b. Write down the residual R(x) = ?U ′′ (x) ? 1 c. Use the orthogonality condition

1

41

R(x) sin πnx dx = 0,

0

n = 1, 2 ,

to determine the constants A and B . d. Plot the error e(x) = u(x) ? U (x). Problem 8. Consider the boundary value problem ?u′′ (x) + u(x) = x, 0 < x < 1, u(0) = u(1) = 0.

a. Verify that the exact solution of the problem is given by u(x) = x ? sinh x . sinh 1

b. Let U (x) be a solution approximation de?ned by U (x) = A sin πx + B sin 2πx + C sin 3πx, where A, B , and C are unknown constants. Compute the residual function R(x) = ?U ′′ (x) + U (x) ? x. c. Use the orthogonality condition

1

R(x) sin πnx dx = 0,

0

n = 1, 2 , 3 ,

to determine the constants A, B , and C . Problem 9. Let U (x) = ξ0 φ0 (x) + ξ1 φ1 (x) be a solution approximation to ?u′′ (x) = x ? 1, 0 < x < π, u′ (0) = u(π ) = 0,

where ξi , i = 0, 1, are unknown coe?cients and x φ0 (x) = cos , 2 φ1 (x) = cos 3x . 2

42

CHAPTER 1. POLYNOMIAL APPROXIMATION IN 1D

a. Find the analytical solution u(x). b. De?ne the approximate solution residual R(x). c. Compute the constants ξi using the orthogonality condition

1

R(x) φi (x) dx = 0,

0

i = 0, 1 ,

i.e., by projecting R(x) onto the vector space spanned by φ0 (x) and φ1 (x). Problem 10. Use the projection technique of the previous exercises to solve ?u′′ (x) = 0, 0 < x < π, u(0) = 0, u (π ) = 2 ,

2 2 x. π2

assuming that U (x) = A sin x + B sin 2x + C sin 3x +

Chapter 2 Polynomial Interpolation in 1d

2.1 Preliminaries

We recall the idea of polynomial interpolation. Consider a real-valued function f , de?ned on an interval I = [a, b], and a partition Th : a = x0 < x1 < . . . < xM +1 = b, of I into M + 1 subintervals Ij = [xj ?1 , xj ], j = 1, . . . , M + 1. De?nition 8. An interpolant πq f of f on the partition Th is a piecewise polynomial function of degree ≤ q , having the nodal values at xj , j = 1, . . . , M +1, coinciding with those of f : πq f (xj ) = f (xj ) . Here are some simple examples: Linear interpolation on an interval. We start with the unit interval I = [0, 1], without any partitions, and a function f : [0, 1] → R, which is Lipschitz continuous. We let q = 1 and seek the linear interpolation of f on I , i.e. π1 f ∈ P 1 , such that π1 f (0) = f (0) and π1 f (1) = f (1). Thus we seek the constants C0 and C1 in the following representation of π1 f ∈ P 1 , π1 f (x) = C0 + C1 x, where π1 f (0) = f (0) π1 f (1) = f (1) =? =? C0 = f (0), C0 + C1 = f (1) =? C1 = f (1) ? f (0). 43 (2.1.2) x ∈ I, (2.1.1)

44

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

Inserting C0 and C1 in (2.1.1) it follows that π1 f (x) = f (0)+ f (1)?f (0) x = f (0)(1?x)+f (1)x := f (0)λ0 (x)+f (1)λ1(x). In other words π1 f (x) is represented in two di?erent basis: π1 f (x) = C0 +C1 x = C0 ·1+C1 ·x, and π1 f (x) = f (0)(1?x)+f (1)x, with {1?x, x} as the set of basis functions. with {1, x} as the set of basis functions

Note that the functions λ0 (x) = 1 ? x and λ1 (x) = x are linearly independent. Since we can easily see that, if 0 = α0 (1 ? x) + α1 x = α0 + (α1 ? α0 )x, then x=0 x=1 =? =? α0 = 0 α1 = 0 =? α0 = α1 = 0. (2.1.4) x ∈ I, (2.1.3)

f (x) π1 f (x)

λ0 (x) = 1 ? x λ1 (x) = x

1

1

Figure 2.1: Linear interpolation and basis functions for q = 1.

Remark 7. Note that if we de?ne a scalar product on P k (a, b) by

b

(p, q ) =

a

p(x)q (x) dx,

?p, q ∈ P k (a, b),

(2.1.5)

then neither {1, x} nor {1 ? x, x} are orthogonal in the interval [0, 1]. For 2 1 1 example, (1, x) := 0 1 · x dx = [ x ]= 2 = 0. 2

2.1. PRELIMINARIES Now it is natural to ask the following question.

45

Ouestion 1. What will be the error, if one approximates f (x) by π1 f (x)? In other words: what is f (x) ? π1 f (x) =? To answer this question, we need to have a measuring instrument to quantify the di?erence. Gra?cally (geometrically), the deviation of f (x) from π1 f (x) (from at being linear) depends on the curvature of f (x), i.e. on how curved f (x) is. In other words how large is f ′′ (x) on say (a, b)? To this end below we introduce some measuring environments for vectors and functions: De?nition 9 (Vector norm). Let x = (x1 , . . . , xn ), y = (y1 , . . . , yn ) ∈ Rn be two column vectors. We de?ne the scalar product of x and y by x, y = xT y = x1 y1 + · · · + xn yn , and the vector norm for x: x as x = x, x =

2 x2 1 + · · · + xn .

Lp (a, b)-norm: Assume that f is a real valued function such that the integrals as well as the max on the right hand sides below are well-de?ned. Then we de?ne the Lp -norm (1 ≤ p ≤ ∞) of f as

b

Lp -norm L∞ -norm

f f

Lp (a,b)

=

a

|f (x)|p dx

1/p

,

1 ≤ p < ∞,

L∞ (a,b)

= max |f (x)|.

x∈[a,b]

Now, we want to see how far can one answer the question 1 in the Lp -norm environment? Theorem 1. (L∞ -error estimates for the linear interpolation in an interval) Assume that f ′′ ∈ L∞ [a, b]. Then, for q = 1, i.e. only 2 interpolation nodes (the end-points of the interval), there are interpolation constants, ci , independent of the function f (x) and the interval [a, b] such that (1) π1 f ? f (2) π1 f ? f

L∞ (a,b) L∞ (a,b)

≤ ci (b ? a)2 f ′′ ≤ ci (b ? a) f ′

L∞ (a,b)

L∞ (a,b) L∞ (a,b) .

(3) (π1 f )′ ? f ′

L∞ (a,b)

≤ ci (b ? a) f ′′

46

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

Proof. Note that every linear function on [a, b] can be written as a linear combination of the basis functions λa (x) and λb (x) where x?b x?a and λb (x) = . (2.1.6) a?b b?a We point out that linear combinations of λa (x) and λb (x) gives the basis functions {1, x}: λa (x) = λa (x) + λb (x) = 1, aλa (x) + bλb (x) = x. (2.1.7) Now, π1 f (x) is a linear function connecting the two points (a, f (a)) and (b, f (b)), π1 f (x) f (x) λa (x) + λb (x) = 1

1

λa (x) =

b?x b?a

λb (x) =

x ?a b?a

a

b

a

x b

Figure 2.2: Linear Lagrange basis functions for q = 1. which can be represented by π1 f (x) = f (a)λa (x) + f (b)λb (x). By the Taylor’s expansion for f (a) and f (b) about x we can write 1 f (a) = f (x) + (a ? x)f ′ (x) + (a ? x)2 f ′′ (ηa ), ηa ∈ [a, x] 2 1 f (b) = f (x) + (b ? x)f ′ (x) + (b ? x)2 f ′′ (ηb ), ηb ∈ [x, b]. 2 Inserting f (a) and f (b) from (2.1.9) in (2.1.8), it follows that 1 π1 f (x) =[f (x) + (a ? x)f ′ (x) + (a ? x)2 f ′′ (ηa )]λa (x)+ 2 1 +[f (x) + (b ? x)f ′ (x) + (b ? x)2 f ′′ (ηb )]λb (x). 2 (2.1.9) (2.1.8)

2.1. PRELIMINARIES

47

Rearranging the terms and using (2.1.7) (the identity (a ? x)λa (x) + (b ? x)λb (x) = 0) we get π1 f (x) = f (x)[λa (x) + λb (x)] + f ′ (x)[(a ? x)λa (x) + (b ? x)λb (x)]+ 1 1 + (a ? x)2 f ′′ (ηa )λa (x) + (b ? x)2 f ′′ (ηb )λb (x) = 2 2 1 1 = f (x) + (a ? x)2 f ′′ (ηa )λa (x) + (b ? x)2 f ′′ (ηb )λb (x), 2 2 and consequently |π1 f (x) ? f (x)| = 1 1 (a ? x)2 f ′′ (ηa )λa (x) + (b ? x)2 f ′′ (ηb )λb (x) . (2.1.10) 2 2

Now, for a ≤ x ≤ b we have (a ? x)2 ≤ (a ? b)2 and (b ? x)2 ≤ (a ? b)2 , furthermore λa (x) ≤ 1, λb (x) ≤ 1. Finally, by the de?nition of the maximum norm f ′′ (ηa ) ≤ f ′′ L∞ (a,b) , f ′′ (ηb ) ≤ f ′′ L∞ (a,b) . Thus (2.1.10) can be estimated as 1 |π1 f (x) ? f (x)| ≤ (a ? b)2 · 1 · f ′′ 2 and hence |π1 f (x) ? f (x)| ≤ (a ? b)2 f ′′

L∞ (a,b) L∞ (a,b) +

1 (a ? b)2 · 1 · f ′′ 2

L∞ (a,b) ,

(2.1.11)

with

ci = 1 .

(2.1.12)

The other two estimates are proved similarly. Theorem 1 can be proved in a more general setting, for an arbitrary subinterval of I = (a, b), in Lp -norm, 1 ≤ p ≤ ∞. This, general version ( concisly stated below as Theorem 2), is the mainly used Lp -interpolation error estimate. The proof is however, just a scaling of the argument used in the proof of Theorem 1. Remark 8. For a uniform partition Th : a = x0 < x1 < x2 < . . . < xn < xn+1 = b of the interval [a, b] with mesh parameter h = |xj +1 ? xj |, j = 0, 1, . . . , n, it is customary to denote the interpolation function by πh v (x) rather than πq v (x). Here the subindex h refers to the mesh size h, and not to the degree of interpolating polynomial q . The degree q and the meaning of the notation used for the interpolant will be clear from the context.

48

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

Theorem 2. Let πh v (x) be the piecewise linear interpolant of the function v (x) on the partition Th . That is πh v (xj ) = v (xj ), for j = 0, 1, . . . , N + 1. Then, assuming that v is su?ciently regular (v ∈ C 2 (a, b)), there are interpolation constants ci such that πh v ? v Lp ≤ ci h2 v ′′ Lp , (πh v )′ ? v ′ Lp ≤ ci hv ′′ Lp , πh v ? v Lp ≤ ci hv ′ Lp . p = 1, 2 , . . . , ∞ , (2.1.13) (2.1.14) (2.1.15)

For p = ∞ this is just the previous theorem, applied to each subinterval. Note that for a uniform mesh we have h constant and therefore in this case h can be written outside the norms on the right hand sides.. Proof. This is a straightforward generalization of the proof of the Theorem 1 and left as an excercise. Below we review a simple piecewise linear interpolation procedure on a partition of an interval: Vector space of piecewise linear functions on an interval. Given I = [a, b], let Th : a = x0 < x1 < x2 < . . . < xn < xn+1 = b be a partition of I into subintervals Ij = (xj ?1 , xj ) of length hj = |Ij | := xj ? xj ?1 ; j = 1, 2, . . . , N . Let Vh := {v |v is continuous piecewise linear function onTh }, (2.1.16) then Vh is a vector space with the, previously introduced hat functions: { ?j } N j =0 as basis functions. Note that ?0 (x) and ?N (x) are left and right half-hat functions, respectively. It is easy to show that every function in Vh is a linear combination of ?j s: Lemma 1. We have that

N

?v ∈ Vh ;

v (x) =

j =0

v (xj )?j (x),

=? dimVh = N + 1 .

(2.1.17)

Proof. Both left and right hand side are continuous piecewise linear functions. Thus it su?ces to show that they have the same nodal values: Let x = xj then RHS =v (x0 )?0 (xj ) + v (x1 )?1 (xj ) + . . . + v (xj ?1 )?j ?1(xj ) + v (xj )?j (xj ) + v (xj +1 )?j +1(xj ) + . . . + v (xN )?N (xj ) (2.1.18) =v (xj ) = LHS,

2.1. PRELIMINARIES where we have used the fact that ? is piecewise linear and ?i (xj ) = δij .

49

De?nition 10. Assume that f is a Lipschitz continuous function in [a, b]. Then the continuous piecewise linear interpolant of f is de?ned by πh f (xj ) = f (xj ), Or, alternatively:

N

j = 0, 1, . . . , N.

(2.1.19)

πh f (x) =

j =0

f (xj )?j (x),

x ∈ [a, b].

Note that for each interval Ij , j = 1, . . . , N , we have that (i) πh f (x) is linear on Ij , =? πh f (x) = c0 + c1 x on Ij .

(ii) πh f (xj ?1 ) = f (xj ?1 ) and πh f (xj ) = f (xj ). f (x)

πh f (x)

x0

x1

x2

xj

xN ?1

xN

x

Figure 2.3: Piecewise linear interpolant πh f (x) of f (x). Now using (i) and (ii) we get the equation system ? ? ? π f (x ) = c + c x ? c = = f ( x ) h j ?1 0 1 j ?1 j ?1 1 =? ? π f (x ) = c + c x = f (x ) ? c =

h j 0 1 j j 0

f (xj )?f (xj ?1 ) x j ? x j ?1 ?xj ?1 f (xj )+xj f (xj ?1 ) , x j ? x j ?1

Consequently we may write ? ? c = f ( x ) x j + f ( x ) ? x j ?1 j x j ? x j ?1 0 j ? 1 x j ? x j ?1 ? c x = f (x ) ? x + f (x ) x . 1 j ? 1 x j ? x j ?1 j x j ? x j ?1

(2.1.20)

50

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D j = 1, 2 , . . . , N

Hence for xj ?1 ≤ x ≤ xj ,

πh f (x) = c0 + c1 x = f (xj ?1)

xj ? x x ? xj ?1 + f (xj ) xj ? xj ?1 xj ? xj ?1 = f (xj ?1 )λj ?1(x) + f (xj )λj (x),

where λj ?1(x) and λj (x) are the piecewise linear basis functions on Ij :

1 λj ?1 (x) =

x j ?x x j ? x j ?1

λj (x) =

x ? x j ?1 x j ? x j ?1

a

xj ?1

xj

x b

Figure 2.4: Linear Lagrange basis functions for q = 1 on the subinterval Ij .

We generalize the above procedure and introduce Lagrange interpolation bases functions:

2.2

Lagrange interpolation

Consider P q (a, b); the vector space of all polynomials of degree ≤ q on the interval (a, b), and the basis functions 1, x, x2 , . . . , xq for P q . De?nition 11 (Cardinal functions). Lagrange basis is the set of polynomials q { λi } q i=0 ? P (a, b) associated to the (q + 1) distinct points, a = x0 < x1 < . . . < xq = b, in (a, b) and determined by the requirement that λi (xj ) = 1 for i = j , and 0 otherwise, i.e. for x ∈ (a, b), λi (x) = (x ? x0 )(x ? x1 ) . . . (x ? xi?1 ) ↓ (x ? xi+1 ) . . . (x ? xq ) . (2.2.1) (xi ? x0 )(xi ? x1 ) . . . (xi ? xi?1 ) ↑ (xi ? xi+1 ) . . . (xi ? xq )

2.2. LAGRANGE INTERPOLATION By the arrows ↓ , ↑ in (2.2.1) we want to emphasize that λi (x) = x ? xi . Hence does not contain the singular factor xi ? xi λi (xj ) =

51 x ? xj xi ? xj

j =i

(xj ? x0 )(xj ? x1 ) . . . (xj ? xi?1 )(xj ? xi+1 ) . . . (xj ? xq ) = δij , (xi ? x0 )(xi ? x1 ) . . . (xi ? xi?1 )(xi ? xi+1 ) . . . (xi ? xq )

and λi (x), i = 1, 2, . . . , q , is a polynomial of degree q on (a, b) with ? ? 1 i=j λi (xj ) = δij = ? 0 i = j. i = 1, j = 2 ? δ12 = λ1 (x2 ) = (x2 ? x0 )(x2 ? x2 ) =0 (x1 ? x0 )(x1 ? x2 ) (x1 ? x0 )(x1 ? x2 ) = 1. (x1 ? x0 )(x1 ? x2 )

(2.2.2)

Example 13. Let q = 2, then we have a = x0 < x1 < x2 = b, where

i = j = 1 ? δ11 = λ1 (x1 ) =

A polynomial p ∈ P q (a, b) with the values pi = p(xi ) at the nodes xi , i = 0, 1, . . . , q , can be expressed in terms of the corresponding Lagrange basis as p(x) = p0 λ0 (x) + p1 λ1 (x) + . . . + pq λq (x). (2.2.3) Using (2.2.2), p(xi ) = p0 λ0 (xi ) + p1 λ1 (xi ) + . . . pi λi (xi ) + . . . + pq λq (xi ) = pi . Recall that in the previous chapter, introducing examples of ?nite dimensional linear spaces, we did construct Lagrange basis functions for q = 1: λ0 (x) = (x ? ξ1 )/(ξ0 ? ξ1 ) and λ1 (x) = (x ? ξ0 )/(ξ1 ? ξ0 ), for an arbitrary subinterval (ξ0 , ξ1 ) ? (a, b). Below we want to compare the Lagrange polynomial of degree q with another well-known polynomial: namely the Taylor polynomial of degree q . De?nition 12 (Taylor’s Theorem). Suppose that the function f is q +1-times continuously di?erentiable at the point x0 ∈ (a, b). Then, f can be expressed by a Taylor expansion about x0 as f (x) = Tq f (x0 ) + Rq f (x0 ), (2.2.4)

52 where

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

1 1 Tq f (x) = f (x0 ) + f ′(x0 )(x ? x0 ) + f ′′ (x0 )(x ? x0 )2 + . . . + f (q) (x0 )(x ? x0 )q , 2 q! is the Taylor polynomial of degree q , approximating f about x = x0 and Rq f (x) = 1 f (q) (ξ )(x ? x0 )q+1 , (q + 1)! (2.2.5)

is the remainder term, where ξ is a point between x0 and x. For a continuous function f (x) on [a, b], we de?ne the Lagrange interpolation polynomial πq f ∈ P q (a, b), corresponding to the Taylor poynomial Tq f (x). De?nition 13. Let a ≤ ξ0 < ξ1 < . . . < ξq ≤ b, be q +1 distinct interpolation nodes on [a, b]. Then πq f ∈ P q (a, b) interpolates f (x) at the nodes ξi , if πq f (ξi ) = f (ξi), i = 0, 1 , . . . , q (2.2.6)

and the Lagrange’s formula (2.2.3) for πq f (x) reads as πq f (x) = f (ξ0 )λ0 (x) + f (ξ1)λ1 (x) + . . . + f (ξq )λq (x), a ≤ x ≤ b.

Example 14. For q = 1, and considering the whole interval we have only x?a x?b and λb (x) = , thus as we the nodes a and b. Recall that λa (x) = a?b b?a see in the introduction to this chapter π1 f (x) = f (a)λa (x) + f (b)λb (x). Theorem 3. We have the following interpolation error estimates |f (x) ? Tq f (x0 )| = Rq (f ) ≤ 1 (x ? x0 )q+1 · | max f (q+1) (x)|, x∈[a,b] (q + 1)!

q

(2.2.7)

for the Taylor interpolation which is of order q + 1 near x = x0 ; and 1 |f (x) ? πq f (x)| ≤ (q + 1)! |x ? xi | · max |f (q+1) (x)|,

a≤x≤b

i=0

for the Lagrange interpolation error which is of order 1 at each node point x0 , x1 , . . . , xq .

2.2. LAGRANGE INTERPOLATION

53

Proof. The Taylor part is well known from elementary calculus. As for the Lagrange interpolation error we note that at the node points xi we have that f (xi ) ? πq f (xi ) = 0, for i = 0, 1, . . . , q . Since f (xi ) ? πq f (xi ) has q + 1 zeros in [a, b], hence there is a function g (x) de?ned on [a, b] such that f (x) ? πq f (x) = (x ? x0 )(x ? x1 ) . . . (x ? xq )g (x). To determine g (x), we de?ne an auxiliary function ? by setting ?(t) := f (t) ? πq f (t) ? (t ? x0 )(t ? x1 ) . . . (t ? xq )g (x). (2.2.9) (2.2.8)

Note that g (x) is independent of t. Further, the function ?(t) = 0 at the nodes, xi , i = 0, . . . q as well as for t = x, i.e. ?(x0 ) = ?(x1 ) = . . . = ?(xq ) = ?(x) = 0. Thus ?(t) has (q + 2) roots in the interval [a, b]. Now by the Generalized Rolle’s theorem (see below), there exists a point ξ ∈ [a, b] such that ?(q+1) (ξ ) = 0. Taking the (q + 1)-th derivative of the function ?(t), with respect to t, we get ?(q+1) (t) = f (q+1) (t) ? 0 ? (q + 1)!g (x), (2.2.10)

where we use the fact that deg(πq f (x)) = q , (t ? x0 )(t ? x1 ) . . . (t ? xq ) = tq+1 + αtq + . . ., (for some constant α), and g (x) is independent of t. Thus 0 = ?(q+1) (ξ ) = f (q+1) (ξ ) ? (q + 1)!g (x), and we have g (x) = f (q+1) (ξ ) . (q + 1)! (2.2.11)

(2.2.12)

Hence, inserting g from (2.2.12) in (2.2.8), we get the error in Lagrange interpolation as f (q+1) (ξ ) E (x) = f (x) ? πq f (x) = (q + 1)! and the proof is complete. Theorem 4 (Generalized Rolle’s theorem). If a function u(x) ∈ C q+1 (a, b) has (q + 2) roots, x0 , x1 , . . . , xq , x, in a closed interval [a, b], then there is a point ξ ∈ (a, b), generated by x0 , x1 , . . . , xq , x, such that u(q+1) (ξ ) = 0.

q

i=0

(x ? xi ),

(2.2.13)

54

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

In the approximation procedure of solving di?erential equations we encountered matrices with entries being the integrals of products of basis functions. Except some special cases (see calculations for A and Aunif in the previous chapter), these integrals are not elementary and can only be computed approximately. Below we brie?y review some of these numerical integration techniques.

2.3

Numerical integration, Quadrature rule

b

We want to approximate the integral I = a f (x)dx where, on each subinterval, we approximate f using piecewise polynomials of degree d. We denote the approximate value by Id . To this end we assume, without loss of generality, that f (x) > 0 is continuous on the interval [a, b], then the integral b I = a f (x)dx can be interpreted as the area of the domain under the curve y = f (x); limited by x-axis and the lines x = a and x = b. we shall approximate this area using the values of f at certain points as follows.

b b Simple midpoint rule. Uses midpoint a+ of [a, b], i.e. only f a+ . This 2 2 means that f is approximated by the constant function (polynomial of degree b and the area under the curve y = f (x) by 0) P0 (x) = f a+ 2 b

I=

a

f (x)dx ≈ (b ? a)f

a+b . 2

(2.3.1)

This is the general idea of the simple midpoint rule. To prepare for the generalizations, if we denote x0 = a and x1 = b and assume that the length of interval is h, then h (2.3.2) I ? I0 = h · f (a + ). 2 Simple trapezoidal rule. We use two endpoints a and b, i.e, f (a) and f (b). Here f is approximated by the linear function (polynomial of degree 1) P1 (x) passing through the points a, f (a) and b, f (b) and the area under the curve y = f (x) by

b

I=

a

f (x)dx ≈ (b ? a)

f (a) + f (b) . 2

(2.3.3)

2.3. NUMERICAL INTEGRATION, QUADRATURE RULE

55

f (b) f (x + h/2) f (a) a = x0 a + h/2 x b = x1

x1 x0

P0 (x)

Figure 2.5: Midpoint approximationI0 of the integral I =

f (x)dx.

This is the area of the trapezoidal between the lines y = 0, x = a and x = b and under P1 (x), and therefore is refereed as the simple trapezoidal rule. Once again, for the purpose of generalization, we let x0 = a, x1 = b and assume that the length of interval is h, then (2.3.3) can be written as I ? I1 = h · f (a) + h[f (a + h) ? f (a)] f (a) + f (a + h) =h . 2 2 (2.3.4)

P1 (x)

f (b)

f (a) a = x0 x b = x1 = a + h

x1 x0

Figure 2.6: Trapezoidal approximation I1 of the integral I =

f (x)dx.

Simple Simpson’s rule. Here we use two endpoints a and b, and the midpoint

56

a+b , 2

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

b i.e. f (a), f (b), and f a+ . In this case the area under y = f (x) is 2 approximated by the area under the graph of the second degree polynomial b b = f a+ , and P2 (b) = f (b). To deP2 (x) with P2 (a) = f (a), P2 a+ 2 2 termine P2 (x) we may use Lagrange interpolation functions for q = 2: let x0 = a, x1 = (a + b)/2 and x2 = b, then

P2 (x) = f (x0 )λ0 (x) + f (x1 )λ1 (x) + f (x2 )λ2 (x), where ? ? ? λ (x) = ? ? 0 λ1 (x) = ? ? ? ? λ (x) =

2 b b (x?x1 )(x?x2 ) , (x0 ?x1 )(x0 ?x2 ) (x?x0 )(x?x2 ) , (x1 ?x0 )(x1 ?x2 ) (x?x0 )(x?x1 ) . (x2 ?x0 )(x2 ?x1 )

(2.3.5)

(2.3.6)

Thus I=

a

2

b

f (x)dx ≈

f (x)P2 (x) dx =

a i=0

f (xi )

a

λi (x) dx.

(2.3.7)

Now we can easily compute the integrals

b b

λ0 (x) dx =

a a b

b?a λ2 (x) dx = , 6

b

λ1 (x) dx =

a

4(b ? a) . 6

(2.3.8)

Hence I=

a

f (x)dx ≈

b?a [f (x0 ) + 4f (x1 ) + f (x2 )]. 6

(2.3.9)

This is the basic idea behind the simple Simpson’s rule. For the generalization purpose, due to the fact that in this approximation we are using 2 intervals b b ) and ( a+ , b), it is more appropriate to assume an interval of length (a, a+ 2 2 2h. Then, (2.3.9) can be written as

b

I=

a

f (x)dx ≈

h [f (x0 ) + 4f (x1 ) + f (x2 )]. 3

(2.3.10)

2.3. NUMERICAL INTEGRATION, QUADRATURE RULE

57

f (x) f (a) P2 (x)

f (b)

a = x0

a + h/2

x b = x1

x1 x0

Figure 2.7: Simpson’s rule I2 approximating the integral I =

f (x)dx.

Obviously these approximations are less accurate for large interval [a, b] and/or oscillatory functions f . Following the Riemann’s idea we can use these rules, instead of for the whole interval [a, b], for an appropriate partition of [a, b] on each subinterval. Then we get the generalized versions composite rules based on the following algorithm: General algorithm. To approximate the integral

b

I=

a

f (x)dx,

we use the following steps (i) Divide the interval [a, b], uniformly, into N subintervals a = z0 < z1 < z2 < . . . < zN ?1 < zN = b. (ii) Write the integral as

b z1 zN N zk

(2.3.11)

f (x)dx =

a z0

f (x) dx + . . . +

z N ?1

f (x) dx =

k =1 z k ?1

f (x) dx. (2.3.12)

(ii) For each subinterval Ik := [zk?1 , zk ], k = 1, 2, . . . , N , apply the very same integration rule. Then we have the following generalizations.

58

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

(M) Composite midpoint rule: approximates f by the constants (that are the values of f at the midpoint of the subinterval) on each subinterval: Let b?a zk?1 + zk h = |IN | = , x ?k = , k = 1, 2, . . . , N, N 2 then using the simple midpoint rule for the interval Ik := [zk?1 , zk ], we have zk zk f (? xk ) dx = hf (? xk ). (2.3.13) f (x) dx ≈

z k ?1 z k ?1

Summing over k = 1, 2, . . . , N , we get

b a N

f (x)dx ≈

hf (? xk ) = h[f (? x1 ) + . . . + f (? xN )] := mN .

k =1

(2.3.14)

(T) Composite trapezoidal rule: Simple trapezoidal rule on Ik yields

zk z k ?1 zk

f (x) dx ≈

N

f (? xk ) dx =

z k ?1

h [f (zk?1 ) + f (zk )]. 2

(2.3.15)

Summing over k = 1, 2, . . . , N , we have

b a

f (x)dx ≈

k =1

h [f (zk?1) + f (zk )] 2

(2.3.16)

h = [f (z0 ) + 2f (z1 ) + . . . + 2f (zk?1) + f (zk )] := tN . 2 (S) Composite Simpson’s rule: (Quadratic approximation on each subinterval). Recall that this corresponds to a quadrature rule based on piecewise quadratic interpolation using the endpoints and midpoint of each subinterval. Thus, this time we use the simple Simpson’s rule on z +z k and zk : each subinterval Ik = [zk?1 , zk ] using the points zk?1 , k?1 2

zk z k ?1

f (x) dx ≈

zk?1 + zk h f (zk?1 ) + 4f + f (zk ) . 3 2

(2.3.17)

To proceed we use the following identi?cation on each subinterval Ik , k = 1, . . . , k , x2k?2 = zk?1 , x2k?1 = zk?1 + zk := z ?k , 2 x2k = zk . (2.3.18)

2.3. NUMERICAL INTEGRATION, QUADRATURE RULE

59

z0

z ?1

z1

z2k?1

z ?k

zk

zN = b

a = x0

x1

x2

x2k?2

x2k?1

x2k

x2N

Figure 2.8: Identi?cation of subintervals for composite Simpson’s rule

Thus summing (2.3.17) over k and using the above identi?cation, we ?nally obtain the composite Simpson’s rule viz,

b a N

f (x)dx ≈ =

k =1 N

h zk?1 + zk f (zk?1 ) + 4f + f (zk ) 3 2 h f (x2k?2 ) + 4f (x2k?1 ) + f (x2k ) 3 (2.3.19)

k =1

=

h f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) + 2f (x4 ) 3 + . . . + 2f (x2N ?2 ) + 4f (x2N ?1 ) + f (x2N ) := SN .

Here is the starting procedure where the numbers in the brackets indicate the actual coe?cient on each subinterval. For instance the end of interval 1 :z1 = x2 coincides with the start of interval 2, yielding to the add-up [1] + [1] = 2 as the coe?cient of f (x2 ). A resonance which is repeated for each interior node zk . k = 1, . . . , N ? 1. Remark 9. The rules (M), (T) and (S) use values of the function at equally spaced points. These are not always the best approximation methods. Below we introduce an optimal method:

60

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D [4] [4]

[1] + [1] [1] x0 x1 x2 x3 [1] x4

Figure 2.9: Identi?cation of subintervals for composite Simpson’s rule

2.3.1

Gauss quadrature rule

This is to choose the points of evolution in an optimal manner, not at equally spaced points. We demonstrate this rule through an example viz: Problem: Choose the nodes xi ∈ [a, b], and coe?cients ci , 1 ≤ i ≤ n such that, for an arbitrary function f , the following error is minimal:

b a n

f (x)dx ?

ci f (xi ) for an arbitrary function f (x).

i=1

(2.3.20)

Solution. There are 2n unknowns in (2.3.20) consisting of n nodes xi and n coe?cients ci . Therefore we need 2n equations. Thus if f were a polynomial, optimal choice of our parameters yields a quadrature rule (2.3.20) which is exact for polynomials of degree ≤ 2n ? 1. Example 15. Let n = 2 and [a, b] = [?1, 1]. Then the coe?cients are c1 and c2 and the nodes are x1 and x2 . Thus optimal choice of these 4 parameters should yield that the approximation

1 ?1

f (x)dx ≈ c1 f (x1 ) + c2 f (x2 ),

(2.3.21)

is indeed exact for f (x) replaced by polynomials of degree ≤ 3. Thus we replace f by a polynomial of the form f (x) = Ax3 + Bx2 + Cx + D and

2.3. NUMERICAL INTEGRATION, QUADRATURE RULE

61

require equality in (2.3.21). Thus to determine the coe?cients c1 , c2 and the nodes x1 , x2 , in an optimal way, it su?ces to use the above approximation as equality when f is replaced by the basis functions for polynomials of degree 3: 1, x, x2 and x3 . Consequently we get the equation system

1 ?1 1

1dx = c1 + c2 and we get [x]1 ? 1 = 2 = c1 + c2 xdx = c1 · x1 + c2 · x2 and

x2 1 = 0 = c1 · x1 + c2 · x2 2 ?1 ?1 1 x3 1 2 2 2 x2 dx = c1 · x2 + c · x and = = c1 · x2 2 1 2 1 + c2 · x2 3 3 ? 1 ?1 1 x4 1 3 3 3 3 = 0 = c1 · x3 x dx = c1 · x1 + c2 · x2 and 1 + c2 · x2 , 4 ?1 ?1

(2.3.22)

Thus the approximation

1

which, although nonlinear, has the solution presented below: ? ? ? ? ? ? c + c = 2 c1 = 1 1 2 ? ? ? ? ? ? ? ? ? c x +c x =0 ? c =1 1 1 2 2 2 =? √ 2 ? ? ? ? c1 x2 + c2 x2 =3 x1 = ? 33 ? ? 1 2 ? ? ? ? √ ? ? ? c x3 + c x3 = 0 ? x = 3 1 1 2 2 2 3 3 f (x)dx ≈ c1 f (x1 ) + c2 f (x2 ) = f ? +f 3 ?1 is exact for polynomials of degree ≤ 3.

1

(2.3.23)

√

√

3 , 3

(2.3.24)

Example 16. Let f (x) = 3x2 + 2x + 1. Then ?1 (3x2 + 2x + 1)dx = √ √ [x3 + x2 + x]1 ?1 = 4, and we can easily check that f (? 3/3) + f ( 3/3) = 4, which is also the exact value of the integral. Generalized Gauss quadrature. To generalize Gauss quadrature rule Legendre polynomials are used: Choose {Pn }∞ n=0 such that (1) For each n, Pn is a polynomial of degree n.

62

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

1 ?1

(2) Pn ⊥ Pm if m < n ??

Pn (x)Pm (x)dx = 0

The Legendre polynomial can be obtained through formula: Pk (x) = (?1)k dk k (x (1 ? x)k ), dxk or Pn (x) = 2 2n n!dxn (x2 ? 1)n ,

Here are a few ?rst Legendre polynomials: P0 (x) = 1, P1 (x) = x, 3 1 P2 (x) = x2 ? , 2 2 5 3 P3 (x) = x3 ? x, . . . , 2 2

The roots of Legendre polynomials are distinct, symmetric and the correct choices as quadrature points, i.e. they are giving the points xi , 1 ≤ i ≤ n, as the roots of the n-th Legendre polynomial. (p0 = 1 is an exception). Example 17. Roots of the Legendre polynomial as quadrature points: P1 (x) = x = 0. √ 3 gives x1,2 = ± . (compare with the result above). 3 gives x1 = 0, x2,3 = ± 3 . 5

1 3 P2 (x) = x2 ? = 0, 2 2

5 3 P3 (x) = x3 ? x = 0, 2 2

Theorem 5. Suppose that xi , i = 1, 2, . . . , n, are roots of n-th Legendre polynomial Pn and that

1 n

ci =

?1 j =1 j =i

x ? xj dx, where xi ? xj

n

j =1 j =i

x ? xj xi ? xj

1

is the Lagrange basis.

n

If f (x) is a polynomial of degree < 2n, then

?1

f (x)dx ≡

ci f (xi ).

i=1

Proof. Consider a polynomial R(x) of degree < n. Rewrite R(x) as (n ? 1) Lagrange polynomials with nodes at the roots of the n-th Legendre polynomial Pn . This representation of R(x) is exact, since the error is E (x) = 1 (x?x1 )(x?x2 ) . . . (x?xn )R(n) (ξ ), n! where R(n) (ξ ) ≡ 0. (2.3.25)

2.3. NUMERICAL INTEGRATION, QUADRATURE RULE

n n

63

Further we have R(x) =

i=1 j =1 j =i 1 1 n

x ? xj R(xi ), so that xi ? xj

n

R(x)dx =

?1 ?1 i=1 j =1 j =i 1 n ?1 j =1 j =i

x ? xj R(xi ) dx xi ? xj x ? xj dx R(xi ). xi ? xj

n

(2.3.26)

=

i=1

Moreover

1

R(x)dx =

?1 i=1

ci R(xi )

(2.3.27)

Now consider a polynomial, P (x), of degree < 2n. Dividing P (x) by the n-th Legendre polynomial Pn (x), we get P (x) = Q(x) × Pn (x) + R(x), and

1 1

degQ(x) < n,

1

degR(x) < n, (2.3.28) (2.3.29)

P (x)dx =

?1 ?1

Q(x)Pn (x)dx +

?1

R(x)dx.

Since Q(x) ⊥ Pn (x), ?Q(x) with degree< n, thus using (2.3.28) it follows that

1 ?1 1 1

Q(x)Pn (x)dx = 0 =?

P (x)dx =

?1 ?1

R(x)dx.

(2.3.30)

Then xi ’s are roots of Pn (x), thus Pn (xi ) = 0 and we can use (2.3.28) to write P (xi ) = Q(xi )Pn (xi ) + R(xi ) = R(xi ). (2.3.31) Now using (2.3.27) we obtain that

1 1 n n

P (x)dx =

?1 ?1

R(x)dx =

i=1 1

ci R(xi ) =

i=1 n

ci P (xi ).

(2.3.32)

Summing up: P (x)dx =

?1

ci P (xi ),

i=1

(2.3.33)

and the proof is complete.

64

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

2.4

Exercises

b?x b?a

Problem 11. Use the expressions λa (x) = λa (x) + λb (x) = 1,

and λb (x) =

x ?a b?a

to show that

aλa (x) + bλb (x) = x.

Give a geometric interpretation by plotting, λa (x), λb (x), λa (x) + λb (x), aλa (x), bλb (x) and aλa (x) + bλb (x). Problem 12. Let f : [0.1] → R be a Lipschitz continuous function. Determine the linear interpolant πf ∈ P (0, 1) and plot f and πf in the same ?gure, when (a) f (x) = x2 , (b) f (x) = sin(πx). Problem 13. Determine the linear interpolation of the function f (x) = π 1 (x ? π )2 ? cos2 (x ? ), 2 π 2 ?π ≤ x ≤ π.

where the interval [?π, π ] is divided to 4 equal subintervals. Problem 14. Assume that w ′ ∈ L1 (I ). Let x, x ? ∈ I = [a, b] and w (? x) = 0. Show that |w (x)| ≤ |w ′|dx. (2.4.1)

I

Problem 15. Assume that v interpolates ?, at the points a and b.

? v

a

h

b

x

2.4. EXERCISES Show, using (2.4.1) that (i) (ii) (iii) (iv)

I

65

|(? ? v )(x)| ≤ |(? ? v )′ (x)| ≤

I

I

|(? ? v )′ | dx, |(? ? v )′′ | dx =

I

I I

|?′′ | dx,

max |? ? v | ≤ max |h2 ?′′ |, |? ? v | dx ≤

I 2 I ′′

|h2 ?′′ | dx,

I 1/2

(v) where

??v w

I

≤ h?

I

and

h?2 (? ? v )

I

≤ ?′′ I ,

=

w 2 dx

is the L2 (I )-norm.

Problem 16. Use, in the above problem ?(b) ? ?(a) 1 v = = h h

′ b

?′ dx (?′ is constant on I ),

a

and show that (vi) (vii)

I

|(? ? v )(x)| ≤ 2

I

|?′ | dx,

I

h?1 |? ? v | dx ≤ 2

|?′ | dx

and

h?1 (? ? v ) ≤ 2 ?′ I .

Problem 17. Let now v (t) be the constant interpolant of ? on I .

? v

a

x

b

66 Show that

CHAPTER 2. POLYNOMIAL INTERPOLATION IN 1D

I

h?1 |? ? v | dx ≤

I

|?′ | dx.

(2.4.2)

Problem 18. Show that P q (a, b) := {p(x)|p(x) is a polynomial of degree ≤ q }, is a vector space but P q (a, b) := {p(x)|p(x) is a polynomial of degree = q }, is not! a vector space. Problem 19. Compute formulas for the linear interpolant of a continuous function f through the points a and (b+a)/2. Plot the corresponding Lagrange basis functions. Problem 20. Prove the following interpolation error estimate: 1 ||Π1 f ? f ||L∞ (a,b) ≤ (b ? a)2 ||f ′′ ||L∞ (a,b) . 8 Hint: Use Theorem 5.1 from PDE Lecture Notes. Problem 21. Compute and graph π4 e?8x e

?8x2

2

Problem 22. Write down a basis for the set of piecewise quadratic polyno(2) mials Wh on a partition a = x0 < x1 < x2 < . . . < xm+1 = b of (a, b) into subintervals Ii = (xi?1 , xi ), where Wh = {v : v |Ii ∈ P q (Ii ), i = 1, . . . , m + 1} Problem 23. Prove that

x1 (q )

at 5 equally spaced points in [?2, 2].

on [?2, 2], which interpolates

f′

x0

x1 + x0 2

x?

x1 + x0 dx = 0. 2

Problem 24. Prove that x1 x1 + x0 f (x) dx ? f (x ? x0 ) 2 x0 x1 1 x1 + x0 ′′ ≤ max |f | x? 2 [x0 ,x1 ] 2 x0

2

dx ≤

1 (x1 ? x0 )3 max |f ′′ |. [x0 ,x1 ] 24

Hint: Use Taylor expansion of f about x =

x 1 +x 2 . 2

Chapter 3 Linear System of Equations

This chapter is devoted to numerical solution of the linear system of equations of type Ax = b ? x = A?1 b. To this approach we shall review the wellknown direct method of Gauss elimination and then continue with some more e?cient iterative methods. Numerical linear algebra is undoubtedly the most applied tool in the computational aspects of almost all disciplines.

3.1

Direct methods

? ? ? a11 x1 + a12 x2 + . . . + a1n xn = b1 ? ? ? ? ? a x + a x + ...+ a x = b 21 1 22 2 2n n 2 i = 1, . . . , n, or ? ? ... ? ? ? ? ? a x + a x + ...+ a x = b . n1 1 n2 2 nn n n

Consider the general form of an n × n linear system of equations given by

n

Ax = b ?

aij xj = bi ,

j =1

We introduce the enlarged n × (n + 1) coe?cient matrix A by ? ? a11 a12 . . . a1n b1 ? ? ? ? ? a21 a22 . . . a2n b2 ? ?, ? A := ? ? ? ... ... ... ... ... ? ? ? an1 an2 . . . ann bn 67

(3.1.1)

68

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS

where the coe?cient matrix A is the ?rst n columns in A. Note that to solve the equation system Ax = b t is a bad idea to calculate A?1 and then multiply by b. However, if A is an upper (or lower) triangular, i.e. aij = 0 for i > j (or i < j ), and A is invertible, then we can solve x using the back substitution method: ? ? ? a11 x1 + a12 x2 + . . . + a1n xn = b1 ? ? ? ? ? ? a22 x2 + . . . + a2n xn = b2 ? ? ? ? ? ... ... ... (3.1.2) ? ? ... ... ... ? ? ? ? ? ? an?1,n?1xn?1 + an?1,n xn = bn?1 ? ? ? ? ? ann xn = bn , ? 1 ? ? x1 = [b1 ? a12 x2 ? . . . ? a1n xn ] ? ? a ? 11 ? ? ? ? ... ... ... ? ? ? ? ? ? ... ... ... (3.1.3) 1 ? [ b ? a x ] x = ? n ? 1 n ? 1 ,n n n ? 1 ? an?1,n?1 ? ? ? ? ? ? ? ? ? ? bn ? ? xn = . ann ? Number of operations. Additions and subtractions are not considered as time consuming operations, therefore we shall count only the number of multiplications and divisions. ? The number of multiplications to solve xn from (3.1.3) are zero and the number of divisions is one. ? To solve xn?1 we need one multiplication and one division. ? To solve x1 we need (n ? 1) multiplication and one division. Thus to solve the linear system of equations given by (3.1.2) we shall need 1 + 2 + . . . + (n ? 1) = yields

n(n ? 1) n2 := + Q(n), 2 2 multiplications, where Q(n) is a remainder of order n, and n divisions.

3.1. DIRECT METHODS

69

?Gaussian elimination method. The Gauss elimination method is based on the following obvious facts expressing that: a linear system is not changed under elementary row operations. These are (i) interchanging two equations (ii) adding a multiple of one equation to another (iii) multiplying an equation by a nonzero constant. Before continuing with the Gauss elimination procedure we recall the simple 3 × 3 dimensional uper triangular matrix U , lower triangular matrix L and diagonal matrix D . ? ? ? ? ? ? a 0 0 a 0 0 a b c ? ? ? ? ? ? ? ? ? ? ? ? D = ? 0 d 0 ?. L = ? g d 0 ?, U = ? 0 d e ?, ? ? ? ? ? ? 0 0 f h i f 0 0 f

The Gauss elimination procedure relay on the elementary row operations and converts the coe?cient matrix of the linear equation system to an upper triangular matrix. To this end, we start from the ?rst row of the coe?cient matrix of the equation system and using elementary row operations eliminate the elements ai1 , i > 1, under a11 (make ai1 = 0). ? with The equation system corresponding to this newly obtained matrix A elements a ?ij , a ?i1 = 0, i > 1, has the same solution as the original one. We repeat the same procedure of the elementary row operations to eliminate the ?. Continuing in this way, we thus elements ai2 , i > 2, from the matrix A obtain an upper triangular matrix U with corresponding equation system equivalent to the original system (has the same solution). Below we shall illustrate this procedure through an example: Example 18. Solve the equation system: ? ? ? 2x + x2 + x3 = 2 ? ? 1 4x1 ? x2 + 3x3 = 0 ? ? ? ? 2x + 6x ? 2x = 10. 1 2 3

(3.1.4)

70

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS

In the coe?cient matrix: 2 1 1 | 2 ? ? A = ? 4 ?1 3 | 0 ? 2 6 ?2 | 10 ? ? ? ? ?, ?

(3.1.5)

we have that a11 = 2, a21 = 4, and a31 = 2. We introduce the multipliers mi1 , i > 1 by letting a21 4 a31 2 m21 = = =2 m31 = = = 1. (3.1.6) a11 2 a11 2 Now we multiply the ?rst row by m21 and then subtract it from row 2 and replace the result in row 2: ? ? ? ? 2 1 1 | 2 ·(?2) 2 1 1 | 2 ? ? ? ? ? ? ? ? (3.1.7) =? ? 0 ?3 1 | ?4 ? ? 4 ?1 3 | 0 ? ? ? ? ? 2 6 ?2 | 10 2 6 ?2 | 10

Similarly we multiply the ?rst row by m31 = 1 and subtract it from row 3 to get ? ? 2 1 1 | 2 ? ? ? ? (3.1.8) ? 0 ?3 1 | ?4 ? . ? ? 0 5 ?3 | 8 In this setting we have a ?22 = ?3 and a ?32 = 5, and ? ? 2 1 1 ? ? ? ?=? A ? 0 ?3 1 ? . ? ? 0 5 ?3

(3.1.9)

? by m32 Now let m32 = a ?32 /a ?22 = ?5/3, then multiplying the second row in A and subtracting the result from row 3 yields ? ? 2 1 1 | 2 ? ? ? ? (3.1.10) 0 ? 3 1 | ? 4 ? ?, ? ? 4 4 0 0 ? | 3 3

3.1. DIRECT METHODS where we have obtained the upper triangular matrix ? ? 2 1 1 ? ? ? ? U = ? 0 ?3 1 ? . ? ? 4 0 0 ? 3 The new equivalent equation system is now ? ? 2x1 + x2 + x3 = 2 ? ? ? ?3x2 + x3 = ?4 ? ? 4 ? 4 ? ? x3 = 3 3

71

(3.1.11)

(3.1.12)

with the obvious solution x1 = 1, x2 = 1 and x3 = ?1 which, as we can verify is also the solution of the original equation system (3.1.4) De?nition 14. We de?ne the lower triangular matrices: ? ? ? ? ? ? 1 0 0 1 0 0 1 0 0 ? ? ? ? ? ? ? ? ? ? ? ? and L = L1 = ? ?m21 1 0 ? , L2 = ? 0 ? m21 1 0 ? . 1 0 ? ? ? ? ? ? ? m31 m32 1 0 ?m32 1 ?m31 0 1

The matrices L1 , L2 and L3 are unite (ones on the diagonal) lower triangular 3 × 3-matrices with the property that

1 ?1 L = (L2 L1 )?1 = L? 1 L2 ,

and

A = LU.

(3.1.13)

? LU factorization of the matrix A We generalize the above procedure for the 3 × 3 linear system of equations to n × n. We have then A = LU , where L is a unite lower triangular matrix and U is an upper triangular matrix obtained from A by Gauss elimination. To solve the system Ax = b we let now y = Ux, and ?rst solve Ly = b by forward substitution (from the ?rst row to the last) and obtain the vector y , then using y as the known right hand side ?nally we solve Ux = y by backward substitution (from the last row to the ?rst) and get the solution x.

72

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS

Example 19. Following our previous example m21 = 2, m31 = 1 and m32 = ?5/3, consequently ? 1 0 0 ? ? 1 0 0 ? ? ? ? ? ? ? 1 0 0 ?

Thus

? ? L1 = ? ?2 1 0 ? ?1 0 1 ?

? ? ?, ?

? ? L2 = ? 0 1 0 ? 5 0 1 3 ??

and

? ? ? ? L = ? 2 1 0 ?. ? ? 5 1 ? 1 3 ?

which corresponds to the last (third) elementary row operation performed in the previous example. In general case we have the following result: Proposition 2. The n × n unit lower triangular L is given by L = (Ln?1 Ln?2 . . . L1 )?1 , where Li , i = 1, . . . , n ? 1 are the corresponding n × n row-operation matrices, viz example above. For n = 3 we have (L2 L1 )?1 = L, where ? ? 1 0 0 ? ? ? ? L = ? m21 1 0 ? , ? ? m31 m32 1 and mij are the multipliers de?ned above.

which corresponds to the elimination. Further ? 1 0 0 ? ? L2 L1 A = ? 0 1 0 ? 5 0 1 3

2 1 1 1 0 0 ?? ? ?? ? L1 A = ? ?2 1 0 ? ? 4 ?1 3 ?? ? 2 6 ?2 ?1 0 1 ??

?rst two elementary row operations in Gaussian 2 1 1 ? ? 2 1 1 ?

2 1 1 ? ? ? ? ? = ? 0 ?3 1 ? ? 0 5 ?3

?

? ? ? ? = A, ?

?? ?? ? ? 0 ?3 1 ?? 0 5 ?3

? ? ? ? ? = ? 0 ?3 1 ? ? 4 0 0 ? 3

? ? ? = U, ?

3.1. DIRECT METHODS

73

Thus Ax = b ?? (LU )x = b ?? L(Ux) = b. As we outlined above we let y = Ux and solve Ly = b to obtain y . Then with such obtained y we solve x from Ux = y . We illustrate this procedure through our example: Example 20. In our example ? 1 0 ? ? L=? 2 1 ? 5 1 ? 3 we have that ? 0 ? ? and 0 ? ? 1

? Ly = b yields ? 1 0 ? ? 1 ? 2 ? 5 1 ? 3

Using backward substitution, we get the solution viz x1 = 1, x2 = 1, x3 = ?1.

? Ux = y yields ? 2 1 1 ? ? 1 ? 0 ?3 ? 4 0 0 ? 3

Using forward substitution we get y1 = 2, y2 = ?4, y3 = 4/3. Further with ? ? ? ? 2 1 1 2 ? ? ? ? ? ? ? ? U = ? 0 ?3 and y = 1 ? ? ?4 ? , ? ? ? ? 4 4 0 0 ? 3 3 ? ? 2x1 + x2 + x3 = 2 ? ? ? ? ? ?? ? ? ? ?? ? ? 3x2 + x3 = ?4 ? ? x2 ? = ? ?4 ? ?? ? ? ? ?? ? ? 4 4 4 ? ? ? x3 = . x3 3 3 3 ?? x1 ? ? 2 ?

the system of equations ? ? ? ? ?? ? y = 2 0 ? 2 y ? ? ? 1 ?? 1 ? ? ? ? ? ?? 2y 1 + y 2 = 0 0 ? ? y2 ? = ? 0 ? ?? ? ? ? ? ?? ? ? ? y1 ? 5 y2 + y3 = 10. 1 10 y3 3

? ? ? ? b = ? 0 ?. ? ? 10

?

2

?

74

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS

Theorem 6 (Cholesky’s method). Let A be a symmetric matrix, (aij = aji ), then the following statements are equivalent: (i) A is positive de?nite. (ii) The eigenvalues of A are positive. (iii) Sylvester’s criterion det(?k ) > 0 for k = 1, 2, . . . , n, where ? ? a . . . a1k ? ? 11 ? ? ?k = ? . . . . . . . . . ? . ? ? ak1 . . . akk (iv) A = LLT where L is lower triangular and has positive diagonal elements. (Cholesky’s factorization) We do not give a proof of this theorem. The interested reader is referred to literature in algebra and matrix theory.

3.2

Iterative method

Instead of solving Ax = b directly, we consider iterative solution methods based on computing a sequence of approximations x(k) , k = 1, 2, . . . such that

k →∞

lim x(k) = x

or

k →∞

lim x(k) ? x = 0,

for some norm.

Thus consider the general n × n linear system of equations Ax = b where both the coe?cient matrix A and the vector b have real entries, ? ? ? a11 x1 + a12 x2 . . . +a1n xn = b1 ? ? ? ? ? ... ... ... ... ... (3.2.1) Ax = b ?? ? ? a x + . . . . . . + a x = b n ? 1 , 1 1 n ? 1 ,n n n ? 1 ? ? ? ? ? a x+ . . . . . . +a x =b .

n1 1 nn n n

3.2. ITERATIVE METHOD

75

For the system (3.2.1) we shall introduce the two main iterative methods. ? Jacobi iteration: Assume that aii = 0, then ? 1 ? ? x1 = ? [a12 x2 + a13 x3 + . . . + a1n xn ? b1 ] ? ? a11 ? ? 1 xn?1 = ? [an?1,1 x1 + an?1,2 x2 + . . . + an?1,n xn ? bn?1 ] an?1,n?1 ? ? ? 1 ? ? ? xn = ? [an1 x1 + an2 x2 + . . . + an,n xn ? bn ]. ann Given an initial approximation for the solution:

(0) (0)

x(0) = (x1 = c1 , x2 = c2 , . . . , x(0) n = c n ), the iteration steps are given by ? 1 (k +1) (k ) (k ) ? k) ? x1 = ? [a12 x2 + a13 x3 + . . . + a1n x( ? n ? b1 ] ? a ? 11 ? ? ? ? ? ? ? 1 (k ) (k ) (k +1) k) x2 = ? [a21 x1 + a23 x3 + . . . + a2n x( n ? b2 ] ? a 22 ? ? ? ? ... ? ? ? ? 1 ? (k ) (k ) (k ) k +1) ? ? x( =? [an1 x1 + an2 x2 + . . . + an,n?1 xn?1 ? bn ] n ann ? ? ?

n j =1 aij xj

Or in compact form in Jacobi coordinates:

= bi ?? aii xi = ?

(k ) n j =1 aij xj j =i

n j =1 j =i

aij xj + bi , (3.2.2)

Convergence criterion Jacobi gives convergence to the exact solution if A is diagonally dominant.

n

k +1) ? ? aii x( =? i

+ bi .

|aii | >

j =1 j =i

|aij | i = 1, 2, . . . , n.

76

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS ? ?

Note, the Jacobi method needs less operations than Gauss elimination. ? ? ? ? 2 ?1 1 ? and b = ? ?. Example 21. Solve Ax = b where A = ? ?1 2 1 A is diagonally dominant and the matrix equation Ax = b is equivalent to the linear equation system ? ? 2x ? x = 1 1 2 (3.2.3) ? ?x + 2x = 1.

1 2 (0) (0)

4 2 1 ? ? ? ? Problem 25. Show that A = ? 1 5 1 ? is diagonally dominant. ? ? 0 1 3

We choose zero initial values for x1 and x2 , i.e. x1 = 0 and x2 = 0 and build the Jacobi iteration system ? ? 2x(k+1) = x(k) + 1 1 2 (3.2.4) ( k +1) ( k) ? 2x = x + 1,

2 1

where k is the iteration step. Then we have ? ? 2x(1) = x(0) + 1 1 2 ? 2x(1) = x(0) + 1 2 1

with the solution

Continuing we have obviously lim xi

k →∞

In the next iteration step: ? ? ? ? x(2) = 3 ? 2x(2) = x(1) + 1 ? 2x(2) = 1 + 1 1 1 1 2 2 4 ? ? ? x(2) = 3 ? 2x(1) = x(1) + 1 ? 2x(2) = 1 2 1 2 2 2 4

(k )

? 1 ? ? x(1) 1 = 2 1 ? ? x(1) 2 = . 2

(3.2.5)

(3.2.6)

= xi , i = 1, 2, where x1 = x2 = 1.

(k )

Below we have a few ?rst iterations giving the corresponding x1 values

and x2

(k )

3.2. ITERATIVE METHOD k 0 1 2 3 x1

(k )

77 x2

(k )

0 1/2 3/4 7/8

0 1/2 3/4 7/8

∞

Now if we use the maximum norm: ek e0 e1

∞ ∞ (0) (0)

:= max |xi ? xi |, then

i=1,2

(k )

= max(|x1 ? x1 |, |x2 ? x2 |) = max = max(|x1 ? x1 |, |x2 ? x2 |) = max = max(|x1 ? x1 |, |x2 ? x2 |) = max = max(|x1 ? x1 |, x2 ? x2 |) = max

(3) (3) (2) (2) (1) (1)

0?1 , 0?1

=1 = 1 2 1 4 1 8

1 1 ?1 , ?1 2 2 3 3 ?1 , ?1 4 4

e2

∞

=

e3

∞

7 7 ?1 , ?1 8 8

=

1 ek ∞ , where ek is the error for step k, k ≥ 0. In this way ek+1 ∞ = 2 Iterating we get that for the k -the Jacobi iteration the convergence rate is 1 k : 2 ek

∞

=

1 ek?1 2

∞

=

1 2

2

ek?2

∞

= ... =

1 2

k

e0

∞

=

1 2

k

.

78

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS

? Gauss-Seidel Method Give an initial approximation of the solution x = x1 , x2 , . . . , x(0) , n then using the fact that the ?rst row in the k -th Jacobi iteration gives x1 (k +1) (k +1) and in the i +1-th row we have already computed values for x1 , . . . , xi on the right hand sides of the ?rst i rows. The idea with the Gauss-Seidel method is that, in the same iteration step, simultaneously use this computed values. More speci?cally the Gauss-Seidel iteration steps are given by: ? ?1 (k ) (k ) (k +1) (k ) ? [a12 x2 + a13 x3 + . . . + a1n xn ? b1 ] x = ? 1 ? ? a 11 ? ? ? ? ? ? ? ? (k+1) ?1 ? (k +1) (k ) (k ) ? x2 = [a21 x1 + a23 x3 + . . . + a2n xn ? b2 ] ? ? ? a 22 ? ... ? ? ?1 ? (k +1) (k +1) (k +1) k) ? ? xn?1 = [a(n?1),1 x1 + . . . + a(n?1),n?2 xn?2 + a(n?1),n x( ? n ? bn?1 ] ? a ? n?1,n?1 ? ? ? ? ? ? ? ? ? ? x(k+1) = ?1 [an1 x(k+1) + an2 x(k+1) + . . . + an,n?1 x(k+1) ? bn ]. 1 2 n?1 n ann Or in a compact way in Gauss-Seidel coordinates.

n i n (k +1) (0) (0)

Ax = b ??

j =1

aij xj = bi ??

aij xj +

j =1 j =1+1

aij xj = bi .

(3.2.7)

Therefore the iterative forms for the Gauss-Seidel is given by ? (k ) ? i a x(k+1) = ? n j =1 ij j j =i+1 xj + bi ?? (k ) ? a x(k+1) = ? i?1 a x(k+1) ? n a x +b.

ii i j =1 ij j j =i+1 ij j i

(3.2.8)

Example 22. We consider the same example as above: Ax = b with ? ? ? ? ? ? 2 ?1 x 1 ?, x = ? 1 ?, b = ? ?. A=? ?1 2 x2 1

3.2. ITERATIVE METHOD Recall the Jacobi iteration system ? ? 2x(k+1) = x(k) + 1 1 2 ? 2x(k+1) = x(k) + 1. 2 1

79

(3.2.9)

We choose the same initial values for x1 and x2 as in the Jacobi iterations, (0) (0) i.e. x1 = 0, and x2 = 0. Now the ?rst equation in (3.2.10): 1 (1) (0) (1) 2x1 = x2 + 1 =? x1 = . 2 Inserting this value of x1 =

(1) (1) (1) 1 2

The corresponding Gauss-seidel iteration system reads as follows ? ? 2x(k+1) = x(k) + 1 1 2 ( k +1) ( k +1) ? 2x =x +1

2 1

(3.2.10)

in the second equation in (3.2.10) yields

(1)

2x2 = x1 + 1 =? 2x2 =

3 1 (1) + 1 =? x2 = . 2 4

Below we list a few ?rst iteration steps for this Gauss-Seidel approach: k 0 1 2 3 Obviously lim x1 = lim x2

k →∞ k →∞ (k ) (k )

x1

(k )

x2

(k )

0 1/2 7/8

0 3/4 15/16

31/32 63/64 = 1. Now with ek

∞

get the successive iteration errors: e1

∞

= max |xi ? xi |, we

i=1,2

(k )

= max(|x1 ? x1 |, |x2 ? x2 |) = max = max 7 15 ?1 , ?1 8 16 1 == , 8 e3

(1)

(1)

3 1 ?1 , ?1 2 4 = max

=

1 2

e2

∞

∞

1 1 1 = . , 32 64 32

80

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS

1 Thus for the Gauss-Seidel iteration ek+1 ∞ = ek ∞ , where ek is the error 4 for step k , and hence we can conclude that the Gauss-Seidel method converges faster than the Jacobi method: ek

∞

=

1 ek?1 4

∞

=

1 4

2

ek?2

∞

= ··· =

1 4

k

e0

∞

=

1 4

k

.

? The successive over-relaxation method (S.O.R.). The S.O.R. method is a modi?ed version of the Gauss-Seidel iteration. The iteration procedure is given by

(k +1) xi

= (1 ?

(k ) ω )xi

ω bi ? + aii

i?1 j =1

n (k +1) aij xj

?

aij xj

j =i+1

(k )

(3.2.11)

For ω > 1 the method is called an over-relaxation method and if 0 < ω < 1, it is referred as an under-relaxation method. S.O.R. coordinates:

(k +1) aii xi

=

(k ) aii xi

?ω

i?1 j =1

n (k +1) aij xj

+

j =i+1

aij xj ? bi

(k )

(3.2.12)

? Abstraction of iterative methods In our procedures we have considered Ax = b and x = Bx + C as equivalent linear equation systems, where B is the iteration matrix and xk+1 = Bxk + C . Potential advantages of iteration methods over direct methods: they are (i) Faster (depends on B , accuracy is required) (ii) Less memory is required (Sparsity of A can be preserved.) Questions: (Q1) For a given A, what is a good choice for B ? (Q2) When does xk → x? (Q3) What is the rate of convergence? The error at step k is ek = xk ? x and that of step (k + 1) is ek+1 = xk+1 ? x. Then we have ek+1 = xk+1 ? x = (Bxk + C ) ? (Bx ? C ) = B · (xk ? x) = Bek .

ek

Iterating, we have ek = Bek?1 = B · Bek?2 = B · B · Bek?3 = B 4 ek?4 = . . . = B k ek?k = B k e0 .

3.2. ITERATIVE METHOD Thus we have shown that ek ? 0 ... ... ? ? ? a21 0 ... L=? ? ? ... ... ... ? an1 . . . an,n?1 = B k e0 . Let now ? ? 0 0 ? ? ? ? ? ... 0 ? ? ?, U = ? ? ? ... ? ... ? ? 0 0 a11 0 ... 0

81

a12 . . . ... ... ... 0

a1n ... an?1,n 0

? ? ? ? ? ? ? ?

... ...

and

then A = L + D + U , which is a splitting of A. Now we can rewrite Ax = b as (D + D + U )x = b then Dx = ?(L + U )x + b, and we may reformulate the iterative methods as follows: Jacobi’s method where BJ is the Jacobi’s iteration matrix. Dxk+1 = ?(L + U )xk + b ? BJ = ?D ?1 (L + U ),

? ? ? ? ? 0 a22 0 . . . ? ?, D=? ? ? ? ... ... ... ... ? ? ? 0 . . . 0 ann

?

?

Example 23. Write the linear system in the matrix form x = BJ x + C ? 1 1 ? ? ? x1 = x2 + ? ? 2 2 ? 2x ? x = 1 ? 1 2 ? which in the matrix form is ? ?x + 2x = 1 ? ? 1 2 ? 1 1 ? ? x2 = x1 + 2 2 ? ? x1 x2 ? ? ?=? x1 x2 ? ? 0

1 2 1 2

0

?? ??

x1 x2

1 2

?

?+? ?

?

1 2 1 2

?

? , where ?

1 2 1 2

x=?

? , BJ = ?

?

0

1 2

0

? and C = ?

?

?.

82

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS

Example 24. Determine the same BJ by the formula BJ = ?D ?1 (L + U ), A=? ? 2 ?1 ?1 2 ? ?,L = ? ? 0 0 ? ?,U = ?

1 2

?1 0

?

0 ?1 0 0

?

?,D = ?

?

2 0 0 2

? ?

We can easily see that D ?1 = ? ? ? 0

1 2

0

?

and thus

?, ? 1 2 ? ?. 0

?? ? ? 1 0 ? 0 0 ?1 ? ? ? ?=? BJ = ?D ?1 (L + U ) = ? 2 ? ? 1 1 ?1 0 0 ? 2 2

Gauss-Seidel’s method As above we write Ax = b as (L + D + U )x = b but now we choose (L + D )x = ?Ux + b. Similar to the previous procedure we have (L + D )xk+1 = ?Uxk + b, and then BGS = ?(L + D )?1 U , where BGS is Gauss-Seidel’s iteration matrix. Relaxation Gauss-Seidel gives that (L + D )x = ?Ux + b, thus the iteration procedure is: Dxk+1 = Dxk ? [Lxk+1 + (D + U )xk ? b]. where ω is the Relaxation parameter, ω = 1 gives the Gauss-seidel iteration. Now we have that (ωL + D )xk+1 = [(1 ? ω )D ? ωU ]xk + ωb, thus the Relaxation iteration matrix is: Bω = (ωL + D )?1 [(1 ? ω )D ? ωU ].

3.3. EXERCISES

83

3.3

Exercises

Problem 28. Find the unique the LDU factorization for the matrix ? ? 1 1 ?3 ? ? ? ? A=? 0 1 1 ?. ? ? 3 ?1 1

Problem 27. Solve A4 x = b for ? ? ?1 2 ? A=? 2 ?3

Problem 26. Illustrate the LU factorization for the matrix ? ? 1 3 2 ? ? ? ? A = ? ?2 ?6 1 ? . ? ? 2 5 7 ? 144 ?233 ? ?

b=?

where c2 + s2 = 1

Problem 29. Show that every orthogonal 2 × 2 matrix is of the form ? ? ? ? c s c s ? ?, A1 = ? or A2 = ? ?s c s ?c

using 3 iterations of the following methods using a starting value of u0 = (0, 0)T . (a) Jacobi Method. (b) Gauss-Seidel Method. (c) Optimal SOR (you must compute the optimal value of ω = ω0 ?rst).

Problem 30. Solve the following system ? ?? ? ? ? 4 ?1 u 1 ? ?? 1 ? = ? ? ?1 4 u2 ?3

84

CHAPTER 3. LINEAR SYSTEM OF EQUATIONS

Chapter 4 Two-points BVPs

In this chapter we focus on ?nite element approximation procedure for the two point boundary value problems (BVPs) of Dirichlet, Neumann and mixed type. For each PDE we formulate a corresponding variational formulation(VF) and a minimization problem (MP) and prove that to solve the boundary value problem is equivalent to the (VF) which in turn is equivalent to solve a minimization problem (MP), i.e, (BV P ) ?? (V F ) ?? (MP ). We also prove a priori and a posteriori error estimates for BVPs.

4.1

A Dirichlet problem

Assume a horizontal elastic bar that occupies the interval I := [0, 1], is ?xed at the end-points. Let u(x) denote the displacement at a point x ∈ I , and a(x) be the modulus of elasticity and f (x) a load function, then one can show that u satis?es the following boundary value problem for the Poisson’s equation ? ′ ? ? a(x)u′ (x) = f (x), ? u(0) = u(1) = 0. 0 < x < 1,

(BV P )1

(4.1.1)

We assume that a(x) is piecewise continuous in (0, 1), bounded for 0 ≤ x ≤ 1 and a(x) > 0 for 0 ≤ x ≤ 1. 85

86

CHAPTER 4. TWO-POINTS BVPS

Let v (x) and its derivative v ′ (x), x ∈ I , be square integrable functions, that is: v, v ′ ∈ L2 (0, 1), and set

1 1 H0 = v (x) : 0

[v (x)2 + v ′ (x)2 ]dx < ∞,

v (0) = v (1) = 0 .

(4.1.2)

The variational formulation (VF) for (BVP)1 is obtained by multiplying the 1 equation by a function v (x) ∈ H0 (0, 1) and integrating over (0, 1):

1 1

?

[a(x)u (x)] v (x)dx =

0 0

′

′

f (x)v (x)dx.

(4.1.3)

By partial integration we get ? a(x)u (x)v (x)

′ 1 0 1 1

+

0

a(x)u (x)v (x)dx =

0

′

′

f (x)v (x)dx.

(4.1.4)

Now since v (0) = v (1) = 0 we have

1 1

a(x)u′ (x)v ′ (x)dx =

0 0

f (x)v (x)dx.

(4.1.5)

Thus the variational formulation for the problem (4.1.1) is as follows: 1 Find u(x) ∈ H0 such that

1 1

(VF)1

a(x)u′ (x)v ′ (x)dx =

0 0

f (x)v (x)dx,

1 ?v (x) ∈ H0 .

(4.1.6)

In other words we have that if u satis?es (BVP)1 (4.1.1), then u also satis?es the (VF)1 , (4.1.5) above. We write this as (BVP)1 =? (VF)1 . Now the question is whether the reverse implication is through, i.e. if which conditions can we deduce the reverse implication (VF)1 =? (BVP)1 ? It appears that this question has an a?rmative answer and the two problems are indeed equivalent. We prove this observation in the following theorem. Theorem 7. u satis?es (BVP)1 ?? u satis?es (VF)1 . Proof. We have already shown (BVP)1 =? (VF)1 . It remains to show that (VF)1 =? (BVP)1 . Integrating by parts on the left hand side in (4.1.5) and using v (0) = v (1) = 0 we come back to

1 1

?

[a(x)u′ (x)]′ v (x)dx =

0 0

f (x)v (x)dx,

1 ?v (x) ∈ H0

(4.1.7)

4.1.

A DIRICHLET PROBLEM

87

which can also be written as

1 0

? a(x)u′ (x)

′

? f (x) v (x)dx = 0,

′

1 ?v (x) ∈ H0 .

(4.1.8)

We claim that (4.1.8) implies ? a(x)u′ (x) ? f (x) ≡ 0,

′

?x ∈ (0, 1).

(4.1.9)

If our claim is not true!, then there exists at least one ξ ∈ (0, 1), such that ? a(ξ )u′ (ξ ) ? a(ξ )u′(ξ )

′

? f (ξ ) = 0 ,

(4.1.10)

where we may assume, without loss of generality, that ? f (ξ ) > 0 (or < 0). (4.1.11)

Thus, assuming that f ∈ C (0, 1) and a ∈ C 1 (0, 1), by continuity, ?δ > 0, such that in a δ -neighborhood of ξ , g (x) := ? a(x)u′ (x)

′

? f (x) > 0,

for

x ∈ (ξ ? δ, ξ + δ ).

(4.1.12)

Now we take v (x) in (4.1.8) as a hat function, v ? (x) > 0 on (ξ ? δ, ξ + δ ) and y 1 v ? (x) x

g (x)

0

ξ?δ

ξ

ξ+δ

1

Figure 4.1: The hat function v ? (x) over the interval (ξ ? δ, ξ + δ ).

1 1 we have v ? (x) ∈ H0 and 0

? a(x)u′ (x)

>0

′

? f (x) v ? (x) dx > 0, which

>0

contradicts (4.1.8), thus our claim is true and the proof is complete.

88

CHAPTER 4. TWO-POINTS BVPS

Corollary 1. (i) If both f (x) and a(x) are continuous and a(x) is di?erentiable, i.e. f ∈ C (0, 1) and a ∈ C 1 (0, 1), then (BVP) and (VF) have the same solution. (ii) If a(x) is discontinuous, then (BVP) is not always well-de?ned but (VF) has meaning. Therefore (VF) covers a larger set of data than (BVP).

4.2

Minimization problem

0 < x < 1,

we formulate a minimization problem (MP) as follows:

1 1 Problem 31. Find u ∈ H0 such that F (u) ≤ F (w ), ?w ∈ H0 , where F (w ) is the total energy of w (x) given by 1 1 1 ′ 2 F (w ) = a(w ) dx ? f wdx 2 0 0 Internal energy Load potential

For the problem (4.1.1), i.e. our (BVP)1 ? ′ ? ? a(x)u′ (x) = f (x), ? u(0) = u(1) = 0.

(4.2.1)

(4.2.2)

Theorem 8. The minimization problem above is equivalent to variational formulation (V F )1 , (MP ) ?? (V F ) i.e.

1 1 1 F (u) ≤ F (w ), ?w ∈ H0 ??

au′ v ′ dx =

0 0

f vdx,

1 ?v ∈ H0 .

(4.2.3)

1 1 Proof. (?=) For w ∈ H0 we let v = w ? u, then v ∈ H0 and

1 F (w ) = F (u + v ) = 2 1 = 2

1

a (u + v )

0 1 0 ′ ′

′

2

1

dx ?

1 0

f (u + v )dx =

0 1 ′ 2

1 2au v dx + 2

(i) 1 1

1 a(u ) dx + 2

(ii)

a(v ′ )2 dx

0

?

0

f udx ?

(iii)

f vdx.

0 (iv)

4.3.

A MIXED BOUNDARY VALUE PROBLEM

1 1

89 f vdx. Further by

0

but (i) + (iv ) = 0, since by (VF)1

au′ v ′ dx =

0

de?nition of the functional F we have (ii) + (iii) = F (u). Thus F (w ) = F (u ) + 1 2

1

a(x)(v ′ (x))2 dx,

0

(4.2.4)

and since a(x) > 0 we have F (w ) > F (u). (=?) Let now F (u) ≤ F (w ) and set g (ε, w ) = F (u + εv ), then g has a = 0. We have minimum at ε = 0. In other words g ′ (ε, w )

ε=0 2

g (ε, w ) = F (u + εv ) = 1 = 2

1 0 ′ 2

1 2

1

1

a (u + εv )′

0 2 ′ 2

dx ?

f (u + εv )dx =

0 1 1

{a(u ) + aε (v ) + 2aεu v }dx ?

′ ′

0

f udx ? ε

f vdx.

0

′ Now we compute the derivative gε (ε, w ).

1 ′ gε (ε, w ) = {2aε(v ′ )2 + 2au′ v ′ }dx ? 2

′ and gε (ε=0)

1

f vdx

0

(4.2.5)

= 0, yields

1 0 1

au′ v ′ dx ?

f vdx = 0,

0

(4.2.6)

which is the variational formulation. Thus we conclude that F (u) ≤ F (w ) =? (VF)1 and the proof is complete. We summarize the two theorems in short as Corollary 2. (BV P )1 ?? (V F )1 ?? (MP ).

4.3

A mixed Boundary Value Problem

Obviously changing the boundary conditions would require changes in the variational formulation. This can be, e.g. seen in formulating the (VF)

90

CHAPTER 4. TWO-POINTS BVPS

corresponding to the following mixed boundary value problem ? ′ ? ? a(x)u′ (x) = f (x), 0 < x < 1 (BVP)2 ? u(0) = 0, a(1)u′ (1) = g1 = 0.

1 1

(4.3.1)

We multiply the equation by a suitable function v (x) with v (0) = 0 and integrate over the interval (0, 1) to obtain ? [a(x)u (x)] v (x)dx =

0 0 ′ ′

f (x)v (x)dx.

(4.3.2)

By partial integration we get, as before, that

1 1

?[a(x)u′ (x)v (x)]1 0 +

a(x)u′ (x)v ′ (x)dx =

0 0

f (x)v (x)dx

(4.3.3)

Using the boundary data v (0) = 0 and a(1)u′(1)v (1) = g1 v (1) we get

1 1

a(x)u′ (x)v ′ (x)dx =

0 0

f (x)v (x)dx + g1 v (1),

? 1, ?v ∈ H 0

(4.3.4)

where

1 ?0 H = {v (x) : 1 0

[v (x)2 + v ′ (x)2 ]dx < ∞, such that v (0) = 0}.

(4.3.5)

? 1 such that Then (4.3.4) yields the variational formulation: Find u ∈ H 0

1 1

(VF)2

a(x)u′ (x)v ′ (x)dx =

0 0

f (x)v (x)dx + g1 v (1),

?1 ?v ∈ H 0

Now we want to show that Theorem 9. (BVP)2 ?? (VF)2 Proof. (=?) This part is trivial and already proved along the above lines. (?=) To prove that a solution of the variational problem (VF)2 is also a solution of the two-point boundary value problem (BVP)2 we have to show (i) the solution satis?es the di?erential equation

4.3.

A MIXED BOUNDARY VALUE PROBLEM

91

(ii) the solution satis?es the boundary conditions We start with (V F )2 and perform a reversed order partial integration to get

1 1

a(x)u (x)v (x)dx = [a(x)u

0

′

′

′

(x)v (x)]1 0

?

[a(x)u′ (x)]′ v (x) dx.

0

(4.3.6)

Since v (0) = 0, we get

1 0 1

a(x)u′ (x)v ′ (x)dx = a(1)u′ (1)v (1) ?

[a(x)u′ (x)]′ vdx

0

(4.3.7)

Thus the variational formulation (VF)2 can be written as

1 1

?

[a(x)u′ (x)]′ vdx + a(1)u′(1)v (1) =

0 0

f (x)v (x)dx + g1 v (1).

(4.3.8)

? 1 (0, 1), including a test funcThe equation (4.3.8) is valid for every v (x) ∈ H 0 tion v (x) with v (0) = v (1) = 0 as in the Dirichlet problem: ?(au′ )′ = 1 ? 1 (0, 1). Consef, u(0) = u(1) = 0. This is simply because H0 (0, 1) ? H 0 quently choosing v (1) = 0 (4.3.8) is reduced to

1 1

?

[a(x)u (x)] vdx =

0 0

′

′

f (x)v (x)dx,

1 ?v (x) ∈ H0

(4.3.9)

Now as in the case of the Dirichlet problem (4.3.9) gives the di?erential equation in (4.3.1) and hence claim (i) is through. On the other hand (4.3.9) is just the equation in (4.3.1) multiplied by a test function v and integrated over (0, 1), so (4.3.9) is equally valid for v ∈ 1 ?0 H (0, 1). Now inserting (4.3.9) in (4.3.8) we also get g1 v (1) = a(1)u′ (1)v (1), which choosing v (1) = 0, e.g. v (1) = 1, gives that g1 = a(1)u′ (1) and the proof is complete. Remark 10. i) The Dirichlet boundary conditions is called the essential boundary condition and is strongly imposed in the test function space: Enforced explicitly to the trial and test functions in (VF). ii) The Neumann and Robin Boundary conditions are called the natural boundary conditions and are automatically satis?ed in (VF), therefore are weakly imposed.

92

CHAPTER 4. TWO-POINTS BVPS

4.4

The ?nite element method. (FEM)

Let Th = {0 = x0 < x1 < . . . < xM < xM +1 = 1} be a partition of 0 ≤ x ≤ 1 into subintervals Ik = [xk?1 , xk ] and hk = xk ? xk?1 De?ne the piecewise x0 = 0 x1 x2 xk?1 hk xk xM xM +1 = 1 x

constant function h(x) = xk ? xk?1 = hk for x ∈ Ik . Let now C I, P1 (Ik ) denote the set of all continuous piecewise linear functions on Th (continuous in whole I , linear on each subinterval Ik ), and de?ne Vh Note that Vh

(0) (0)

= {v : v ∈ C I, P1 (Ik ) ,

v (0) = v (1) = 0}

(4.4.1)

is a subspace of

1

1 H0 = {v (x) :

0

[v (x)2 + v ′ (x)2 ]dx < ∞,

and

v (0) = v (1) = 0}. (4.4.2)

A cG(1) (continuous Galerkin of degree 1) ?nite element formulation for our (0) Dirichlet boundary value problem (BVP) is given by: Find uh ∈ Vh such that

1 1

(FEM)

0

a(x)u′h (x)v ′ (x)dx =

f (x)v (x)dx,

0

?v ∈ Vh .

(0)

(4.4.3)

Now the purpose is to make estimate of error arising in approximating the (0) solution for (BV P ) by the functions in Vh . To this end we need to introduce some measuring environment for the error. Recall the de?nition of Lp -norms

1

Lp -norm L∞ -norm Weighted L2 -norm Energy-norm Note that

v v v v v

Lp L∞

=

0

|v (x)|p dx

1/p

,

1≤p<∞

= sup |v (x)|

x∈[0,1] 1 0

a

=

1

a(x)|v (x)|2 dx a(x)|v ′ (x)|2 dx

1/2

,

a(x) > 0

E E

=

0

1/2

= v′ a.

4.5. ERROR ESTIMATES IN THE ENERGY NORM

93

v E describes the elastic energy for an elastic string modeled for our Dirichlet boundary value problem (BVP).

4.5

Error estimates in the energy norm

There are two types of error estimates: An a priori error estimate depends on the exact solution u(x) and not on the approximate solution uh (x). In such estimates the error analysis are performed theoretically and before computations. An a posteriori error estimate where the error depends on the residual,i.e, the di?erence between the left and right hand side in the equation when the exact solution u(x) is replaced by the approximate solution uh (x). A posteriori error estimates can be derived after that the approximate solution is computed. Below ?rst we shall prove a general theorem which shows that the ?nite element solution is the best approximate solution for either of our Dirichlet problem in the energy norm. Theorem 10. Let u(x) be a solution of the Dirichlet boundary value problem ? ′ ? ? a(x)u′ (x) = f (x), ? u(0) = 0 u(1) = 0. 0<x<1

BVP

(4.5.1)

and uh (x) its ?nite element element approximation given by (4.4.3). Then we have u ? uh

E

≤ u?v

E , ?v (x)

∈ Vh .

(0)

(0)

(4.5.2)

This means that the ?nite element solution uh ∈ Vh (0) tion of the solution u by functions in Vh .

is the best approxima-

Proof. Recall the variational formulation associated to the problem (4.4.1):

1 1

(VF)

0

a(x)u′ (x)v ′ (x)dx =

0

f (x)v (x)dx,

1 ?v ∈ H0 .

(4.5.3)

94

(0)

CHAPTER 4. TWO-POINTS BVPS

We take an arbitrary v ∈ Vh , then by the de?nition of the energy norm

1

u ? uh

2 E

=

0 1

a(x)(u′ ? u′h )2 (x)dx a(x) u′ (x) ? u′h (x)

1

=

0

u′(x) ?v ′ (x) + v ′ (x) ?u′h (x) dx

=0

=

0

a(x) u′ (x) ? u′h (x)

1 0

u′(x) ? v ′ (x) dx v ′ (x) ? u′h (x) dx (4.5.4)

+

a(x) u′(x) ? u′h (x)

(0)

Now since v ? uh ∈ Vn

1 0

1 ? H0 , we have by the variational formulation (4.5.3) 1

a(x)u′ (x) v ′ (x) ? u′h (x) dx =

0

f v (x) ? uh (x) ,

(4.5.5)

with its ?nite element counterpart, see (4.4.3),

1 0 1

a(x)u′h (x) v ′ (x) ? u′h (x) dx =

0

f v (x) ? uh (x) .

(4.5.6)

Subtracting these two relations the last line of the estimate (4.5.4) above vanishes, so we end up with

1

u ? uh

2 E

=

0 1

a(x)[u′ (x) ? u′h (x)][u′ (x) ? v ′ (x)]dx a(x) 2 [u′ (x) ? u′h (x)]a(x) 2 [u′ (x) ? v ′ (x)]dx

1 0

1 1

=

0

≤

a(x)[u (x) ?

E

′

u′h (x)]2 dx

E,

1 2

1 0

a(x)[u (x) ? v (x)] dx

′

′

2

1 2

= u ? uh

· u?v

(4.5.7)

where, in the last estimate, we used Cauchy-Schwartz inequality. Thus u ? uh and the proof is complete.

E

≤ u?v

E,

(4.5.8)

4.5. ERROR ESTIMATES IN THE ENERGY NORM

(0)

95

Next step is to show that there exists a function v (x) ∈ Vh such that u ? v E is not too large. The function that we shall study is v (x) = πh u(x): the piecewise linear interpolant of u(x), introduced in chapter 2. Recall the interpolation error estimate in Lp -norms: Theorem 11. (i) Let 0 = x0 < x1 < x2 < . . . < xM < xM +1 = 1 be a partition of [0, 1] and h = (xj +1 ? xj ), j = 0, 1, . . . , M (ii) Let πh v (x) be the piecewise linear interpolant of v (x). Then there is an interpolation constant ci such that πh v ? v Lp ≤ ci h2 v ′′ Lp 1 ≤ p ≤ ∞ (πh v )′ ? v ′ Lp ≤ ci hv ′′ Lp πh v ? v Lp ≤ ci hv ′ Lp . (4.5.9) (4.5.10) (4.5.11)

Theorem 12 (An apriori error estimate). Let u and uh be the solutions of the Dirichlet problem (BVP) and the ?nite element problem (FEM), respectively. Then there exists an interpolation constant Ci , depending only on a(x), such that u ? uh E ≤ Ci hu′′ a . (4.5.12) Proof. According to our general above we have u ? uh

(0) E

≤ u?v

E,

?v ∈ Vh .

(0)

(4.5.13)

Now since πh u(x) ∈ Vh , we may take v = πh u(x) in (4.5.13) and use, e.g. the second estimate in the interpolation theorem to get u ? uh

E

≤ u ? πh u ≤ Ci hu′′

a

E

= u′ ? (πh u)′

1 0

a 1/2

= Ci

a(x)h2 (x)u′′ (x)2 dx

,

(4.5.14)

which is the desired result. Remark 11. Now if the objective is to divide (0,1) into a ?xed, ?nite, number of subintervals, then one can use the proof of theorem 8.3: to obtain an optimal (a best possible) partition of (0,1); in the sense that: whenever a(x)u′′ (x) gets large we compensate by making h(x) smaller. This, however, requires that the exact solution u(x) is known. Now we want to study a posteriori error analysis, where instead of the unknown value of u(x), we use the known computed values of the approximate solution.

96

CHAPTER 4. TWO-POINTS BVPS

Theorem 13 (a posteriori error estimate). There is an interpolation constant ci depending only on a(x) such that the error in ?nite element approximation of the Dirichlet boundary value problem (BVP) (4.5.1) satis?es

1

e(x)

E

≤ ci

0

1 2 h (x)R2 [uh (x)]dx a(x)

1 2

,

(4.5.15)

1 where e(x) = u(x) ? uh (x), note that e ∈ H0 . 1 1

Proof. By the de?nition of the energy norm we have e(x)

2 E

=

0 1

a(x)[e′ (x)]2 dx =

0

a(x)[u′ (x) ? u′h (x)]e′ (x)dx

1 0

(4.5.16)

=

0

a(x)u (x)e (x)dx ?

1

′

′

a(x)u′h (x)e′ (x)dx

Since e ∈

1 H0

the variational formulation (VF) gives that

1

a(x)u (x)e (x)dx =

0 0

′

′

f (x)e(x)dx.

(4.5.17)

Thus we have

1 1

e(x)

2 E

=

0

f (x)e(x)dx ?

0

a(x)u′h (x)e′ (x)dx.

(4.5.18)

Adding and subtracting the interpolant πh e(x) and πh e′ (x) to e and e′ in the integrands above yields

1 1

e(x)

2 E

=

0

f (x)[e(x) ? πh e(x)]dx +

1

f (x)πh e(x)dx

0 (i) 1

?

0

a(x)u′h (x)[e′ (x) ? πh e′ (x)]dx ?

0

a(x)u′h (x)πh e′ (x)dx .

(ii)

Since uh (x) is a solution of the (FEM) (4.4.3) and πh e(x) ∈ Vh ?(ii) + (i) = 0. Hence

1 1

(0)

we have

e(x) =

2 E 1 0

=

0

f (x)[e(x) ? πh e(x)]dx ?

M xk

0

a(x)u′h (x)[e′ (x) ? πh e′ (x)]dx a(x)u′h (x)[e′ (x) ? (πh e′ (x)]dx.

f (x)[e(x) ? πh e(x)]dx ?

k =1

x k ?1

4.5. ERROR ESTIMATES IN THE ENERGY NORM

97

Now, for the integrals in the sum above, we integrate by parts over each subinterval (xk?1 , xk ):

xk

?

x k ?1

a(x)u′h (x) (e′ (x) ? πh e′ (x)) dx = [P.I.] =

g (x) F ′ (x) xk x k ?1 xk

= ? a(x)u′h (x) (e(x) ? πh e(x))

g (x) F (x)

+

x k ?1

(a(x)u′h (x))′ (e(x) ? πh e(x)) dx

g ′ (x) F (x)

Since e(xk ) = πh e(xk ), k = 0, 1 . . . , M , where xk :s are the interpolation nodes we have F (xk ) = F (xk?1 ) = 0, and thus

xk

?

x k ?1

a(x)u′h (x)(e′ (x)

? πh e (x))dx =

′

xk

x k ?1

a(x)u′h (x) (e(x) ? πh e(x))dx.

′

Hence summing over k , we get

1 1

?

0

a(x)u′h (x)[e′ (x) ? πh e′ (x)]dx =

1

0

[a(x)u′h (x)]′ (e(x) ? πh e(x))dx,

1

and therefore e(x)

2 E

=

0 1

f (x)[e(x) ? πh e(x)]dx +

0

[a(x)u′h (x)]′ (e(x) ? πh e(x))dx

=

0

{f (x) + [a(x)u′h (x)]′ }(e(x) ? πh e(x))dx,

Let now R(uh (x)) = f (x) + (a(x)u′h (x))′ , i.e. R(uh (x)) is the residual error, which is a well-de?ned function except in the set {xk }, since (a(xk )u′x (xk ))′ are not de?ned. Thus we can get the following estimate

1

e(x)

2 E

=

0 1

R(uh (x))(e(x) ? πh e(x))dx = 1

1

=

0

≤

0

a(x) 1 2 h (x)R2 (uh (x))dx a(x)

h(x)R(uh (x)) ·

a(x)

1 2

e(x) ? πh e(x) dx h(x)

1 0

e(x) ? πh e(x) a(x) h(x)

2

dx

1 2

,

where we have used Cauchy Schwarz inequality. Now recalling the de?nition of the weighted L2 -norm we have, e(x) ? πh e(x) h(x)

1 a

=

0

a(x)

e(x) ? πh e(x) h(x)

2

dx

1 2

.

(4.5.19)

98

CHAPTER 4. TWO-POINTS BVPS

To estimate (4.5.19) we use the third interpolation estimate for e(x) in a subinterval and get e(x) ? πh e(x) h(x) Thus

1

a

≤ ci e′ (x)

a

= ci e(x)

E.

(4.5.20)

e(x)

2 E

≤

0

1 2 h (x)R2 (uh (x))dx a(x)

1 2

· ci e(x)

E,

(4.5.21)

and the proof is complete. Adaptivity Below we brie?y outline the adaptivity procedure based on the a posteriori error estimate which uses the approximate solution and which can be used for mesh-re?nements. Loosely speaking this predicts local mesh re?nement, i.e. indicates changing the length of the interval h(x) in the regions (subintervals) which is necessary. More concretely the idea is as follows: Assume that one seeks an error bound less that a given error tolerance TOL: e(x)

E

≤ TOL.

(4.5.22)

Then one may use the following steps as a mesh re?nement strategy: (i) Make an initial partition of the interval (ii) Compute the corresponding FEM solution uh (x) and residual R(uh (x)). (iii) If e(x)

E

> TOL re?ne the mesh in the places for which

1 R2 (uh (x)) a(x)

is large and perform the steps (ii) and (iii) again.

4.6. EXERCISES

99

4.6

Exercises

?u′′ = f, 0 < x < 1; u(0) = u(1) = 0. (4.6.1)

Problem 32. Consider the two-point boundary value problem

Let V = {v : v + v ′ < ∞,

a. Use V to derive a variational formulation of (4.6.1). b. Discuss why V is valid as a vector space of test functions. c. Classify whether the following functions are admissible test functions or not: sin πx, x2 , x ln x, ex ? 1, x(1 ? x). Problem 33. Assume that u(0) = u(1) = 0, and that u satis?es

1 1

v (0) = v (1) = 0}.

u′ v ′ dx =

0 0

f v dx, v (0) = v (1) = 0}.

1

for all v ∈ V = {v : v + v ′ < ∞, 1 2

1 0

a. Show that u minimizes the functional F (v ) =

(v ′ )2 dx ?

f v dx.

0

(4.6.2)

Hint: F (v ) = F (u + w ) = F (u) + . . . ≥ F (u). ?u′′ = f, ?u′′ = 1, 0 < x < 1;

b. Prove that the above minimization problem is equivalent to u(0) = u(1) = 0.

Problem 34. Consider the two-point boundary value problem 0 < x < 1; u(0) = u(1) = 0. (4.6.3)

j , j = 0, 1, . . . , 4, denote a partition of the interval 0 < x < 1 Let Th : xj = 4 into four subintervals of equal length h = 1/4 and let Vh be the corresponding space of continuous piecewise linear functions vanishing at x = 0 and x = 1.

b. Prove that U ∈ Vh is unique.

a. Compute a ?nite element approximation U ∈ Vh to (4.6.3).

100

CHAPTER 4. TWO-POINTS BVPS

Problem 35. Consider once again the two-point boundary value problem ?u′′ = f, 0 < x < 1; u(0) = u(1) = 0.

a. Prove that the ?nite element approximation U ∈ Vh to u satis?es (u ? U ) ′ ≤ (u ? v ) ′ , b. Use this result to deduce that (u ? πh u)′ ≤ C hu′′ , where C is a constant and πh u a piecewise linear interpolant to u. Problem 36. Consider the two-point boundary value problem ?(au′ )′ = f, u(0) = 0, 0 < x < 1, a(1)u′ (1) = g1 , (4.6.5) (4.6.4) for all v ∈ Vh .

where a > 0 is a positive function and g1 is a constant. a. Derive the variational formulation of (4.6.5). b. Discuss how the boundary conditions are implemented. Problem 37. Consider the two-point boundary value problem ?u′′ = 0, 0 < x < 1; u(0) = 0, u′ (1) = 7. (4.6.6)

1 Divide the interval 0 ≤ x ≤ 1 into two subintervals of length h = 2 and let Vh be the corresponding space of continuous piecewise linear functions vanishing at x = 0.

a. Formulate a ?nite element method for (4.6.6). b. Calculate by hand the ?nite element approximation U ∈ Vh to (4.6.6). Study how the boundary condition at x = 1 is approximated. Problem 38. Consider the two-point boundary value problem ?u′′ = 0, 0 < x < 1; u′ (0) = 5, u(1) = 0. (4.6.7)

4.6. EXERCISES

101

Let Th : xj = jh, j = 0, 1, . . . , N, h = 1/N be a uniform partition of the interval 0 < x < 1 into N subintervals and let Vh be the corresponding space of continuous piecewise linear functions. a. Use Vh to formulate a ?nite element method for (4.6.7). b. Compute the ?nite element approximation U ∈ Vh assuming N = 3. Problem 39. Consider the problem of ?nding a solution approximation to ?u′′ = 1, 0 < x < 1; u′ (0) = u′ (1) = 0. (4.6.8)

Let Th be a partition of the interval 0 < x < 1 into two subintervals of equal length h = 1 and let Vh be the corresponding space of continuous piecewise 2 linear functions. a. Find the exact solution to (4.6.8) by integrating twice. b.Compute a ?nite element approximation U ∈ Vh to u if possible. Problem 40. Consider the two-point boundary value problem ?((1 + x)u′ )′ = 0, 0 < x < 1; u(0) = 0, u′ (1) = 1. (4.6.9)

1 and Divide the interval 0 < x < 1 into 3 subintervals of equal length h = 3 let Vh be the corresponding space of continuous piecewise linear functions vanishing at x = 0.

a. Use Vh to formulate a ?nite element method for (4.6.9). b. Verify that the sti?ness matrix A and the load vector b are given by ? ? ? ? 0 16 ?9 0 ? ? ? 1? ? ? ? ? b = ? 0 ?. A = ? ?9 20 ?11 ? , 2? ? ? ? 1 0 ?11 11 c. Show that A is symmetric tridiagonal, and positive de?nite. d. Derive a simple way to compute the energy norm U

1 2 E,

de?ned by

U

2 E

=

0

(1 + x)U ′ (x)2 dx,

where U ∈ Vh is the ?nite element solution approximation.

102

CHAPTER 4. TWO-POINTS BVPS

Problem 41. Consider the two-point boundary value problem ?u′′ = 0, 0 < x < 1; u(0) = 0, u′ (1) = k (u(1) ? 1). (4.6.10)

2 1 and x1 = 3 be a partition Let Th : 0 = x0 < x1 < x2 < x3 = 1, where x1 = 3 of the interval 0 ≤ x ≤ 1 and let Vh be the corresponding space of continuous piecewise linear functions, which vanish at x = 0.

a. Compute a solution approximation U ∈ Vh to (4.6.10) assuming k = 1. Problem 42. Consider the ?nite element method applied to ?u′′ = 0, 0 < x < 1; u(0) = α, u′ (1) = β,

b. Discuss how the parameter k in?uence the boundary condition at x = 1.

where α and β are given constants. Assume that the interval 0 ≤ x ≤ 1 is divided into three subintervals of equal length h = 1/3 and that {?j }3 0 is a nodal basis of Vh , the corresponding space of continuous piecewise linear functions. a. Verify that the ansatz U (x) = α?0 (x) + ξ1 ?1 (x) + ξ2 ?2 (x) + ξ3 ?3 (x), yields the following system of equations ? ? ? ? 0 ?1 2 ?1 0 ? ? ?? ξ ? ? 1 ? ?=? 0 ?1 2 ?1 ? ? ? ? 0 ?? ? ξ2 ? ? ? β 0 0 ?1 1 ? ξ3 ? ? α ? ?

1? ? ? h?

? ? ?. ?

(4.6.11)

b. If α = 2 and β = 3 sgow that (4.6.11) can be reduced to ? ? ? ?? ? ?1 ?2h ξ 2 ?1 0 ? ?? 1 ? ? 1? ? ? ? ?? ? ?. ? ?1 0 2 ?1 ? ? ξ2 ? = ? ? ? ? ?? h? 3 ξ3 0 ?1 1 c. Solve the above system of equations to ?nd U (x).

4.6. EXERCISES Problem 43. Compute a ?nite element solution approximation to ?u′′ + u = 1; 0 ≤ x ≤ 1, u(0) = u(1) = 0,

103

(4.6.12)

using the continuous piecewise linear ansatz U = ξ1 ?1 (x) + ξ2 ?2 (x) where ? ? 1 1 ? ? ? ? 3x, 0<x< 3 0, 0<x< 3 ? ? ? ? 1 1 2 2 ?1 (x) = ?2 (x) = 2 ? 3x, 3 3x ? 1, 3 <x< 3 , <x< 3 . ? ? ? ? ? ? 2 ? 0, ? 3 ? 3x, 2 < x < 1 <x<1 3 3 Problem 44. Consider the following eigenvalue problem ?au′′ + bu = 0; 0 ≤ x ≤ 1, u(0) = u′ (1) = 0, (4.6.13)

where a, b > 0 are constants. Let Th : 0 = x0 < x1 < . . . < xN = 1, be a non-uniform partition of the interval 0 ≤ x ≤ 1 into N intervals of length hi = xi ? xi?1 , i = 1, 2, . . . , N and let Vh be the corresponding space of continuous piecewise linear functions. Compute the sti?ness and mass matrices.

104

CHAPTER 4. TWO-POINTS BVPS

Chapter 5 Scalar Initial Value Problem

Consider the following ordinary di?erential equation (ODE) ? (DE) ? u ˙ (t) + a(t)u(t) = f (t), 0 < t ≤ T (IV) ? u(0) = u

0

(5.0.1)

du . Here a(t) is a bounded function. dt For a(t) ≥ 0 (5.0.1) is called a parabolic problem, while a(t) > 0 yields a dissipative problem. Below ?rst we give a few analytic aspects where f (t) is the source term and u ˙ (t) =

5.1

Fundamental solution and stability

t

Theorem 14 (Fundamental solution). The solution for the ODE (5.0.1) is given by u(t) = u0 · e?A(t) + where A(t) =

t 0

e?(A(t)?A(s)) f (s)ds,

0

(5.1.1)

a(s)ds is the integrating factor. d [u(t)eA(t) ] = eA(t) f (t), dt

t

Proof. Multiplying the (DE) by the integrating factor eA(t) we get ˙ (t)eA(t) u(t) = eA(t) f (t), u ˙ (t)eA(t) + A i.e.

˙ (t). Integrating over (0, t) yields where we used that a(t) = A

t 0

d [u(s)eA(s) ]ds = ds

t 0

eA(s) f (s)ds ?? u(t)eA(t) ?u(0)eA(0) = 105

eA(s) f (s)ds.

0

106

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

Now since A(0) = 0 and u(0) = u0 we get the desired result

t

u(t) = u0 · e?A(t) +

e?(A(t)?A(s)) f (s)ds.

0

(5.1.2)

Theorem 15 (Stability estimates). Using the fundamental solution we can derive the following stability estimates: 1 (i) If a(t) ≥ α > 0, then |u(t)| ≤ e?αt |u0 | + (1 ? e?αt ) max |f (s)| 0≤s≤t α (ii) If a(t) ≥ 0 (i.e. α = 0 the parabolic case), then

t

|u(t)| ≤ |u0 | +

0

|f (s)|ds or |u(t)| ≤ |u0 | + f

L1 t

(5.1.3)

Proof. (i) For a(t) ≥ 0, ?t > 0, we have that A(t) =

t t

a(s)ds is non0

decreasing and A(t) ? A(s) ≥ 0, ?t > s. For a(t) ≥ α > 0 we have A(t) =

0

a(s)ds ≥

0

α · ds = αt. Further

t

A(t) ? A(s) =

s

a(r ) dr ≥ α(t ? s).

t

(5.1.4)

Thus e?A(t) ≤ e?αt and e?(A(t)?A(s)) ≤ e?α(t?s) . Hence using (5.1.2) we get u(t) ≤ u0 · e

?αt

+

0

e?α(t?s) max |f (s)|ds,

0≤s≤t

(5.1.5)

which after integration gives that |u(t)| ≤ e?αt |u0 | + max |f (s)|

0≤s≤t

1 ?α(t?s) e α

s =t s=0

1 |u(t)| ≤ e?αt |u0 | + (1 ? e?αt ) max |f (s)|. 0≤s≤t α

t

complete.

(ii) Let α = 0 in (5.1.5) then |u(t)| ≤ |u0| +

0

|f (s)|ds, and the proof is

Remark 12. Recall that we refer to the set of functions where we seek the approximate solution as the trial space and the space of functions used for the orthogonality condition, as the test space.

5.2. GALERKIN FINITE ELEMENT METHODS (FEM) FOR IVP

107

5.2

Galerkin ?nite element methods (FEM) for IVP

To start, for discretization in time we shall introduce some general class of piecewise polynomial test and trial functions. However, in most of our studies in this notes we shall restrict ourselves to two simple cases: ? cG(1), continuous Galerkin of degree 1: In this case the trial functions are piecewise linear and continuous while the test functions are piecewise constant and discontinuous, i.e. unlike the cG(1) for BVP, here the trial and test functions are indi?erent spaces. ? dG(0), Discontinuous Galerkin of degree 0: Here both the trial and test functions are piecewise constant and discontinuous, i.e. like the cG(1) for BVP they are in the same space of functions, however, they are of one lower degree (piecewise constant) and discontinuous. Generally we have ? gG(q), Global Galerkin of degree q : Formulated for our initial value problem (5.0.1) as follows: Find U ∈ P q (0, T ) with U (0) = u0 such that

T 0

˙ + aU )vdt = (U

0

T

f v dt,

?v ∈ P q (0, T ), with v (0) = 0,

(5.2.1)

where v := {t, t2 , . . . , tq } := span[t, t2 , . . . , tq ]. ? cG(q), Continuous Galerkin of degree q : Find U ∈ P q (0, T ) with U (0) = u0 such that

T 0

˙ + aU )vdt = (U

0

T

f vdt,

?v ∈ P q?1 (0, T ),

(5.2.2)

where now v := {1, t, t2 , . . . , tq?1 }. Note the di?erence between the two test function spaces above. Example 25. Consider cG(q) with q = 1 then tq?1 = t0 = 1 and v ≡ 1, thus

T 0

˙ + aU )vdt = (U

0

T

˙ + aU )dt = U (T ) ? U (0) + (U

T

aU (t)dt

0

(5.2.3)

But U (t) is a linear function through U (0) and the unknown quantity U (T ), thus t T ?t U (t) = U (T ) + U (0) , (5.2.4) T T

108

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

inserting U (t) in (5.2.3) we get

T

U (T ) ? U (0) +

0

t T ?t dt = a U (T ) + U (0) T T

T

f dt.

0

(5.2.5)

which gives us U (T ) and consequently, through (5.2.4) and a given U (0), U (t). Using this idea we can formulate: ? The cG(1) Algorithm for the partition Tk of [0, T ] to subintervals Ik = (tk?1 , tk ]. (1) Given U (0) = U0 , apply (5.2.5) to (0, t1 ] and compute U (t1 ). Then using (5.2.4) one gets automatically U (t), ?t ∈ [0, t1 ]. (2) Assume that U is computed in all the successive intervals (tk?1 , tk ], k = 0, 1, n ? 1. (3) Compute U (t) for t ∈ (tn?1 , tn ]. This is done through applying (5.2.5) to the interval (tn?1 , tn ], instead of (0, T ]: i.e. with Un := U (tn ) and Un?1 := U (tn?1 ),

tn

Un ? Un?1 +

a

tn?1

t ? tn?1 tn ? t Un + Un?1 dt = tn ? tn?1 tn ? tn?1

tn

f dt.

tn?1

Now since Un?1 is known we can calculate Un and then U (t), t ∈ (tn?1 , tn ] is determined by the nth-version of the relation formula (5.2.4): U (t) = Un t tn ? t + Un?1 . tn tn

Global forms ?Continuous Galerkin cG(q): Find U (t) ∈ Vk , such that U (0) = U0 and tn tn (q ?1) ˙ f wdt, ?w ∈ Wk , (5.2.6) (U + aU )wdt =

0 0 (q )

Vk Wk

(q )

= {v : v continuous piecewise polynomials of degree q on Tk },

(q ?1)

= {w : w discontinuous piecewise polynomials of degree q ?1 on Tk }.

5.2. GALERKIN FINITE ELEMENT METHODS (FEM) FOR IVP ?Discontinuous Galerkin dG(q): Find U (t) ∈ P q (0, T ) such that

T 0

109

˙ + aU )vdt + a(U (0) ? u(0))v (0) = (U

T

f vdt,

0

?v ∈ P q (0, T ). (5.2.7)

This approach gives up the requirement that U (t) satis?es the initial condition. Instead, the initial condition is represented by U (0) ? u(0) = 0. In the sequel we shall use the following notation: + ? v (tn ± s) and [vn ] = vn ? vn is the jump in v (t) at time t. Let v ±n = lim +

s→0 +? vn

? [vn ]

? kn tn?1

? ? vn

tn

tn+1

t

Figure 5.1: The jump [vn ] and the right and left limits v ± Then dG(q) reads as follows: For n = 1, . . . , N ?nd U (t) ∈ P q (tn?1 , tn ) such that

tn tn?1 + + ˙ + aU )vdt + Un (U ?1 vn?1 = tn tn?1 q ? + f vdt + Un ?1 vn?1 , ?v ∈ P (tn?1 , tn ).

(5.2.8) Let q = 0, then v ≡ 1 is the only base function and we have U (t) = Un = + ? ˙ Un ?1 = Un on In = (tn?1 , tn ] and U ≡ 0. Thus for q = 0 (5.2.8) gives the dG(0) formulation: For n = 1, . . . , N ?nd piecewise constants Un such that

tn tn

aUn dt + Un =

tn?1 tn?1

f dt + Un?1 .

(5.2.9)

Finally summing over n in (5.2.8), we get the global dG(q) formulation: Find (q ) ? U (t) ∈ Wk , with U0 = u0 such that

N tn tn?1 N

˙ + aU )wdt + (U

n=1

n=1

+ [Un?1 ]wn ?1 =

tN

f wdt,

0

?w ∈ Wk . (5.2.10)

(q )

110

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

5.3

An a posteriori error estimate for cG(1)

u ˙ (t) + a(t)u(t) = f (t), ?t ∈ (0, T ), u(0) = u0 . (5.3.1)

?The continuous problem Recall the initial value problem

Let us rewrite (5.3.1) in a general variational form

T T

(u ˙ + au)vdt =

0 0

f vdt,

for all test functions v . Integrating by parts we get the equivalent equation

T T

u(T )v (T ) ? u(0)v (0) +

0

u(t) ? v ˙ (t) + av (t) dt =

f vdt.

0

(5.3.2)

If we now choose v to be the solution of the dual problem: ?v ˙ + av = 0, then (5.3.2) is simpli?ed to

T

in (0, T ),

(5.3.3)

u(T )v (T ) = u(0)v (0) +

0

f vdt,

?v (t) ∈ P q (0, T ).

(5.3.4)

In other words choosing v to be the solution of the dual problem (5.3.3) we may get the ?nal value u(T ) of the solution directly coupled to the initial value u(0) and the data f . This type of representation will be crucial in, e.g. a posteriori error analysis as in the proof of the next theorem. The Dual problem for (5.3.1) is formulated as follows: Find ?(t) such that ? ? ?? ˙ (t) + a(t)?(t) = 0, tN > t ≥ 0 (5.3.5) ? ? (t ) = e , e = u ? U = u (t ) ? U (t ).

N N N N N N N

Note that (5.3.5) runs backward in time starting at time t = tN .

5.3.

y

AN A POSTERIORI ERROR ESTIMATE FOR CG(1)

111

u(t) problem x 0 ? (t) problem T

Theorem 16. For N = 1, 2, . . . the cG(1) solution U (t) satis?es |eN | ≤ S (tN ) · max |k r (U )|,

[0,tN ]

(5.3.6)

where k = kn = |In | for t ∈ In = (tn?1 , tn ) is the time step and r (U ) = ˙ + aU ? f is the residual error. Further S (tN ), speci?ed below, is the U stability factor satisfying the quantitative bound ? tN if |a(t)| ≤ λ, ?t |? ˙ |dt ? eλ tN , S (tN ) := 0 ≤ (5.3.7) ? eN 1, if a(t) ≥ 0, ?t Proof. Let e(t) = u(t) ? U (t). Using the dual problem ?? ˙ (t) + a(t)?(t) = 0 we can write

2 2 e2 N = eN + 0 = eN + tN 0

e(?? ˙ + a?) dt,

tN tN

(5.3.8)

and by partial integration we get

tN 0

N e(?? ˙ + a(t)?)dt = [?e(t)?(t)]t 0 +

e? ˙ dt +

0 0

ea? dt

tN

tN

= ? e(tN ) ?(tN ) +

=e N =e N

(e ˙ + ae)? dt =

0

?e2 N

+

0

(e ˙ + ae)? dt,

where evaluating the boundary term we used e(0) = 0. Note that ˙ (t) + a(t)u(t) ? a(t)U (t), e ˙ (t) + a(t)e(t) = u ˙ (t) ? U and since f (t) = u ˙ (t) + a(t)u(t) we observe that ˙ + aU ? where the last equality is just the de?nition of the residual: r (U ) = U f . Consequently we get the error representation formula: e2 N = ?

tN

˙ (t) ? a(t)U (t) := ?r (U ), e ˙(t) + a(t)e(t) = f (t) ? U

(5.3.9)

r (U (t))?(t)dt.

0

(5.3.10)

112

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

1 kn In tN

To continue we use the interpolant πk ? = e2 N

tN

?(s)ds of ? and write r (U )πk ?(t)dt. (5.3.11)

=?

0

r (U )(?(t) ? πk ?(t))dt +

0

Now from the discrete variational formulation:

tN 0

˙ + aU )πk ?(t)dt = (U

0

tN

f πk ?(t)dt

(5.3.12)

we have the Galerkin orthogonality relation

tN

r (U )πk ?(t)dt = 0.

0

(5.3.13)

Thus the ?nal form of the error representation formula is e2 N = ?

tN 0

r (U )(?(t) ? πk ?(t))dt.

(5.3.14)

Now applying the interpolation error to the function ? in the interval In , |In | = kn we have

In

|? ? πk ?|dt ≤ kn

In

|? ˙ |dt.

(5.3.15)

This would yield the estimate

tN 0 N N In

|? ? πk ?|dt =

t∈J

n=1

|? ? πk ?|dt ≤

kn

n=1 In

|? ˙ |dt

(5.3.16)

Let now |v |J = max |v (t)|, then using (5.3.16) and the ?nal form of the error representation formula (5.3.14) we have that

N tN

|eN | ≤ Now since get

2

n=1

| r (U )| I n · k n

In

|? ˙ |dt ≤ max (kn |r (U )|In )

1≤n≤N

0

|? ˙ |dt.

tN 0

|?|dt = |eN | · S (tN ), (see the de?nition of S (tN )), we ?nally |eN |2 ≤ |eN |S (tN ) max(k |r (U )|).

[0,tN ]

(5.3.17)

This completes the proof of the ?rst assertion of the theorem.

5.3.

AN A POSTERIORI ERROR ESTIMATE FOR CG(1)

113

To prove the second assertion, we claim that: |a(t)| ≤ λ, |a(t)| ≥ 0, 0 ≤ t ≤ tN =? |?(t)| ≤ eλ tN |eN |, 0 ≤ t ≤ tN ?t =? |?(t)| ≤ |eN |, ?t ∈ [0, tN ]. (5.3.18) (5.3.19)

To prove this claim let s = tN ? t, (t = tN ? s) and de?ne ψ (s) = ?(tN ? s), then using the chain rule dψ dψ dt = · = ?? ˙ (tN ? s). ds dt ds The dual problem is now reformulated as ?nd ?(t) such that ?? ˙ (tN ? s) + a(tN ? s)?(tN ? s) = 0. The corresponding problem for ψ (s): ? ? dψ (s) + a(t ? s)ψ (s) = 0, t > s ≥ 0 N N ds ? ψ (0) = ?(tN ) = eN , eN = uN ? UN = u(tN ) ? U (tN ), (5.3.21) (5.3.20)

has the fundamental solution ψ (s) = CeA(tN ?s) , where ψ (0) = eN implies that C = e?A(tN ) eN and thus ψ (s) = eN e?A(tN ) eA(tN ?s) = eN eA(t)?A(tN ) . Now inserting back in the relation ψ (s) = ?(t), tN ? s = t, we get ?(t) = eN · eA(t)?A(tN ) , and ? ˙ (t) = eN · a(t)eA(t)?A(tN ) . (5.3.22)

Now the proof of both assertion in the claims are easily followed: (a) For |a(t)| ≤ λ, we have |?(t)| = |eN |e

Rt

tN

a(s)ds

≤ |eN |emaxt |a(t)|(tN ?t) ≤ |eN |eλ·tN

(5.3.23)

(b) For |a(t)| ≥ 0, we have |?(t)| = |eN |e

R tN

0

a(s)ds

≤ |eN |emint a(t)(t?tN )

(5.3.24)

and since (t ? tN ) < 0 we get that |?(t)| ≤ |eN |.

114

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

Now we return to the estimates for S (tN ). Note that for a(t) ≥ 0 we have using the second relation in (5.3.22) that

tN 0

˙t)|dt = |eN | | ?(

tN 0

N a(t)eA(t)?A(tN ) dt = |eN | · [eA(t)?A(tN ) ]t 0

= |eN | · 1 ? eA(0)?A(tN ) ≤ 1, which gives that S (tN ) = ˙t)|dt | ?( ≤ 1. |eN | As for the case |a(t)| ≤ λ, we use again (5.3.22): ? ˙ (t) = a(t)eN · eA(t)?A(tN ) and write

tN 0

|? ˙ (t)| ≤ λ|eN |e

tN 0

A(t)?A(tN )

= λ|eN |e

Rt

tN

a(s)ds

≤ λ|eN |eλ(tN ?t) .

(5.3.25)

Integrating over (0, tN ) we get

tN

|? ˙ (t)|dt ≤ |eN |

0

λeλ(tN ?t) dt = |eN | ? eλ(tN ?t)

tN 0

= |eN |(?1 + eλtN ),

which gives that S (tN ) ≤ (?1 + eλ·tN ) ≤ eλ·tN , and completes the proof of the second assertion. Theorem 17 ( Convergence order O(k 2 )). For N = 1, 2, . . . and with SN as in the previous theorem, the error for cG(1) solution U (t) satis?es |eN | ≤ S (tN ) max k 2 (aU ? f ) .

[0,tN ]

(5.3.26)

Proof. Using the orthogonality [g (t) ? πk g (t)] ⊥ (constants) ?g (t), and since ˙ (t) is constant on IN we have that tN U ˙ (? ? πk ?)dt = 0. Thus using error U 0 representation formula (5.3.14) yields e2 N = ? =

0 tN 0 tN tN

r (U )[?(t) ? πk ?(t)]dt =

tN

0

˙ )(? ? πk ?)dt (f ? aU ? U

(f ? aU )(? ? πk ?)dt ?

tN 0

0

˙ (? ? πk ?)dt U

=?

(aU ? f )(? ? πk ?)dt.

5.3.

AN A POSTERIORI ERROR ESTIMATE FOR CG(1)

115

g (t)

g (t) ? πk g (t)

Vh

πk g (t)

Figure 5.2: Orthogonality: (g (t) ? πk g (t)) ⊥ (constants) ?g (t). Similarly using the fact that πk (aU ? f ) is a constant we get

tN 0

πk (aU ? f )(? ? πk ?)dt = 0.

(5.3.27)

Consequently we can write e2 N = ?

tN 0

(aU ? f ) ? πk (aU ? f ) (? ? πk ?)dt.

(5.3.28)

Now using the above theorem and the interpolation error estimate we get |eN | ≤ S (tN ) · k |(aU ? f ) ? πk (aU ? f )| d ≤ S (tN ) · k (aU ? f ) dt

2 [0,tN ] [0,tN ]

(5.3.29)

.

116

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

5.4

A dG(0) a posteriori error estimate

|u(tN ) ? UN | ≤ S (tN )|kR(U )|[0,tN ] , UN = U (tN ) (5.4.1)

Theorem 18. For N = 1, 2, . . ., the dG(0) solution U (t) satis?es

where R (U ) = |UN ? UN ?1 | + |f ? aU | kn for tN ?1 < t ≤ tN . (5.4.2)

Proof. The proof uses similar techniques as in the cG(1) case. Note that here the residual error includes jump terms and since dual problem satis?es ?? ˙ (t) + a(t)?(t) = 0, we can write

N 2 e2 N = eN + n=1 N tn tn?1

e[?? ˙ (t) + a(t)?(t)]dt = [P I ] =

n (e ˙ + ae)?(t)dt ? [e?]t tn?1

tn

= =

e2 N e2 N

+

n=1 N tn?1 tn tn?1

(5.4.3)

N

+

n=1

(f ? aU )?dt ?

n [e?]t tn?1 ,

n=1

˙ + au ? aU = f ? aU and where in the last relation we use e ˙ + ae = u ˙ ?U ˙ also the fact that U = constant U = 0. We rewrite the last sum as follows

N

n (e?)t tn?1 =

N ? + + e(t? n )?(tn ) ? e(tn?1 )?(tn?1 )

n=1

n=1

? + + = {for a given functiong ; g (t? n ) = gn , g (tn?1 ) = gn?1 } N

=

n=1

+ + ? ? + + ? ? + + ? (e? n ?n ? en?1 ?n?1 ) = (e1 ?1 ? e0 ?0 ) + (e2 ?2 ? e1 e1 )

? + + ? ? + + + . . . + (e? N ?1 ?N ?1 ? eN ?2 ?N ?2 ) + (eN ?N ? eN ?1 ?N ?1 ). ? + + To continue for i = 1, . . . N ? 1, we write ?? i = (?i ? ?i + ?i ), then N

?

n=1

? ? + + ? ? + + + + n (e?)t tn?1 = ?eN ?N + e0 ?0 ? e1 (?1 ? ?1 + ?1 ) + e1 ?1 ? + + + + ? e? 2 (?2 ? +?2 + ?2 ) + e2 ?2 . . . ? + + + + ? e? N ?1 (?N ?1 ? ?N ?1 + ?N ?1 ) + eN ?1 ?N ?1 ,

5.4. A DG(0) A POSTERIORI ERROR ESTIMATE where a general i-th term can be rewritten as

? + + + + ? ? ? + ? + + + ? e? i (?i ? ?i + ?i ) + ei ?i = ?ei ?i + ?ei ?i ? ei ?i + ei ?i + ? + + ? ? + = e? i (?i ? ?i ) + ?i (ei ? ei ) = ei [?i ] + ?i [ei ],

117

with [g ] = g + ? g ? representing the jump. Hence we have

N

?

n=1

2 + + n (e?)|t tn?1 = ?eN + e0 ?0 +

N ?1 n=1

[en ]?+ n +

N ?1 n=1

e? n [?n ].

(5.4.4)

Inserting in (5.4.3) we get that

N tn tn?1 tn tn?1 N

e2 N

= =

e2 N e2 N

+

n=1 N

(f ? aU )?dt ? (f ? aU )?dt ?

tn

n [e?]t tn?1

n=1

+

n=1

e2 N

+

+ e+ 0 ?0

+

N ?1 n=1

[en ]?+ n

+

N ?1 n=1

[?n ]e? n =

= {?n , un smooth ? [?n ] = 0, [un ] = 0}

N

= =

+ e+ 0 ?0 N

+

n?1 tn tn?1 tn?1

(f ? aU )?dt +

N ?1 n=1

[en ]?+ n = {[un ] = 0 ? [en ] = [?Un ]}

n=1

(f ? aU )?dt ? [Un?1 ]?+ n?1 =

N tn tn?1

= {Galerkin} =

n=1

{(f ? aU )(? ? πk ?) ? [Un?1 ](? ? πk ?)+ n?1 }dt.

Now to continue we just follow the previous theorem. ?Adaptivity for dG(0) To guarantee that the dG(0) approximation U (t) satis?es |eN | = |u(tn ) ? U (tn )| ≤ T OL, (TOL is a given tolerance) (5.4.5)

we seek to determine the time step kn so that S (tN ) max |kn R(U )| = T OL,

t∈In

n = 1, 2, . . . , N.

(5.4.6)

118

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

?An adaptivity algorithm (i) Compute Un from Un?1 using a predicted step kn , for example

tn tn

aUn dt + Un =

tn?1 tn?1

f dt + Un?1 .

(5.4.7)

(ii) Compute |kR(U )|In := max |kn R(U )| and follow the chart:

n

Is (5.4.6) valid for this kn ? NO! Recompute (5.4.6) with a smaller kn

YES! —?→

Accept the solution Un and go to the next time step

5.5

A priori error analysis

?The discontinuous Galerkin method dG(0). The dG(0) method for u ˙ + au = f , a=constant, is formulated as follows: Find U = U (t), t ∈ In , such that

tn tn?1

˙ dt + a U

tn

Udt =

tn?1 In

f dt.

(5.5.1)

Note that U (t) = Un is constant for t ∈ In . Let Un = U (tn ), Un?1 = U (tn?1 ) and kn = tn ? tn?1 , then

tn tn?1

˙ dt + a U

tn tn?1

Udt = U (tn ) ? U (tn?1 ) + akn Un = Un ? Un?1 + akn Un .

Hence with a given initial data u(0) = u0 , the equation (5.5.1) is written as Un ? Un?1 + akn Un = f dt n = 1, 2, . . .

In

U0 = u0 .

(5.5.2)

For the exact solution u(t) of u ˙ + au = f , the same procedure yields

tn

u(tn ) ? u(tn?1 ) + kn aun (t) =

In

f dt + kn aun (t) ? a

u(t)dt,

tn?1

(5.5.3)

5.5. A PRIORI ERROR ANALYSIS

t

119

u(t)dt to the right hand side and add where we have moved the term a tnn ?1 kn aun (t) to both sides. Thus from (5.5.2) and (5.5.3) we have that (1 + kn a)Un (t) = Un?1 (t) + (1 + kn a)un (t) = un?1(t) + f dt,

In tn In

(5.5.4) u(t)dt.

f dt + kn aun (t) ? a

tn?1

Let now en = un ? Un and en?1 = un?1 ? Un?1 then (5.5.3) ? (5.5.2) yields en = (1 + kn a)?1 (en?1 + ρn )

tn tn?1

(5.5.5)

where ρn := kn aun (t) ? a Lemma 2. We have that

u(t)dt. Thus in order to estimate the error en

we need an iteration procedure and an estimate of ρn . 1 |ρn | ≤ |a||kn |2 max |u ˙ (t)| In 2 Proof. Recalling the de?nition we have ρn = kn aun (t) ? a |ρn | ≤ |a||kn | un ? 1 |kn | udt .

In tn tn?1

(5.5.6) u(t)dt. Thus (5.5.7)

Using a Taylor expansion of the integrand u(t) about tn , viz u(t) = un + u ˙ (ξ )(t ? tn ), we get that |ρn | ≤ |a||kn | un ? 1 [un + u ˙ (ξ )(t ? tn )]dt kn I n 1 (t ? tn )2 tn 1 ˙ (ξ ) ≤ |a||kn | un ? kn un ? u tn?1 kn kn 2 2 2 1 1 1 kn kn = |a||kn | ? u = |a||kn |2 |u = |a||kn | ? u ˙ (ξ ) 0 ? ˙ (ξ ) ˙ (ξ )| . kn 2 kn 2 2 1 ˙ (t)| |ρn | ≤ |a||kn |2 max |u In 2 for some ξ, tn?1 < ξ < tn (5.5.8)

Thus we have the following ?nal estimate for ρn , (5.5.9)

120

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

To simplify the estimate for en we split, and gather, the proof of technical details in the following lemma: Lemma 3. For kn |a| ≤ 1/2, n ≥ 1 we have that (i) (1 ? kn |a|)?1 ≤ e2kn |a| . (ii) Let τn = tN ? tn?1

N N

1 then |eN | ≤ 2

tN 0

n=1

(e2|a|τn |a|kn ) max kn |u ˙ |In .

1≤n≤N

(iii)

n=1

e2|a|τn |a|kn ≤ e

|a|e2|a|τ dτ .

We postpone the proof of this lemma and ?rst show that using these results we can obtain a bound for the error eN (our main result) viz, Theorem 19. If kn |a| ≤ 1 , n ≥ 1 then the error of the dG(0) approximation 2 U satis?es |u(tN ) ? U (tN )| = |eN | ≤ e 2|a|tN e ?1 4

1≤n≤N

max kn |u ˙ (t)|In .

(5.5.10)

Proof. Using the estimates (ii) and (iii) of the above lemma we have that 1 |eN | ≤ 2

N

1 e = e 2 2

n=1 2|a|τ

(e2|a|τn |a|kn ) max kn |u ˙ |In ≤

1≤n≤N tN 0

1 e 2

tN 0

|a|e2|a|τ dτ

1≤n≤N

1≤n≤N

max kn |u ˙ |In

· max kn |u ˙ (t)|In =

1≤n≤N

e 2|a|tN e ?1 4

max kn |u ˙ (t)|In .

e 2|a|tN e ? 1 may grow depending on 4 |a| and tN , and then this result may not be satisfactory at all. Now we return to the proof of our technical results: Note that the stability constant Proof of Lemma 3. (i) For 0 ≤ x := kn |a| ≤ 1/2, we have that 1/2 ≤ 1 ? x < 1 and 0 < 1 ? 2x ≤ 1. We can now multiply both side of the ?rst claim: 1 < e2x by 1 ? x ≥ 1/2 > 0 to obtain the equivalent relation 1?x f (x) := (1 ? x)e2x > 1. (5.5.11)

5.5. A PRIORI ERROR ANALYSIS

121

Note that since f (0) = 1 and f ′ (x) = (1 ? 2x)e2x > 0 the relation (5.5.11) is valid. (ii) We recall that en = (1 + kn a)?1 (en?1 + ρn ). To deal with the coe?cient (1 + kn a)?1 ?rst we note that (1 + kn a)?1 ≤ (1 ? kn a)?1 if a ≥ 0. Thus 1 (1 + kn |a|)?1 ≤ (1 ? kn |a|)?1 , a ∈ R. Further the assumption kn |a| ≤ for 2 n ≥ 1, combined with (i), implies that (1 ? kn |a|)?1 ≤ e2kn |a| , n ≥ 1. Thus 1 1 |eN ?1 | + |ρN | ≤ |eN ?1 | · e2kN |a| + |ρN | · e2kN |a| . 1 ? kN |a| 1 ? kN |a| (5.5.12) Relabeling, e.g. N to N ? 1 we get |eN | ≤ |eN ?1 | ≤ |eN ?2 | · e2kN ?1 |a| + |ρN ?1 | · e2kN ?1 |a| = e2kN ?1 |a| |eN ?2 | + |ρN ?1 | , which, inserting in (5.5.12) gives that |eN | ≤ e2kN |a| e2kN ?1 |a| |eN ?2 | + |ρN ?1 | + |ρN | · e2kN |a| . (5.5.13)

Similarly we have |eN ?2 | ≤ e2kN ?2 |a| |eN ?3 | + |ρN ?2| . Now iterating (5.5.13) and using the fact that e0 = 0 we get, |eN | ≤e2kN |a| e2kN ?1 |a| e2kN ?2 |a| |eN ?3 | + e2kN |a| e2kN ?1 |a| e2kN ?2 |a| |ρN ?2 | + e2kN |a| e2kN ?1 |a| |ρN ?1 | + |ρN | · e2kN |a| ≤ · · · ≤

N N 2|a| PN

n=1

≤e

kn

|e0 | +

e

n=1

2|a|

PN

m=n

km

|ρn | =

e2|a|

n=1

PN

m=n

km

|ρn |.

1 ˙ (t)|. Thus Recalling (5.5.6) (Lemma 2): |ρn | ≤ |a||kn |2 max |u In 2

N

|eN | ≤ Note that

N

e2|a|

n=1

PN

m=n

km 1

2

˙ (t)|. |a||kn |2 max |u

In

(5.5.14)

m=n

km = (tn ? tn?1 )+(tn+1 ? tn )+(tn+2 ? tn+1 )+ . . . +(tN ? tN ?1 ) = tN ? tn?1 .

122

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

Hence we have shown the assertion (ii) of the lemma, i.e.

N

|eN | ≤

e

n=1

2|a|(tN ?tn?1 ) 1

1 ˙ (t)| = |a||kn | max |u In 2 2

2

N

n=1

(e2|a|τn |a|kn ) max kn |u ˙ |In .

1≤n≤N

(iii) To prove this part we note that τn = tN ? tn?1 = (tN ? tn ) + (tn ? tn?1 ) = τn+1 + kn , (5.5.15)

and since |a|kn ≤ 1/2 we have 2|a|τn = 2|a|τn+1 + 2|a|kn ≤ 2|a|τn+1 + 1. Further for τn+1 ≤ τ ≤ τn , we can write

τn τn

e2|a|τn · kn = =

τn+1 τn τn+1

e2|a|τn dτ ≤ e ·e

1 2|a|τn+1

e(2|a|τn+1 +1) dτ

τn+1 τn

(5.5.16) e

2|a|τ

dτ ≤ e

dτ.

τn+1

Multiplying (5.5.16) by |a| and summing over n we get

N N τn τn+1

e

n=1

2|a|τn

|a|kn ≤ e =e

n=1 τ1 τN +1

e2|a|τ dτ |a|

tN 0

(5.5.17) |a|e2|a|τ dτ,

e2|a|τ |a|dτ = e

which is the desired result and the proof is complete.

5.6

The parabolic case (a(t) ≥ 0)

We state and proof the basic estimate of this case Theorem 20. Consider the dG(0) approximation U for u ˙ + au = f , with 1 a(t) ≥ 0. Assume that kj |a|Ij ≤ , ?j , then we have the error estimates 2 ? ? ? 3e2λtN max |k u ˙ | if |a(t)| ≤ λ 0≤t≤tN |u(tN ) ? UN | ≤ (5.6.1) ? ? 3 max |k u ˙| if a(t) ≥ 0.

0≤t≤tN

5.6.

THE PARABOLIC CASE (A(T ) ≥ 0)

123

Sketch of the proof. Let e = u ? U = (u ? πk u) + (πk u ? U ) := e ?+ e ?, where (0) e ? is the interpolation error with πk u being the L2 -projection into Wk . To estimate e ?, we shall use the following discrete dual problem (DDP): (0) Find Φ ∈ Wk , such that for n = N, N ? 1, . . . , 1. ? tn ? ˙ + a(t)Φ)vdt ? [Φn ]vn = 0, ?v ∈ W (0) ? (?Φ k tn?1 (5.6.2) (DDP ) ? ? Φ+ = Φ ?N . N +1 = (πk u ? U )N :≡ e N Let now v = e, then

N tn tn?1

|e ?N | =

2

n=1

˙ + a(t)Φ)? (?Φ edt ?

N ?1 n=1

[Φn ]? en + ΦN e ?N .

(5.6.3)

We now use e ? = (πk u ? U ) = (πk u ? u + u ? U ) and write (5.6.3) as

N tn

|eN | =

2

n=1 tn?1 N ?1

˙ + a(t)Φ](πk u ? u + u ? U )UT [?Φ

?

n=1

[Φn ](πk u ? u + u ? U )n + ΦN (πk u ? u + u ? U )N .

Using Galerkin orthogonality we replace u by U . Therefore the total contribution from the terms with the factor u ? U is identical to zero. Thus due ˙ = 0 on each subinterval, we have the error representation to the fact that Φ formula:

N tn tn?1 tN

|eN | =

2

n=1

˙ + a(t)Φ)(πk u ? u)dt ? (?Φ

N ?1 n=1

N ?1 n=1

[Φn ](πk u ? u)n + ΦN (πk u ? u)N

=

0

(a(t)Φ)(u ? πk u) dt +

[Φn ](u ? πk u)n ? ΦN (u ? πk u)N .

To continue we shall need the following results: Lemma 4. If |a(t)| ≤ λ, ?t ∈ (0, tN ) and kj |a|Ij ≤ 1 , j = 1, 2, . . . , N , then 2 the solution of the discrete dual problem satis?es (i) |Φn | ≤ e2λ(tN ?tn?1 ) |e ?N |.

124 (ii)

N ?1 n=1 N

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM |[Φn ]| ≤ e2λtN |e ?N |.

tn tn?1

(iii)

n=1

a(t)|Φn |dt ≤ e2λtN |e ?N |.

(iv) If a(t) ≥ 0 then Max |Φn |,

N ?1 n=1 N tn tn?1

|[Φn ]|,

n=1

a(t)|Φn | dt ≤ |e ?N |.

Proof. We show the last estimate (iv), (the proofs of (i)-(iii) are similar to that of the stability factor in the previous theorem). Consider the discrete dual problem with v ≡ 1: ? tn ? ˙ + a(t)Φ)dt ? [Φn ] = 0, ? (?Φ tn?1 (5.6.4) (DDP ) ? ? Φ ?N . N +1 = (πk u ? U )N :≡ e For dG(0) this becomes ? tn ? ?Φ n+1 + Φn + Φn tn?1 a(t) = 0, (DDP ) ? Φ ?N , Φn = Φ|In . N +1 = e By iterating we get

N

n = N, N ? 1, . . . , 1 (5.6.5)

Φn =

j =n

1+

Ij ?1

a(t)dt

?1

ΦN +1

(5.6.6)

For a(t) ≥ 0 we have 1 +

Ij

a(t)dt

≤ 1, thus (5.6.6) implies that (5.6.7)

|Φn | ≤ ΦN +1 = |e ?N |. Further we have using (5.6.6) that

N

Φn?1 =

1+

j =n?1 Ij

a(t)dt

?1

ΦN +1 = 1 +

In?1

a(t)dt

?1

Φn ≤ Φn

5.6.

THE PARABOLIC CASE (A(T ) ≥ 0)

125

which implies that

? [Φn ] = Φ+ n ? Φn = Φn+1 ? Φn ≥ 0.

(5.6.8)

Thus

N

n=1

|[Φn ]| = ΦN +1 ? ΦN + ΦN ? ΦN ?1 + . . . + Φ2 ? Φ1 = ΦN +1 ? Φ1 ≤ ΦN +1 ≤ |e ?N |.

(5.6.9)

Finally in the discrete equation:

tn tn?1

˙ + a(t)Φ)vdt ? [Φn ]vn = 0, (?Φ

?v ∈ Wk

(0)

(5.6.10)

˙ ≡ 0 for the dG(0). Hence (5.6.10) can be rewritten as we have v ≡ 1 and Φ

tn

a(t)Φn dt = [Φn ].

tn?1

(5.6.11)

Summing over n, this gives that

N tn tn?1 N

n=1

a(t)Φn dt ≤

n=1

[Φn ] ≤ |e ?n |.

(5.6.12)

Combining (5.6.7), (5.6.9), and (5.6.12) the proof of (iv) is now complete. ?Quadrature rule for f : Assume that a=constant. Then the error representation formula, combining dG(0), with the quadrature role for f is as follows:

N tn tn?1 tn

e2 N =

n=1

(f ? aU )(? ? πk ?)dt ? [Un?1 ](? ? πk ?)+ n?1 (5.6.13)

+

tn?1

f πk ?dt ? (f πk ?)n kn

quadrature error

126

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

where for the endpoint-rule g ?n = g (tn ), whereas for the midpoint-rule g ?n := tN ?(tN ) := 0 |?|dt , where g (t(n?1/2) ). We also de?ne the weak stability factor S |eN | ? is the solution of the dual problem ?? ˙ + a? = 0, for tN > t ≥ 0 ?(tN ) = eN .

Note that πk ? is piecewise constant and

In

|πk ?(t)|dt ≤

In

|?(t)|dt.

We can prove the following relations between the two stability factors: ?(tN ) ≤ tN (1 + S (tN )). S ?(tN ) >> S (tN ). Note that if a > 0 is su?ciently small, then S Theorem 21 (The modi?ed a posteriori estimate for dG(0)). The dG(0) approximation U (t) computed using quadrature on terms involving f satis?es for N = 1, 2, . . . ?(tN )Cj |k j f (j ) |(0,t ) , |u(tn ) ? Un | ≤ S (tn )|kR(U )|(0,tN ) + S N where R (U ) = (5.6.14)

|Un ? Un?1 | + |f ? aU |, on In (5.6.15) kn and j = 1 for the rectangle rule, j = 2 for the midpoint rule, C1 = 1, C2 = 1 ¨. , f (1) = f˙ and f (2) = f 2

5.6.1

Short summary of error estimates

In this part we shall derive some short variants for the error estimates above Lemma 5. Let U be the cG(1) approximation of u satisfying the initial value problem u ˙ + u = f, t > 0, u(0) = u0 . (5.6.16) Then we have that ˙ ? U )| , |(u ? U )(T )| ≤ max |k (f ? U

[0,T ]

(5.6.17)

where k is the time step.

5.6.

THE PARABOLIC CASE (A(T ) ≥ 0)

127

Proof. The error e = u ? U satis?es Galerkin orthogonality:

T

(e ˙ + e)vdt = 0,

0

for all piecewise constants v (t).

(5.6.18)

Let ? satisfy the dual equation ?? ˙ + ? = 0, Then we have that gives t < T, ?(T ) = e(T ). (5.6.19)

?(t) = e(T ) · et?T : Note that integrating ?? ˙ +? = 0 (5.6.20)

? ˙ dt = 1 · dt =? ln ? = t + C. ? Let now C = ln C1 , then (5.6.20) can be written as ln ? ? ln C1 = ln ? = t =? ?(t) = C1 · et . C1

(5.6.21)

Finally, since ?(T ) = e(T ) we have that C1 · eT = e(T ), i.e. C1 = e(T ) · e?T =? ?(t) = e(T ) · et?T .

T T

(5.6.22)

To continue we have using ?? ˙ + ? = 0,

T

|e(T )|2 = e(T ) · e(T ) +

0

e(?? ˙ + ?)dt = e(T ) · e(T ) ?

e? ˙ dt +

0 0

e?dt.

Note that integration by parts gives

T 0 T T

e?dt ˙ = [e ·

?]T t=0

?

0

e?dt ˙ = e(T )?(T ) ? e(0)?(0) ?

e?dt. ˙

0

Using ?(T ) = e(T ), and e(0) = 0, we thus have

T T T

|e(T )|2 = e(T ) · e(T ) ? e(T ) · e(T ) +

T T

e? ˙ dt +

0 0

e? dt =

0

(e ˙ + e)? dt

=

0

(e ˙ + e)(? ? v )dt =

0

˙ ? U (? ? v )dt. u ˙ + u ?U

=f

˙ + U ? f := r (U ), is the residual and We have that U

T T

|e(T )|2 = ?

0

r (U ) · (? ? v )dt ≤ max |k · r (U )|

[0,T ]

0

1 |? ? v |dt. (5.6.23) k

128 Recall that

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

I

h?1 |? ? v |dx ≤

T

I

|?′ |dx.

T

(5.6.24)

Further ?? ˙ + ? = 0 implies ? ˙ = ?, and ?(t) = e(T ) · et?T . Thus |e(T )|2 ≤ max |k · r (U )|

[0,T ] [0,T ] 0

|? ˙ |dt = max |k · r (U )|

[0,T ] T

0

|?(t)| dt

(5.6.25)

≤ max |kr (U )|e(T )| and since

T 0

e

0

t?T

dt,

0 ?T et?T dt = [et?T ]T = 1 ? e?T ≤ 1, 0 = e ?e

T > 0,

we ?nally end up with the desired result |e(T )| ≤ max |k · r (U )|.

[0,T ]

Problem 45. Generalize the Lemma to the problem u ˙ + au = f , with a = positive constant. Is the statement of Lemma 1 valid for u ˙ ? u = f? Problem 46. Study the dG(0)-case for u ˙ + au = f, a>0 Lemma 6. Let u ˙ + u = f, t > 0. Show for the cG(1)-approximation U (t) that |(u ? U )(T )| ≤ max |k 2 u ¨|T. (5.6.26)

[0,T ]

Sketchy proof, via the dual equation. Let ? be the dual solution satisfying ? ˙ + ? = 0, t < T, We compute the error at time T , viz

T

?(T ) = e(T ).

T 0 T 0

|e(T )| = |Θ(T )| = Θ(T )?(T ) +

T

2

2

0 T 0

? ?Φ ˙ + Φ) dt = Θ(

=0

˙ + Θ)Φ ? dt (Θ

=?

0

? dt = ? (ρ ˙ + ρ)Φ

? dt ≤ max |k 2 u ρ·Φ ¨|

[0,T ]

? | dt |Φ

≤ max |k 2 u ¨| · T · |e(T )|.

[0,T ]

5.7. EXERCISES

129

Here ρ = u ? u ?, Θ = u ? ? U and Φ is cG(1)-approximation of φ such that T ˙ v (?Φ + Φ) dt = 0 for all piecewise constant v (t). Furthermore u ? is 0 the piecewise linear interpolant of u and w ? = is the piecewise constant mean value.

5.7

Exercises

Problem 47. (a) Derive the sti?ness matrix and load vector in piecewise polynomial (of degree q ) approximation for the following ODE in population dynamics, ? ? u ˙ (t) = λu(t), for 0 < t ≤ 1, ? u(0) = u .

0

(b) Let λ = 1 and u0 = 1 and determine the approximate solution U (t), for q = 1 and q = 2. Problem 48. Consider the initial value problem u ˙ (t) + a(t)u(t) = f (t), 0 < t ≤ T, u(0) = u0 .

Show that for a(t) > 0, and for N = 1, 2, . . . , the piecewise linear approximate solution U for this problem satis?es the a posteriori error estimate ˙ + aU ? f )|, |u(tN ) ? UN | ≤ max |k (U

[0,tN ]

k = kn , for tn?1 < t ≤ tn .

Problem 49. Consider the initial value problem: u ˙ (t) + au(t) = 0, t > 0, u(0) = 1.

a) Let a = 40, and the time step k = 0.1. Draw the graph of Un := U (nk ), k = 1, 2, . . . , approximating u using (i) explicit Euler, (ii) implicit Euler, and (iii) Crank-Nicholson methods. b) Consider the case a = i, (i2 = ?1), having the complex solution u(t) = e?it with |u(t)| = 1 for all t. Show that this property is preserved in CranckNicholson approximation, (i.e. |Un | = 1 ), but NOT in any of the Euler approximations.

130

CHAPTER 5. SCALAR INITIAL VALUE PROBLEM

Problem 50. Consider the initial value problem u ˙ (t) + au(t) = 0, t > 0, u(0) = u0 , (a = constant).

Assume a constant time step k and verify the iterative formulas for dG(0) ? , respectively: i.e. and cG(1) approximations U and U Un = 1 1 + ak

n

u0 ,

?n = 1 ? ak/2 U 1 + ak/2

n

u0 .

Problem 51. Let U be the cG(1) approximation of u satisfying the initial value problem u ˙ + au = f, t > 0, u(0) = u0 . Let k be the time step and show that for a = 1, ˙ ? U )||L∞ [0,T ] , T ||k 2u |(u ? U )(T )| ≤ min ||k (f ? U ¨||L∞ [0,T ] . Problem 52. Consider the scalar boundary value problem u ˙ (t) + a(t)u(t) = f (t), t > 0, u(0) = u0 .

(a) Show that for a(t) ≥ a0 > 0, we have the stability estimate

t

|u(t)| ≤ e

?a0 t

|u 0 | +

0

ea0 s |f (s)| ds

(b) Formulate the cG(1) method for this problem, and show that the con1 a0 k > ?1, where k is the time step, guarantees that the method is dition 2 operational, i.e. no zero division occurs. (c) Assume that a(t) ≥ 0, f (t) ≡ 0, and estimate the quantity Problem 53. Consider the initial value problem (u = u(x, t)) u ˙ + Au = f, t > 0; u(t = 0) = u0 .

RT

0

|u ˙ | dt . | u0 |

Show that if there is a constant α > 0 such that (Av, v ) ≥ α||v ||2,

t

?v, 1 α

t 0

then the solution u of the initial value problem satis?es the stability estimate ||u(t)||2 + α ||u(s)||2 ds ≤ ||u0||2 + ||f (s)||2 ds.

0

Chapter 6 The heat equation in 1d

In this chapter we focus on some basic stability and ?nite element error estimates for the one-space dimensional heat equation. A general discussion on classical heat equation can be found in our Lecture Notes in Fourier Analysis. We start our study considering an example of an initial boundary value problem with mixed boundary conditions. Higher dimensional case is considered in forthcoming lecture notes based on the present one. Here we consider an example of an initial boundary value problem for the heat equation, viz ? ? ? u ˙ ? u′′ = f (x), 0 < x < 1, t > 0, ? ? (6.0.1) (IBV P ) u(x, 0) = u0 (x), 0 < x < 1, ? ? ? ? u(0, t) = u (1, t) = 0, t > 0, x u ˙ := ut = ?u , ?t u′ := ux = ?u , ?x u′′ := uxx = ?2u . ?x2

where we have used the following di?erentiation notation in the 1 ? D case:

Note that the partial di?erential equation in (6.0.1) containing three derivatives yields three degrees of freedom and therefore, to determine a unique solution, it is necessary to supply three data: here two boundary condition associated to two spatial derivatives (in u′′ ) and an initial condition corresponding to the time derivative (u ˙ ). To have an idea we formulate an example, viz

131

132 u(x, t)

CHAPTER 6. THE HEAT EQUATION IN 1D

u0 (x) x tn?1 tn t Figure 6.1: A decreasing temperature pro?le

Problem 54. Give physical meaning to the IBVP (6.0.1) where f = 20 ? u. solution: Heat conduction with u(x, t) = temperature at x at time t. u(x, 0) = u0 (x), the initial temperature at time t = 0. u(0, t) = 0, u′ (1, t) = 0, f = 20 ? u, ?xed temperature at time x = 0. isolated boundary at x = 1 (no hear ?ux). heat source, in this case a control system to force u ? 20.

6.1

Stability estimates

In this part we shall derive a general stability estimate for the mixed IBVP above, prove a 1 ? D version of the Poincare inequality and then derive some homogeneous stability estimates. Theorem 22. The IBVP (6.0.1) satis?es the stability estimates

t

||u(·, t)|| ≤ ||u0|| + ||ux(·, t)||2 ≤ ||u′0 ||2 +

0

||f (·, s)|| ds,

t 0

(6.1.1) (6.1.2)

||f (·, s)||2 ds.

6.1. STABILITY ESTIMATES

133

Proof. Multiply the equation in (6.0.1) by u and integrate over (0, 1) to get

1 0 1 1

uu ˙ dx ?

u′′ u dx =

0 0

f u dx.

(6.1.3)

Integrating by parts we get 1d 2 dt

1 1 1

u2 dx +

0 0

(u′)2 dx ? u′ (1, t)u(1, t) + u′ (0, t)u(0, t) =

f u dx,

0

and using the boundary conditions and the Cauchy-Schwartz inequality we end up with 1 d ′ 2 f u dx ≤ ||f ||||u||. (6.1.4) ||u|| ||u|| + ||u || = dt 0 Consequently ||u|| d ||u|| ≤ ||f ||||u||, dt and thus d ||u|| ≤ ||f ||. dt (6.1.5)

Integrating over time we get

t

||u(·, t)|| ? ||u(·, 0)|| ≤

0

||f || ds,

(6.1.6)

which gives (6.1.1). To prove (6.1.2) we multiply the di?erential equation by u ˙ and integrate over (0, 1) to obtain

1 0 1 1

(u ˙ )2 dx ?

0

u′′ u ˙ dx = ||u ˙ ||2 +

1

0

u′ u ˙ ′ dx ? u′ (1, t)u ˙ (1, t) + u′(0, t)u ˙ (0, t)

=

0

fu ˙ dx.

The expression above gives ||u ˙ ||2 + Thus and hence 1d ′ 2 ||u || = 2 dt

1 0

fu ˙ dx ≤ ||f ||||u ˙ || ≤

1 ||f ||2 + ||u′||2 . 2

(6.1.7)

1d ′ 2 1 1 ||u ˙ ||2 + ||u || ≤ ||f ||2, 2 2 dt 2 d ′ 2 ||u || ≤ ||f ||2. dt

(6.1.8) (6.1.9)

134

CHAPTER 6. THE HEAT EQUATION IN 1D

Now integrating over (0, t) we get the desired result

t

||u′(·, t)||2 ? ||u′(·, 0)||2 ≤ and the proof is complete.

0

||f (·, s)||2 ds,

(6.1.10)

To continue we prove a one-dimensional version of the one of the most important inequalities in PDE and analysis. Theorem 23 (Poincare inequality in 1 ? D case). Assume that u and u′ are square integrable. There exists a constant C , independent of u but dependent of L, such that if u(0) = u(L) = 0, then there is constant C , independent of u but dependent of L, such that

L 0 L

u(x) dx ≤ C

2

u′ (x)2 dx,

0

i.e. ||u|| ≤

√

C ||u′||.

(6.1.11)

Proof. Note that we can successively write

x x x

u(x) =

0

u′ (y ) dy ≤

L 0

0

|u′ (y )| dy ≤

1/2 L

0

|u′(y )| · 1 dy

1/2

≤ Thus

L 0

|u′ (y )|2 dy

L

·

12 dy

0

=

√

L

L

0

|u′(y )|2 dy

1/2

.

L

L

u(x)2 dx ≤

L

0 0

|u′ (y )|2 dy = L2 ||u|| ≤ L||u′||.

0

|u′(y )|2 dy,

(6.1.12)

and hence (6.1.13)

Remark 13. The constant c = L means that the Poincare inequality is valid for arbitrary bounded intervals, but not! for unbounded intervals. It is also unnecessary to have both boundary values equal zero. For instance if v (0) = 0 and, for simplicity L = 1, then by the same argument as above we get the following version of one-dimensional Poincare’s’ inequality:

2 ′ 2 ||u||2 L2(0,1) ≤ 2 v (0) + ||u ||L2 (0,1) .

(6.1.14)

6.1. STABILITY ESTIMATES

135

Theorem 24 (Stability of the homogeneous heat equation). The homogeneous INBVP for the heat equation ? ? ? u ˙ ? u′′ = 0, 0 < x < 1, t>0 ? ? (6.1.15) u(0, t) = ux (1, t) = 0, t>0 ? ? ? ? u(x, 0) = u (x), 0 < x < 1,

0

satis?es the stability estimates a)

d ||u||2 + 2||u′||2 = 0, dt

1 1

b) ||u(·, t)|| ≤ e?t ||u0 ||.

Proof. a) Multiply the equation by u and integrate over x ∈ (0, 1),

1

0=

0

(u ˙ ? u′′ )u dx =

uu ˙ dx +

0 0

(u′ )2 dx ? u′(1, t)u(1, t) + u′(0, t)u(0, t).

Using integration by parts and the boundary data we get 1d 2 dt

1 1

u2 dx +

0 0

(u′)2 dx =

d ||u||2 + 2||u′||2 = 0. dt

This gives the proof of a). As for b) using a) together with the Poincare inequality with L = 1: ||u|| ≤ ||u′|| we have that d ||u||2 + 2||u||2 ≤ 0. dt Multiplying both sides of (6.1.16) by e2t yields d d ||u||2e2t ≤ ||u||2 + 2||u||2 e2t ≤ 0. dt dt We replace t by s and integrate over s ∈ (0, t) to obtain

t 0

(6.1.16)

(6.1.17)

d ||u||2e2s ds = ||u(·, t)||2e2t ? ||u(·, 0)||2 ≤ 0. ds

(6.1.18)

This yields ||u(·, t)||2 ≤ e?2t ||u0||2 =? ||u(·, t)|| ≤ e?t ||u0 ||, and completes the proof. (6.1.19)

136

CHAPTER 6. THE HEAT EQUATION IN 1D

Remark 14. For the sake of generality and application of this technical argument in higher dimensions we shall use a general notation for the domain and its boundary: namely ? and ? ? respectively. The reader may replace ? by any interval (a, b), for instance I = (0, 1) and ? ? by the corresponding boundary. The proof of the general theorem for the energy estimate in higher dimensions is given in part II. Theorem 25 (An energy estimate). For any small ε > 0 We have that

t ε

u ˙ (s)ds ≤

1 2

ln

t u0 . ε

(6.1.20)

Proof. Multiply the di?erential equation: u ˙ ? u′′ = 0, by ?tu′′ and integrate over ? to obtain ?t uu ˙ ′′ dx + t

? ?

?(u′′ )2 dx = 0.

(6.1.21)

Integrating by parts and using the fact that u = 0 on ? ? we get uu ˙ ′′ dx = ? u ˙ ′ · u′ dx = ? 1d ′ 2 u , 2 dt (6.1.22)

?

?

so that (11.1.11) can be written as t 1d ′ u 2 dt

2

+ t u′′

2

= 0, =

d (t dt

(6.1.23) u′ 2 ) ? u′

2

d u′ and by using the obvious relation t dt

2

we get (6.1.24)

d (t u′ 2 ) + 2t u′′ dt

2

= u′ 2 .

We now change t to s and integrate over s ∈ (0, t) to get

t 0

d (s u′ 2 (s)) ds + 2 ds

t

t

s u′′ 2 (s)ds =

0 0

u′ 2 (s)ds ≤

1 u0 2 , 2

where in the last inequality we just integrate the stability estimate (a) in the previous theorem. Consequently

t

t u′ 2 (t) + 2

0

s u′′ 2 (s) ds ≤

1 u0 2 . 2

(6.1.25)

6.2. FEM FOR THE HEAT EQUATION In particular, we have: (I ) 1 u (t) ≤ √ u0 2t

′ t

137

(II )

0

s u′′ 2 (s) ds

1/2

≤

1 u0 2

(6.1.26)

Analogously we can show that 1 u′′ (t) ≤ √ u0 2t (6.1.27)

Now using the di?erential equation u ˙ = u′′ and integrating (6.1.27) we obtain

t ε

1 u ˙ (s)ds ≤ √ u0 2

t ε

t 1 1 ds = √ ln u0 s 2 ε

(6.1.28)

or more carefully

t t t t

u ˙ (s)ds =

ε ε

u′′ (s)ds =

ε t

1 · u′′ (s)ds =

t

0 1/2

1 √ ε √ · s u′′ (s)ds s ≤ 1 2 ln t u0 , ε

≤

s?1 ds

ε

1/2

·

s u′′ 2 (s) ds

ε

where in the ?rst inequality is just an application of the Cauchy Schwartz inequality and the second is an application of (6.1.26) (II) and we have obtained the desired result. Problem 55. Prove (6.1.27). Hint: Multiply (1) by t2 (u′′ )2 and note that u′′ = u ˙ = 0 on ? ?, or alternatively: di?erentiate u ˙ ? u′′ = 0 with respect to 2 t and multiply the resulting equation by t u ˙.

6.2

FEM for the heat equation

Consider the one-dimensional heat equation with Dirichlet boundary condition ? ? ? u ˙ ? u′′ = f, 0 < x < 1, t > 0, ? ? (6.2.1) u(0, t) = u(1, t) = 0, t > 0, ? ? ? ? u(x, 0) = u (x), 0 < x < 1.

0

138

CHAPTER 6. THE HEAT EQUATION IN 1D

The Variational formulation for the problem (6.2.1) reads as follows: For every time interval In = (tn?1 , tn ] ?nd u(x, t), t ∈ In such that

1 1

(uv ˙ + u′ v ′ )dxdt =

In 0 In 0

f vdxdt,

?v : v (0, t) = v (1, t) = 0. (VF)

A piecewise linear Galerkin approximation: For each time interval In = (tn?1 , tn ], with tn ? tn?1 = k , let U (x, t) = Un?1 (x)Ψn?1 (t) + Un (x)Ψn (t), where Ψn (t) = and Un (x) = Un,1 ?1 (x) + Un,2 ?2 (x) + . . . + Un,m ?m (x), (6.2.4) t ? tn?1 , k Ψn?1 (t) = tn ? t , k k = tn ? tn?1 , (6.2.3) (6.2.2)

with ?(xj ) = δij being the usual ?nite element basis corresponding to a partition of ? = (0, 1), with 0 = x1 < · · · < xk < xk+1 < · · · < xm = 1. In other words U is piecewise linear in both space and time variables and the unknowns are the coe?cients Un,k satisfying the following discrete variational formulation:

1 In 0

˙ ?j + U ′ ?′ ) dxdt = (U j

1

f ?j dxdt,

In 0

j = 1, 2 , . . . , m

(6.2.5)

Note on In = (tn?1 , tn ] and with Un := U (xn ) and Un?1 := U (xn?1 ) we have ˙ (x, t) = Un?1 (x)Ψ ˙ n?1 (t) + Un (x)Ψ ˙ n (t) = Un ? Un?1 . U k (6.2.6)

Further di?erentiating (6.2.2) with respect to x we get

′ ′ U ′ (x, t) = Un ?1 (x)Ψn?1 (t) + Un (x)Ψn (t).

(6.2.7)

6.2. FEM FOR THE HEAT EQUATION

139

U n+1 (x)

y ψ (t ) n

ψn+1(t)

t t n+1 tn t n-1

x i-1

Un(x)

1

? (x)

i

xi

x i+1

x

Inserting (6.2.6) and (6.2.7) in (6.2.5) we get using that Ψn?1 dt = k 2 In

1 0 1

In

dt = k and

In

Ψn dt =

1

Un ?j dx ?

M · Un

0

Un?1 ?j dx +

M ·Un?1 1

In

Ψn?1 dt

k 2

0

′ ′ Un ?1 ?j dx S ·Un?1

1 ′ ′ Un ?j S · Un

(6.2.8)

+

In

Ψn dt

0

k 2

dx =

In 0

f ?j dxdt

F

which can be written in a compact form as the Crank- Nicholson system (CNS) k k (CNS) M + S Un = M ? S Un?1 + F, 2 2 with the solution Un given by k Un = M + S 2

B ?1 ?1

k k M ? S Un?1 + M + S 2 2

A B ?1

?1

F,

(6.2.9)

140 where

CHAPTER 6. THE HEAT EQUATION IN 1D ? ?

Thus with a given source term f we can determine the source vector F and then, for each n = 1, 2, . . . N , given the vector Un?1 we use the CNS to compute the m-dimensional vector Un (m nodal values of U at the time level tn ). Example 26. Derive a corresponding equation system, as above, for the dG(0). The matrices S and M introduced in (6.2.8) are known as the sti?ness matrix and Mass matrix respectively. Below we compute these matrices. Note that di?erentiating (6.2.4): Un (x) = Un,1 ?1 (x) + Un,2 ?2 (x) + . . . + Un,m ?m (x), we get

′ Un (x) = Un,1 ?′1 (x) + Un,2 ?′2 (x) + . . . + Un,m ?′m (x).

U ? n,1 ? ? Un,2 Un = ? ? ? ... ? Un,m

? ? ? ?. ? ? ?

(6.2.10)

(6.2.11)

Thus for j = 1, . . . , m we have

1 1 ′ ′ Un ?j = 0 1 1

SUn =

0

?′j ?′1 Un,1 +

0

?′j ?′2 Un,2 +. . .+

0

?′j ?′m Un,m ,

which can be written in the matrix form as ? 1 ′ ′ 1 ′ ′ 1 ?? ?1 ?2 . . . 0 ?′1 ?′m 0 ? 0 1 1 ? 1 ′ ′ 1 ′ ′ 1 ? 0 ?2 ?1 ?2 ?2 . . . 0 ?′2 ?′m 0 ? SUn = ? ? ... ... ... ... ? 1 ′ 1 1 ?m ?′1 0 ?′m ?′2 . . . 0 ?′m ?′m 0

??

?? ?? ? ? Un,2 ?? ?? ?? ... ?? Un,m

Un,1

?

? ? ? ?. ? ? ?

(6.2.12)

6.2. FEM FOR THE HEAT EQUATION

141

Note that S is just the matrix Aunif that we have already computed in Chapter 1: ? ? 2 ?1 0 0 ... 0 ? ? ? ? ? ?1 2 ?1 0 ... 0 ? ? ? ? ? 0 ?1 2 ?1 . . . 0 ? 1? ?. (6.2.13) S= ? ? h? ? ... ... ... ... ... ... ? ? ? ? ? ? 0 . . . . . . ?1 2 ?1 ? ? ? 0 . . . . . . . . . ?1 2 Similarly, recalling the de?nition for the mass matrix M introduced in (6.2.8), we have that for j = 1, . . . , m

1

MUn =

0

Un ?j .

(6.2.14)

Thus, to compute the mass matrix M one should drop all derivatives from the general form of the matrix for S given by (6.2.13). In other words unlike 1 ′ ′ the form SUn = 0 Un ?j , MUn does not have any derivatives, neither in Un nor in ?j . Consequently ? ? 1 1 1 ?? ?1 ?2 . . . 0 ?1 ?m 0 ? 0 1 1 ? ? 1 ? 1 1 ? 0 ?2 ?1 ? ? . . . ? ? 2 2 2 m ? 0 0 ?. M =? (6.2.15) ? ? ? ? ... ... ... ... ? ? 1 1 1 ? ? ? ? . . . ? ? m 1 m 2 m m 0 0 0

To continue we follow the same procedure as in chapter one recalling that for a uniform partition we have ? ? ? xj ?1 ≤ x ≤ xj ? x ? xj ?1 1? ?j (x) = (6.2.16) xj +1 ? x xj ≤ x ≤ xj +1 h? ? ? ? 0 x∈ / [xj ?1 , xj +1].

142 y 1

CHAPTER 6. THE HEAT EQUATION IN 1D

?j ?1

?j

?j +1

?j +2

xj ?1

xj

xj +1

xj +2

x

Figure 6.2: ?j and ?j +1 .

Thus

1

?j (x)2 dx =

0

1 h2

xj x j ?1

(x ? xj ?1 )2 dx + )3 x j

xj +1 xj

(xj +1 ? x)2

xj +1 xj

=

1 (xj +1 ? x)3 1 (x ? xj ?1 + x j ?1 h2 3 h2 3 3 3 1 h 2 1 h + 2· = h, = 2· h 3 h 3 3

(6.2.17)

and

1 0

1 ?j ?j +1 dx = 2 h = =

xj +1 xj

(xj +1 ? x)(xj ? x) = [P I ]

xj +1 xj

(x ? xj )2 1 ( x ? x ) j +1 h2 2 1 (x ? xj h2 6 )3 xj+1

xj

?

1 h2

xj +1 xj

?

(x ? xj )2 dx 2

1 = h. 6

Obviously we have that

1

?j ?i dx = 0,

0

?|i ? j | > 1.

(6.2.18)

6.3. ERROR ANALYSIS Thus the mass matrix in this case is given by ? 2 1 0 0 6 ? 3 ? 1 2 1 ? 6 0 3 6 ? ? 1 2 1 ? 0 6 3 6 ? M = h? ? ... ... ... ... ? ? 1 ? 0 ... ... 6 ? 0 ... ... ...

143

? ? ... 0 ? ? ? ... 0 ? ?. ? ... ... ? ? 2 1 ? ? 3 6 ?

1 6 2 3

...

0

?

(6.2.19)

6.3

Error analysis

In this section we shall consider a general domain ? with the boundary ? ?. Therefore the analysis are adequate in higher dimensions as well. For our speci?c one dimensional case this means a general interval ? := [a, b] with ? ? = {a, b}. Then a general for of (6.2.5) can be written as ˙ + U ′ v ′ )dxdt = (Uv

In ? In ?

f vdxdt

for all v ∈ Vh ,

(6.3.1)

where Vh = {v (x) : v is continuous, piecewise linear, and v (a) = v (b) = 0}, and in higher dimensional case v (a) = v (b) = 0 is replaced by v |? ? = 0. Note that this variational formulation is valid for the exact solution u and for all v (x, t) such that v (a, t) = v (b, t) = 0: (uv ˙ + u′ v ′ )dxdt =

In ? In ?

f v dxdt,

?v ∈ Vh ,

(6.3.2)

Subtracting (6.3.1) from (6.3.2) we obtain the Galerkin orthogonality relation for the error (ev ˙ + e′ v ′ ) dxdt = 0,

In ?

for all v ∈ Vh .

(6.3.3)

Theorem 26 (A posterirori error estimates). We have the following a posteriori error estimate for the heat conductivity equation given by (6.2.1) e(t) ≤ 2 + ln T max (k + h2 )r (U ) . ε [0,T ] (6.3.4)

144

CHAPTER 6. THE HEAT EQUATION IN 1D

Sketch. To derive error estimates we let ?(x, t) be the solution of the following dual problem ? ? ? ?? ˙ ? ?′′ = 0, ? ? ? = 0, ? ? ? ? ? = e,

in ?

t < T, (6.3.5)

on ? ? t < T, in ? for t = T,

where e = e(t) = e(·, T ) = u(·, T ) ? U (·, T ), T = tN . Note that for w (x, t) = ?(x, T ? t), (t > 0) we can write the backward dual problem (6.3.5) as the following forward problem ? ? ? w ˙ ? w ′′ = 0, ? ? w = 0, ? ? ? ? w = e,

in ?

t > 0, (6.3.6)

on ? ? t > 0, in ? for t = 0.

For this problem we have shown in the energy estimate theorem that

T ε

w ˙ ≤

1 2

ln

T ε

e ,

(6.3.7)

and consequently ( let s = T ? t, then ε → T ? T ? ε → 0, and ds = ?dt) we have for ?:

T ?ε 0

t

s

? ˙ ≤

1 2

ln

T ε

e .

(6.3.8)

Now since ??′′ = ? ˙ we get also

T ?ε 0

?′′ ≤

1 2

ln

T e ε

(6.3.9)

6.3. ERROR ANALYSIS

145

To continue we assume that u0 ∈ Vh then, since (?? ˙ ? ?′′ ) = 0, we can write

T

e(T )

2

=

?

e(T ) · e(T ) dx + e(T ) · e(T ) dx ?

0

?

e(?? ˙ ? ?′′ ) dxdt = [PI in t]

?

=

?

?

e(T ) · e(T ) dx +

e(0) · ?(0) dx

=0

T

+

0 T ?

(e? ˙ + e′ ?′ ) dxdt = {Galerkin Orthogonality (7.1)} e ˙ (? ? v ) + e′ (? ? v )′ dxdt = {PI in x, in 2ed term}

T

=

0 T ?

=

0 T ?

(e ˙ ? e′′ )(? ? v ) dxdt +

0

e′ (? ? v )

=0 T 0 ?

??

dt

=

0 ?

˙ + U ′′ )(? ? v ) dxdt = (f ? U

r (U )(? ? v ) dxdt,

˙ and e′′ = u′′ ? U ′′ to write e ˙ ? U ′′ = where we use e ˙=u ˙ ?U ˙ ? e′′ = u ˙ ? u′′ ? U ′ ˙ ? U := r (U ) which is the residual. Now with mesh variables h = h(x, t) f ?U and k = k (t) in x and t, respectively we can derive an interpolation estimate of the form: ||? ? v ||L2 ≤ k ||? ˙ |L2 + h2 ||?′′ ||L2 ≤ (k + h2 )|? ˙ |L2 + (k + h2 )||?′′ ||L2 , (6.3.10) Summing up we have using maximum principle and the estimates (7.2.2)(7.2.3), basically that

T

e(T )

2

≤

(k + h2 )r (U ) ( ? ˙ + ?′′ )

0 T ?ε 0

≤ max (k + h2 )r (U )

[0,T ]

( ? ˙ + ?′′ ) + 2 max e +2 e .

[T ?ε,T ]

?

≤ max (k + h2 )r (U )

[0,T ]

ln

T ε

This gives our ?nal estimate e(t) ≤ 2 + ln T max (k + h2 )r (U ) . ε [0,T ] (6.3.11)

146

CHAPTER 6. THE HEAT EQUATION IN 1D

?

y

v t

x

The complete proof is given in general form and for higher dimensions in part II. ? Algorithm Starting from the a posteriori estimate of the error e = u ? U for example for ? ? ?u′′ = f, in ? (6.3.12) ? u = 0, on ? ? where r (U ) = |f | + maxIk |[u′]| and [ ] denotes the jump (over the endpoints of a partition interval Ik ), we have the following Algorithm: (1) Choose an arbitrary h = h(x) and a tolerance Tol > 0. (2) Given h, compute the corresponding U . (3) If C hr (U ) ≤ Tol, accept U . Otherwise choose a new (re?ned) h = h(x) and return to step (2) above. ? Higher order elements cG(2), piecewise polynomials of degree 2 is determined by the values of the approximate solution at the end-points of the subintervals. The constructing is through the bases functions of the form: e′ ≤ C h r (U ) , (6.3.13)

i.e.

6.3. ERROR ANALYSIS

147

? Error estimates (a simple case) For ?u′′ = f, 0 < x < 1 associated with Dirichlet (or Neumann) boundary condition we have (u ? U )′ ≤ C h2 D 3 u . (6.3.14) u ? U ≤ C max h h2 D 3 u u ? U ≤ C h2 r (U ) , . (6.3.15) (6.3.16)

where |r (U )| ≤ Ch.

These estimates can be extended to, for example, the space-time discretization of the heat equation. ? The equation of an elastic beam ? ? ? (au′′ )′′ = f, ? ? u(0) = 0, ? ? ? ? u′′ (1) = 0, ? = (0, 1) u′(0) = 0 (Dirichlet) (6.3.17)

(au′′ )′ (1) = 0, (Neumann)

where a is the bending sti?ness, au′ is the moment, f is the load function, and u = u(x) is the vertical de?ection. A variational formulation for this equation can be written as

1 1

au′′ v ′′ dx =

0 0

f vdx,

?v,

such that (0) = v ′ (0) = 0.

(6.3.18)

Here, the piecewise linear ?nite element functions won’t work (inadequate).

148

CHAPTER 6. THE HEAT EQUATION IN 1D

6.4

Exercises

Problem 56. Work out the details with piecewise cubic polynomials having continuous ?rst derivatives: i.e., two degrees of freedom on each node. Hint: A cubic polynomial in (a, b) is uniquely determined by ?(a), ?′ (a), ?(b) and ?′ (b). Problem 57. Prove an a priori and an a posteriori error estimate for a ?nite element method (for example cG(1)) for the problem ?u′′ + αu = f, in I = (0, 1), u(0) = u(1) = 0,

where the coe?cient α = α(x) is a bounded positive function on I , (0 ≤ α(x) ≤ K, x ∈ I ). Problem 58. a) Formulate a cG(1) method for the problem ? ? (a(x)u′ (x))′ = 0, 0 < x < 1, ? a(0)u′ (0) = u , u(1) = 0.

0

and give an a posteriori error estimate.

b) Let u0 = 3 and compute the approximate solution in a) for a uniform partition of I = [0, 1] into 4 intervals and ? ? 1/4, x < 1/2, a(x) = ? 1/2, x > 1/2. c) Show that, with these special choices, the computed solution is equal to the exact one, i.e. the error is equal to 0. Problem 59. Let · denote the L2 (0, 1)-norm. Consider the problem ? ? ?u′′ = f, 0 < x < 1, ? u′ (0) = v , u(1) = 0.

0

a) Show that |u(0)| ≤ u′ and u ≤ u′ . b) Use a) to show that u′ ≤ f + |v0 |.

6.4. EXERCISES

149

Problem 60. Let · denote the L2 (0, 1)-norm. Consider the following heat equation ? ? ? u ˙ ? u′′ = 0, 0 < x < 1, t > 0, ? ? a) Show that the norms: ||u(·, t)|| and ||ux (·, t)|| are non-increasing in time. ||u|| =

1 0

u(0, t) = ux (1, t) = 0, ? ? ? ? u(x, 0) = u (x), 0 < x < 1. 0

1/2

t > 0,

u(x)2 dx

.

b) Show that ||ux (·, t)|| → 0, as t → ∞. Problem 61. Consider the problem ?εu′′ + xu′ + u = f,

c) Give a physical interpretation for a) and b).

in I = (0, 1),

u(0) = u′ (1) = 0,

where ε is a positive constant, and f ∈ L2 (I ). Prove that ||εu′′|| ≤ ||f ||. Problem 62. Give an a priori error estimate for the following problem: (auxx )xx = f, 0 < x < 1, u(0) = u′ (0) = u(1) = u′(1) = 0,

where a(x) > 0 on the interval I = (0, 1). Problem 63. Prove an a priori error estimate for the ?nite element method for the problem ?u′′ (x) + u′(x) = f (x), 0 < x < 1, u(0) = u(1) = 0.

Problem 64. (a) Prove an a priori error estimate for the cG(1) approximation of the boundary value problem ?u′′ + cu′ + u = f where c ≥ 0 is constant. in I = (0, 1), u(0) = u(1) = 0,

(b) For which value of c is the a priori error estimate optimal?

150

CHAPTER 6. THE HEAT EQUATION IN 1D

Problem 65. We modify problem 2 above according to ?εu′′ + c(x)u′ + u = f (x) 0 < x < 1, u(0) = u′ (1) = 0,

where ε is a positive constant, the function c satis?es c(x) ≥ 0, c′ (x) ≤ 0, and f ∈ L2 (I ). Prove that there are positive constants C1 , C2 and C3 such that √ ε||u′ || ≤ C1 ||f ||, ||cu′|| ≤ C2 ||f ||, and ε||u′′|| ≤ C3 ||f ||, where || · || is the L2 (I )-norm. Problem 66. Show that for a continuously di?erentiable function v de?ned on (0, 1) we have that ||v ||2 ≤ v (0)2 + v (1)2 + ||v ′||2 . Hint: Use partial integration for (x ? 1/2) has the derivative 1.

1/2 0

v (x)2 dx and

1 1/2

v (x)2 dx and note that

Chapter 7 The wave equation in 1d

We start with the homogeneous wave equation: Consider the initial-boundary value problem ? ? ? u ¨ ? u′′ = 0, 0<x<1 t>0 (DE ) ? ? u(0, t) = 0, u(1, t) = 0 t>0 (BC ) ? ? ? ? u(x, 0) = u (x), u ˙ (x, 0) = v0 (x), 0 < x < 1 (IC ). 0

(7.0.1)

Below we shall derive the most important property of the wave equation Theorem 27 (Conservation of energy). For the wave equation (7.0.1) we have that 1 1 1 1 ||u ˙ ||2 + ||u′||2 = ||v0 ||2 + ||u′0 ||2 = Constant, 2 2 2 2 where

1

(7.0.2)

||w ||2 = ||w (·, x)||2 =

0

|w (x, t)|2 dx.

(7.0.3)

Proof. We multiply the equation by u ˙ and integrate over I = (0, 1) to get

1 0 1

u ¨ udx ˙ ?

u′′ u u ˙ dx = 0.

0

(7.0.4)

151

152

CHAPTER 7. THE WAVE EQUATION IN 1D

Using partial integration and the boundary data we obtain

1 2 1d u ˙ dx + u ′ (u ˙ )′ dx ? u′ (x, t)u ˙ (x, t) 0 2 dt 0 1 1 2 1d 1d ′ 2 = u ˙ dx + u dx 0 2 dt 0 2 dt 1 1d ||u ˙ ||2 + ||u′||2 = 0. = 2 dt 2 1 1 0

(7.0.5)

Thus, the quantity 1 1 ||u ˙ ||2 + ||u′||2 = Constant, independent of t. 2 2 (7.0.6)

||u ˙ ||2 is the kinetic Therefore the total energy is conserved. We recall that 1 2 1 ′ 2 energy, and 2 ||u || is the potential (elastic) energy. Problem 67. Show that (u ˙ )′ 2 + u′′ 2 = constant, independent of t. Hint: Multiply (DE): u ¨ ? u′′ = 0 by ?(u ˙ )′′ and integrate over I . Alternatively: di?erentiate the equation with respect to x and multiply by u, ˙ .... Problem 68. Derive a total conservation of energy relation using the Robin ?u + u = 0. type boundary condition: ?n

7.1

FEM for the wave equation

solution u(x, t) for the following problem 0<x<1 u′ (1, t) = g (t, ) t>0 t>0 (DE ) (BC ) (7.1.1)

We seek the ?nite element ? ? ? u ¨ ? u′′ = 0, ? ? u(0, t) = 0, ? ? ? ? u(x, 0) = u (x),

0

u ˙ (x, 0) = v0 (x), 0 < x < 1 (IC ).

We let u ˙ = v , and reformulate the problem as a system of PDEs: ? ? u ˙ ?v =0 (Convection) ? v ˙ ? u′′ = 0 (Di?usion).

(7.1.2)

7.1.

FEM FOR THE WAVE EQUATION

153

Remark 15. We can rewrite the above system as w ˙ + Aw = 0 with ? ? ? ? ? ?? ? ? ? u u ˙ a b u 0 ? ? ? = ? ? , (7.1.3) w = ? ? =? w ˙ + Aw = ? ? + ? v v ˙ c d v 0

A

thus we get the following system of equations ? ? au + bv = ?u ˙ ? cu + dv = ?v. ˙

(7.1.4)

? Consequently we have a = 0, b = ?1 and c = ? ?x 2 , d = 0, i.e. ?? ? ? ? ? ? ? 0 ?1 u 0 u ˙ ?? ? = ? ?. ? ?+? 2 ? 0 ? ?x v 0 v ˙ 2 w ˙ A w

Recalling that u ˙ = v and v ˙ = u′′ (7.1.4) can be written as ? ? au + bv = ?v ? cu + dv = ?u′′ .

2

(7.1.5)

(7.1.6)

?The ?nite element discretization procedure For each n we de?ne the piecewise linear approximations as ? ? U (x, t) = U (x)Ψ (t) + U (x)Ψ (t), n?1 n?1 n n 0 < x < 1, t ∈ In , (7.1.7) ? V (x, t) = V (x)Ψ (t) + V (x)Ψ (t),

n?1 n?1 n n

where

? ? U (x) = U (x)? (x) + . . . + U (x)? (x), n n,1 1 n,m m ? V (x) = V n?1 n?1,1 (x)?1 (x) + . . . + Vn?1,m (x)?m (x).

(7.1.8)

154

ψn(t) 1

CHAPTER 7. THE WAVE EQUATION IN 1D

? j (x)

t n-1

tn

t n+1

xj-1

xj

xj+1

Note that since u ˙ ? v = 0, t ∈ In = (tn?1 , tn ] we have

1 In 0 1

u? ˙ dxdt