Math 610 Section 700 - Recitation week 3 week 4 week 6 week 8 - - PowerPoint PPT Presentation

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Math 610 Section 700 - Recitation

TA: Peng Wei Office: Blocker 502E weip@math.tamu.edu

March 6, 2019

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 1 - Getting started with Matlab

• Desktop Basics
• Matrices and Arrays
• Array Indexing
• Calling Functions
• 2-D and 3-D Plots
• Programming and Scripts
• PublishingMATLABCode
• Loop Control Statements

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 2 - finite difference method

Finite difference for 2nd order ODE Example (Dirichlet boundary)

Solve −u′′(x) + u(x) = f (x), u(0) = u(1) = 0, with manufactured solution u = sin(πx). In this case f (x) = (1 + π2) sin(πx). As usual, we are discretizing the domain into uniformly spaced subintervals (xi, xi+1), i = 0, 1, . . . N, where xi = ih, and h = 1/N. First recall the numerical differentiation, u′′(xi) ≈ 1 h2

• u(xi−1) − 2u(xi) + u(xi+1)
• ,

i = 1, · · · , N − 1.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 2 - finite difference method

Finite difference for 2nd order ODE (cont.)

Denote Ui the numerical approximation to u(xi) for i = 0, 1, · · · , N. Notice that due to the Dirichlet boundary condition, U0 = UN = 0 is automatically given. We then have the linear equation system 1 h2 (−✚ ✚ ❃0 U0 + 2U1 − U2) + U1 = F1 := f (x1), 1 h2 (−U1 + 2U2 − U3) + U2 = F2 := f (x2), · · · , 1 h2 (−UN−2 + 2UN−1 − ✚ ✚ ❃0 UN ) + UN−1 = FN−1 := f (xN−1). N − 1 equations corresponds to N − 1 unknowns.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 2 - finite difference method

Finite difference for 2nd order ODE (cont.)

The matrix-vector form of the linear system becomes ( 1 h2 A + I)U = F, with U =      U1 U2 . . . UN−1      , A =          2 −1 . . . −1 2 −1 . . . −1 2 . . . . . . . . . . . . . . . . . . 2 −1 . . . 2          , F =      F1 F2 . . . FN−1      . The matrix A looks familiar! Question: what if we have non-homogeneous boundary condition? For example, u(0) = α, u(1) = β.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 2 - finite difference method

Finite difference for 2nd order ODE (cont.) Example (Neumann boundary)

Solve −u′′(x) + u(x) = f (x), u′(0) = u′(1) = 0. We manufacture a solution u = cos(πx), with right hand side f = (π2 + 1) cos(πx).

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 2 - finite difference method

Choice 1: we approximate the boundary conditions with first order approximation: 0 = u′(0) ≈ 1 h(U1 − U0), 0 = u′(1) ≈ 1 h(UN−1 − UN). So in addition to the N − 1 equations as before, we have two more equations. U0 − U1 = 0, 1 h2 (−U0 + 2U1 − U2) + U1 = F1 := f (x1), 1 h2 (−U1 + 2U2 − U3) + U2 = F2 := f (x2), · · · , 1 h2 (−UN−2 + 2UN−1 − UN) + UN−1 = FN−1 := f (xN−1), UN−1 − UN = 0.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 2 - finite difference method

Choice 2: we approximate the boundary conditions with second

• rder approximation:

0 = u′(0) ≈ 1 2h(−3U0 + 4U1 − U2), 0 = u′(1) ≈ 1 2h(UN−2 − 4UN−1 + 3UN). So we have N + 1 equations. −3U0 + 4U1 − U2 = 0, 1 h2 (−U0 + 2U1 − U2) + U1 = F1 := f (x1), 1 h2 (−U1 + 2U2 − U3) + U2 = F2 := f (x2), · · · , 1 h2 (−UN−2 + 2UN−1 − UN) + UN−1 = FN−1 := f (xN−1), UN−2 − 4UN−1 + 3UN = 0.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 3 - sparse matrix and operations

• A sparse matrix is a matrix with most of its elements zero.

Generally we say a matrix is sparse when at least half if its elements are zeros.

• Examples: the system matrix obtained from finite

difference method (see week 2) and finite element method (see below). In particular, # non-zero elements = O(N).

• Figure. Sparsity pattern of system matrix for a first-order finite

element in two dimensions. Red spots are possible non-zero elements. Source: www.dealii.org/9.0.0/doxygen/deal.II/step_2.html.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 3 - sparse matrix and operations

• In order to store an element in a 2D matrix A, we only need to

know its coordinations and value.

• A natural sparse format is to use a (row id, col id, val)

triplet, and store them in a sorted list first by row id and then by col id. Disadvantage: slow access.

• Compressed Sparse Column (CSC) format. We create 3 arrays -
• ne of floating-point numbers for val, the other two of integers

for row id and col ptr. We traverse the elements of A column-by-column. Advantage: efficient access.

• val contains all non-zero values in matrix A.
• row id stores corresponding row indexes of all entries in

val.

• col ptr has the location in val where a new column
• starts. Conventionally we add nnz + 1 to the end of

row ptr (nnz means number of non-zero elements).

• Compressed Sparse Row (CSR) format is similar to CSC, except

that it traverses the matrix row-by-row, so we store col id and row ptr.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 3 - sparse matrix and operations

An example: compute xTAT with AT =     3 4 1 5 2 6 7     , x =     8 9 10 11     Then the CSC form is (all indexes begin with 1, as in Matlab): val = [1 2 3 4 5 6 7] row id = [2 4 1 1 3 4 4] col ptr = [1 3 4 4 7 8]

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 3 - sparse matrix and operations

The i−th column in AT can be easily obtained from first checking col ptr(i) to col ptr(i + 1) − 1 to obtain the positions, and get:

• values val(col ptr(i)), · · · , val(col ptr(i + 1) − 1),
• row pos. row id(col ptr(i)), · · · , row id(col ptr(i + 1) − 1).

For example, column 4 of A has entries 4 to 7 − 1, therefore, the column has values 4, 5 and 6, positions 1, 3 and 4. val = [1 2 3 4 5 6 7] row id = [2 4 1 1 3 4 4] col ptr = [1 3 4 4 7 8]

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 3 - sparse matrix and operations

• Intuitively, the matrix-vector operation Ax takes the rows of A,

dot product with x. However, in our CSC format, it is expensive to extract rows.

• To bypass the difficulty, we invoke the identity Ax = (xTAT)T.

Now taking rows of A is equivalent to taking columns of AT, which is cheap in CSC format.

• The logic is as follows, extract the i-th column of A, multiply

the values with values in x at the corresponding positions, and sum them up to obtain the i-th value of the result. Example: column 4 val = [1 2 3 4 5 6 7] row id = [2 4 1 1 3 4 4] col ptr = [1 3 4 4 7 8] xT = [8 9 10 11] So the 4-th column of result is 4 ∗ 8 + 5 ∗ 10 + 6 ∗ 11 = 148.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 3 - sparse matrix and operations

The algorithm looks as follows:

function sol = sparse matrix multiplication(val, row id, col ptr, v) numCols = length(col ptr)−1; % number of columns sol = zeros(numCols,1); for i = 1:numCols for j = col ptr(i) : col ptr(i+1)−1 sol(i) = sol(i) + val(j) ∗ v(row id(j)); end end

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 3 - sparse matrix and operations

• A comparison of memory consumption between sparse and

dense matrices. Recall a double precision floating-point value takes 8 bytes memory in Matlab. See an example of N × N matrix in Matlab: simply the identity matrix.

• sparse converts a dense matrix to sparse form.
• nnz returns the number of non-zero elements in a matrix.
• spy visualize the sparsity structure of the matrix.
• whos var lists the variable var’s information in the

workspace, together with its size, bytes, class, etc.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 3 - sparse matrix and operations

Extra readings:

• SuiteSparse, a suite of sparse matrix algorithms developed

by Tim Davis. Matlab x = A\b uses this package.

• Tim Davis also has an online course “direct methods for

sparse linear systems” available at https://youtu.be/1dGRTOwBkQs.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 4 - finite difference method

Finite difference for heat equation

We consider the following 1-dimensional heat equation ∂ ∂t y(x, t) = 1 10 ∂2 ∂x2 y(x, t), x ∈ (0, 1), t ∈ [0, T], with homogeneous Dirichlet boundary condition y(0, t) = y(1, t) = 0, t ∈ [0, T], and initial condition y(x, 0) = y0 = sin(πx). It has a unique solution y(x, t) = e−π2t/10 sin(πx).

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 4 - finite difference method

Finite difference for heat equation (cont.)

Suppose we are discretizing the space domain into uniformly spaced subintervals (xi, xi+1), i = 0, 1, . . . N, where xi = ih, and h = 1/N. Again recall the numerical differentiation, f

′′(x0) ≈ 1

h2

• f (x0 − h) − 2f (x0) + f (x0 + h)
• .

Translated into our PDE, we have the following second derivative midpoint formula ∂2 ∂x2 y(xi, t) ≈ 1 h2

• Yi−1(t) − 2Yi(t) + Yi+1(t)
• ,

i = 1, · · · , N − 1. Notice that due to the Dirichlet boundary condition, Y0(t) = YN(t) = 0 is automatically given.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 4 - finite difference method

Finite difference for heat equation (cont.)

Rewriting the approximated equation into matrix-vector form, d dt Y = AY + F, t ∈ [0, T], Y(0) = Y0, with Y(t) =      Y1(t) Y2(t) . . . YN−1(t)      , A = − 1 10h2          2 −1 . . . −1 2 −1 . . . −1 2 . . . . . . . . . . . . . . . . . . 2 −1 . . . 2          , F =      . . .      . Pay attention to F – what if the b.c. is non-homogeneous? Now the problem becomes an ODE system.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 4 - differential equations

Recall that the Euler’s method for ∂ty = f (y, t) involves approximating the time derivative by d

dt y(tn) ≈ y(tn+1)−y(tn) tn+1−tn

. The different treatment to f (y, t) results in different methods:

• Forward Euler: f (y, t) → f (y(tn), tn),
• Backward Euler: f (y, t) → f (y(tn+1), tn+1),
• Crank-Nicolson: f (y, t) → 1

2(f (y(tn), tn) + f (y(tn+1), tn+1)).

They can be formulated into a general scheme y(tn+1) − y(tn) tn+1 − tn = θf (y(tn), tn) + (1 − θ)f (y(tn+1), tn+1), with θ = 0, 1/2, 1 respectively. Figure out the correspondence! Also recall the stability requirement and convergence rate.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 4 - differential equations

d dt Y = AY + F, t ∈ [0, T], Y(0) = Y0, The general Euler scheme reads Y(tn+1) − Y(tn) k = θAY(tn)+(1−θ)AY(tn+1)+θF(tn)+(1−θ)F(tn+1), with θ = 0, 1/2, 1 respectively. Solving for Y(tn+1) thus is equivalent to solving the linear equation (I − k(1 − θ)A)Y(tn+1) = (I + kθA)Y(tn) + θF(tn) + (1 − θ)F(tn+1). We are already familiar with this linear equation involving sparse matrix.

We are interested in

• the stability condition,
• convergence rate with respect to time, and
• convergence rate with respect to space.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

• assembly of the stiffness matrix
• assembly of the right hand side

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

In one dimensional space, given a domain (a, b), consider a uniform mesh a = x0 < x1 < x2 < . . . < xN−1 < xN = b with the mesh size h = (b − a)/N. For j = 0, . . . , N, the global piecewise linear basis function φj is defined by φj(xj) = 1 and φj(xi) = 0 ∀i = j.

φ0 φ4 φ8

• Figure. Global basis functions.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

The stiffness matrix A = {aij}(N+1)×(N+1), where aij = b

a

a(x)φ′

iφ′ j dx for i, j = 1, . . . , N − 1.

Key property:

• For a fixed index i, for j such that |j − i| > 1, we have

φ′

i(x)φ′ j(x) = 0.

• Thus, for ith row of the matrix A, only ai,i−1, aii, ai,i+1 are

nonzero values.

• This shows that our stiffness matrix is a sparse matrix.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

The stiffness matrix A = {aij}(N+1)×(N+1), where aij = b

a

a(x)φ′

iφ′ j dx for i, j = 1, . . . , N − 1.

Key property:

• For a fixed index i, for j such that |j − i| > 1, we have

φ′

i(x)φ′ j(x) = 0.

• Thus, for ith row of the matrix A, only ai,i−1, aii, ai,i+1 are

nonzero values.

• This shows that our stiffness matrix is a sparse matrix.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

Let us compute the entry aii = b

a a(x)φ′ iφ′ i dx . Since φi is

supported on [xi−1, xi+1] (see for example φ4 in Figure 2), we can write aii = b

a

a(x)φ′

iφ′ i dx =

xi+1

xi−1

a(x)φ′

iφ′ i dx

= xi

xi−1

a(x)φ′

iφ′ i dx +

xi+1

xi

a(x)φ′

iφ′ i dx

We see that the entries aii are composed of integrals on the cells [xi, xi+1] and [xi−1, xi]. Therefore, we loop over all cells, and on each cell, we only focus on the basis functions whose values are non-zero.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

φj-1 φj

• Figure. Local basis functions.

We can write down these two functions restricted in [xj−1, xj] (or the local basis functions), i.e. φj−1(x) = (xj − x) (xj − xj−1) and φj(x) = (x − xj) (xj − xj−1). You can also write down their derivatives, i.e. φ′

j−1(x) = −

1 (xj − xj−1) and φ′

j(x) =

1 (xj − xj−1).

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

φj-1 φj

• Figure. Local basis functions.
• xj

xj−1 a(x)φ′ j−1φ′ j−1 dx is a part of aj−1,j−1.

• xj

xj−1 a(x)φ′ j−1φ′ j dx is a part of aj−1,j.

• xj

xj−1 a(x)φ′ jφ′ j−1 dx is a part of aj,j−1.

• xj

xj−1 a(x)φ′ jφ′ j dx is a part of aj,j.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

local to global

localStiff =      xj

xj−1

a(x)φ′

j−1φ′ j−1 dx

xj

xj−1

a(x)φ′

j−1φ′ j dx

xj

xj−1

a(x)φ′

jφ′ j−1 dx

xj

xj−1

a(x)φ′

jφ′ j dx

    

• A =

     · · · · · · · · · · · · · · · aj−1,j−1 aj−1,j · · · · · · aj,j−1 aj,j · · · · · · · · · · · · · · ·     

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 6 - 1D finite element method

Right hand side assembly

We also compute the global right hand side vector by adding the contributions from the local right hand side vector:      xj

xj−1

f (x)φj−1 dx xj

xj−1

f (x)φj dx      . Here we compute each integral using the two-point Gaussian Quadrature scheme. For example: xj

xj−1

f (x)φx dx ≈

2

• i=1

wif (qi)φ(qi), where w1 = w2 = (xj − xj−1)/2, q1 = xj−1 + xj−xj−1

2

(1 −

1 √ 3) and

q2 = xj−1 + xj−xj−1

2

(1 +

1 √ 3).

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 8 - mesh generation using TRIANGLE

• mesh generation rules and assumptions
• instruction on how to use triangle

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 8 - mesh generation using TRIANGLE

We want to numerically solve the second order elliptic problem in two dimensional space. For example, consider the following weak formulation: find u ∈ H1(Ω) satisfying

• Ω

(∇u·∇v+quv) dx =

• Ω

fv dx+

• ∂Ω

gv ds, for all v ∈ H1(Ω). To begin with, we need to discretize a bounded domain Ω ⊂ R2 into triangles. To simplify our problem, let us assume that our bounded domain Ω is polygonal.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 8 - mesh generation using TRIANGLE

Ω is subdivided into many small closed (including boundary)

• triangles. T denotes the collection of all the triangles.

Since Ω is open (not including boundary), Ω is equal to interior

• f the union of all the triangles, i.e. Ω =
• τ∈T τ
• .
• 1
• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1

• 1
• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1

• Figure. A uniform structured mesh in [−1, 1]2.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 8 - mesh generation using TRIANGLE

Two basic rules:

1 One triangle cannot overlap the other triangles. 2 Given two distinct triangles τi and τj, τi ∩ τj should be

• nly one of the following:

(1) an empty set ∅, (2) a vertex, or (3) a complete edge of both two triangles.

This indicates that the following case is not allowed.

hanging node

• Figure. A non-conforming mesh.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 8 - mesh generation using TRIANGLE

Two assumptions on the mesh:

1 Shape regularity. Cannot be too “flat”. There exists a

constant c so that hτ rτ ≤ c, for all τ ∈ T . Here: hτ = the radius of the smallest ball inscribing τ, rτ = the radius of the largest ball inscribed in τ.

2 Quasi-uniformity. Should be of “similar size”. There

exists a constant ρ so that ρ max

τ∈T hτ ≤ hτ

for all τ ∈ T .

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 8 - mesh generation using TRIANGLE

Source code for triangle: www.cs.cmu.edu/˜quake/triangle.html. In Linux OS or macOS, you can compile the source code in terminal by typing

make all

In mac, you need to install Xcode, command tools and X11. Then modify line 76 of makefile by changing the SWITCHS flag to

CSWITCHES = −O −DNO TIMER −I/usr/X11/include −L/usr/X11/lib

• 1
• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1

• 1
• 0.8
• 0.6
• 0.4
• 0.2

0.2 0.4 0.6 0.8 1

• Figure. A quasi-uniform mesh for a L-shaped domain.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 8 - mesh generation using TRIANGLE

After compiling the code, we get the executable file triangle. Add the poly-file Lshaped.poly to your triangle folder. In terminal, we can generate an L-shaped mesh by typing

./triangle −a.02 −e −q30 Lshaped.poly

• -a.02 controls the maximum area of generated triangles,
• -e outputs (to an .edge file) a list of edges of the triangulation,
• -q30 shows the smallest angle of each triangle ≥30 degree.

Output files include:

• Lshaped.1.node,
• Lshaped.1.ele, and
• Lshaped.1.edge.

Math 610 Section 700 - Recitation week 1 week 2 week 3 week 4 week 6 week 8

Week 8 - mesh generation using TRIANGLE

The MATLAB function readmesh.m will read the files and store the mesh data in MATLAB. The function call is as follows:

[n node, n ele, node, ele, n edge, edge, is edge at bdy]... = readmesh(’Lshaped.1.node’, ’Lshaped.1.ele’, ’Lshaped.1. edge’);

• n node: the number of vertices;
• n ele: the number of elements;
• n edge: the number of edges;
• node: a n node by 2 matrix listing the coordinates of vertices;
• ele and edge: lists of triangle elements and edges;
• is edge at bdy: a boolean vector to check whether each edge

in the edge list is on the boundary.

For example, in the first row of ele, {8, 21, 77} gives the triangle with 8th, 21st and 77th vertices in the node list.