- Home
- Documents
*Sparse Linear Systems and parallel iterative methods Lesson 8. Sparse Linear Systems and parallel...*

prev

next

out of 69

View

0Download

0

Embed Size (px)

DISIM - Università dell’Aquila

Sparse Linear Systems and parallel iterative methods Lesson 8.

Adriano FESTA

Dipartimento di Ingegneria e Scienze dell’Informazione e Matematica, L’Aquila

DISIM, L’Aquila, 07.05.2019

adriano.festa@univaq.it

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Class outline

Poisson Equation

Tridiagonal Systems

General banded matrices

Domain Decomposition Methods

Heat Equation

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Motivation: Poisson Equation

As a typical example of an elliptic partial differential equation we consider the Poisson equation with Dirichlet boundary conditions.

This equation is often called the model problem since its structure is simple but the numerical solution is very similar to many other more complicated partial differential equations,

The two-dimensional Poisson equation has the form

−∆u(x , y) = f (x , y) for all (x , y) ∈ Ω

with domain Ω ∈ R2. The function u : R2 → R is the unknown solution function and the function f : R2 → R is the right-hand side, which is continuous in Ω and its boundary. The operator ∆ is the two-dimensional Laplace operator

∆ = ∂2

∂x2 +

∂2

∂y2

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Motivation: Poisson Equation

Using this notation, the Poisson equation can also be written as

− ∂2

∂x2 u(x , y)−

∂2

∂y2 u(x , y) = f (x , y) for all (x , y) ∈ Ω.

A model problem uses the unit square Ω = (0, 1)× (0, 1)and assumes a Dirichlet boundary condition

u(x , y) = φ(x , y)for all (x , y) ∈ ∂Ω

where φ is a given function and ∂Ω is the boundary of domain Ω, which is ∂Ω = {(x , y)|0 ≤ x ≤ 1, y = 0 or y = 1} ∪ {(x , y)|0 ≤ y ≤ 1, x = 0 or x = 1}. The boundary condition uniquely determines the solution u of the model problem. An example of the Poisson equation from electrostatics is the equation

∆u = − ρ

ε0

where ρ is the charge density, ε0 is a constant, and u is the electrical potential created by the charge.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: Discretization

For the numerical solution of equation −∆u(x , y) = f (x , y), the method of finite differences can be used, which is based on a discretization of the domain Ω ∪ ∂Ω The discretization is given by a regular mesh with N + 2 mesh points in x-direction and in y-direction, where N points are in the inner part and 2 points are on the boundary. The distance between points in the x- or y-direction is h = 1/(N + 1). The mesh points are

(xi , yj ) = (ih, jh) for i, j = 0, 1, ...,N + 1.

The points on the boundary are the points with x0 = 0, y0 = 0, xN+1 = 1, or yN+1 = 1.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: Discretization

The unknown solution function u is determined at the points (xi , yj ) of this mesh, which means that values uij := u(xi , yj ) for i, j = 0, 1, ...,N + 1 are to be found.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: FD Discretization

For the inner part of the mesh, these values are determined by solving a linear equation system with N2 equations.

For each mesh point (xi , yj ), i, j = 1, ...,N, a Taylor expansion is used for the x or y-direction.

The Taylor expansion in x-direction is

u(xi + h, yj ) = u(xi , yj ) + h · ux (xi , yj ) + h2

2 uxx (xi , yj ) +

h3

6 uxxx (xi , yj ) + O(h

4),

u(xi − h, yj ) = u(xi , yj )− h · ux (xi , yj ) + h2

2 uxx (xi , yj )−

h3

6 uxxx (xi , yj ) + O(h

4).

where ux denotes the partial derivative in x -direction (i.e., ux = ∂u/∂x ) and uxx denotes the second partial derivative in x-direction (i.e., uxx = ∂2u/∂x2).

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: FD Discretization

Adding these two Taylor expansions results in

u(xi + h, yj ) + u(xi − h, yj ) = 2u(xi , yj ) + h2uxx (xi , yj ) + O(h4).

Analogously, the Taylor expansion for the y-direction can be used to get

u(xi , yj + h) + u(xi , yj − h) = 2u(xi , yj ) + h2uyy (xi , yj ) + O(h4).

From the last two equations, an approximation for the Laplace operator ∆u = uxx + uyy at the mesh points can be derived

∆u(xi , yj ) = − 1 h2

(4uij − ui+1,j − ui−1,j − ui,j+1 − ui,j−1),

where the higher order terms O(h4) are neglected.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: Discretization

This pattern is known as five-point stencil. Using the approximation of ∆u and the notation fij := f (xi , yj ) for the values of the right-hand side, the discretized Poisson equation or five-point formula results:

1 h2

(4uij − ui+1,j − ui−1,j − ui,j+1 − ui,j−1) = fij for 1 ≤ i, j ≤ N.

For the points on the boundary, the values of uij result from the boundary condition and are given by

uij = φ(xi , yj ) for (i, j) ∈ ∂Ω.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: linear system

The five-point formula including boundary values represents a linear equation system with N2 equations, N2 unknown values, and a coefficient matrix A ∈ RN2×N2 . In order to write the equation system with boundary values in matrix form Az = d, the N2 unknowns uij , i, j = 1, ...,N, are arranged in row-oriented order in a one-dimensional vector z of size n = N2 which has the form

z = (u11,u21, ...,uN1,u12,u22, ...,uN2, .....,u1N ,u2N , ...,uNN).

The mapping of values uij to vector elements zk is

zk := uij with k = i + (j − 1)N for i, j = 1, ...,N.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: linear system

Using the vector z, the five-point formula has the form

1 h2

(4zi+(j−1)N − zi+(j−1)N−1 − zi+(j−1)N+1 − zi+jN − zi+(j−2)N) = di+(j−1)N

for 1 ≤ i, j ≤ N. with di+(j−1)N = fij

one-dimensional vector resulting in the corresponding mapping of f .

Replacing the indices by k = i + (j − 1)N we obtain the easier form

1 h2

(4zk − zk−1 − zk+1 − zk+N − zk−N) = dk

for 1 ≤ k ≤ N2. Thus, the entries in row k of the coefficient matrix contain five entries which are akk = 4 and ak,k+1 = ak,k−1 = ak,k+N = ak,k−N = −1.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

The building of the vector d and the coefficient matrix A = (aij ), i, j = 1, ...,N2 , can be performed by the following algorithm. The loops over i and j , i, j = 1, ...,N, visit the mesh points (i, j) and build one row of the matrix A of size N2 × N2. When (i, j)is an inner point of the mesh, i.e., i, j 6= 1,N, the corresponding row of A contains five elements at the position k, k + 1, k − 1, k + N, k − N for k = i + (j − 1)N. When (i, j) is at the boundary of the inner part, i.e., i = 1, j = 1, i = N, or j = N, the boundary values for φ are used.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: linear system

The linear equation system resulting from this algorithm has the structure

1 h2

B −I 0

−I B . . .

. . . . . . −I

0 −I B

· z = d

where I denotes the N × N unit matrix, which has the value 1 in the diagonal elements and the value 0 in all other entries. The matrix B has the structure

B =

4 −1 0

−1 4 . . .

. . . . . . −1

0 −1 4

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Poisson Equation: sparse structure

In summary, this formula represent a linear equation system with a sparse coefficient matrix, which has non-zero elements in the main diagonal and its direct neighbours as well as in the diagonals in distance N . Thus, the linear equation system resulting from the Poisson equation has a banded structure, which should be exploited when solving the system.

http://adrianofesta.altervista.org/ A. FESTA, Sparse matrices

DISIM - Università dell’Aquila

Tridiagonal Systems

For the solution of a linear equation system Ax = y with a banded or tridiagonal coefficient matrix A ∈ Rn×n, specific solution methods can exploit the sparse matrix structure.

A matrix A = (aij )i,j=1,...,n ∈ Rn×n is called banded when its structure takes the form of a band of non-zero eleme