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 Math 610 Section 700 - Recitation week 3 week 4 week 6 week 8 TA: Peng Wei Office: Blocker 502E weip@math.tamu.edu March 6, 2019

Math 610 Week 1 - Getting started with Matlab Section 700 - Recitation week 1 week 2 week 3 • Desktop Basics week 4 • Matrices and Arrays week 6 week 8 • Array Indexing • Calling Functions • 2-D and 3-D Plots • Programming and Scripts • PublishingMATLABCode • Loop Control Statements

Math 610 Week 2 - finite difference method Section 700 - Recitation week 1 Finite difference for 2nd order ODE week 2 Example (Dirichlet boundary) week 3 Solve week 4 − u ′′ ( x ) + u ( x ) = f ( x ) , u (0) = u (1) = 0 , week 6 week 8 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 ( x i , x i +1 ) , i = 0 , 1 , . . . N , where x i = ih , and h = 1 / N . First recall the numerical differentiation , u ′′ ( x i ) ≈ 1 � � u ( x i − 1 ) − 2 u ( x i ) + u ( x i +1 ) i = 1 , · · · , N − 1 . , h 2

Math 610 Week 2 - finite difference method Section 700 - Recitation week 1 week 2 Finite difference for 2nd order ODE (cont.) week 3 Denote U i the numerical approximation to u ( x i ) for i = 0 , 1 , · · · , N . week 4 Notice that due to the Dirichlet boundary condition, U 0 = U N = 0 is week 6 automatically given. We then have the linear equation system week 8 ❃ 0 1 ✚ h 2 ( − ✚ U 0 + 2 U 1 − U 2 ) + U 1 = F 1 := f ( x 1 ) , 1 h 2 ( − U 1 + 2 U 2 − U 3 ) + U 2 = F 2 := f ( x 2 ) , · · · , ❃ 0 1 ✚ h 2 ( − U N − 2 + 2 U N − 1 − ✚ U N ) + U N − 1 = F N − 1 := f ( x N − 1 ) . N − 1 equations corresponds to N − 1 unknowns.

Math 610 Week 2 - finite difference method Section 700 - Recitation Finite difference for 2nd order ODE (cont.) week 1 The matrix-vector form of the linear system becomes week 2 week 3 ( 1 week 4 h 2 A + I ) U = F , week 6 week 8 with  − 1  2 0 . . . 0   − 1 − 1   U 1 2 . . . 0 F 1     0 − 1 2 . . . 0 U 2 F 2       U =  , A = , F =  .  .   . . . . .   .  . . . . . . .       . . . . . . .       U N − 1 0 0 0 2 − 1 F N − 1   0 0 0 . . . 2 The matrix A looks familiar! Question: what if we have non-homogeneous boundary condition? For example, u (0) = α , u (1) = β .

Math 610 Week 2 - finite difference method Section 700 - Recitation week 1 week 2 week 3 week 4 Finite difference for 2nd order ODE (cont.) week 6 week 8 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 Week 2 - finite difference method Section 700 - Recitation Choice 1: we approximate the boundary conditions with first order approximation: week 1 0 = u ′ (0) ≈ 1 week 2 h ( U 1 − U 0 ) , week 3 0 = u ′ (1) ≈ 1 week 4 h ( U N − 1 − U N ) . week 6 So in addition to the N − 1 equations as before, we have two more week 8 equations. U 0 − U 1 = 0 , 1 h 2 ( − U 0 + 2 U 1 − U 2 ) + U 1 = F 1 := f ( x 1 ) , 1 h 2 ( − U 1 + 2 U 2 − U 3 ) + U 2 = F 2 := f ( x 2 ) , · · · , 1 h 2 ( − U N − 2 + 2 U N − 1 − U N ) + U N − 1 = F N − 1 := f ( x N − 1 ) , U N − 1 − U N = 0 .

Math 610 Week 2 - finite difference method Section 700 - Recitation Choice 2: we approximate the boundary conditions with second order approximation: week 1 week 2 0 = u ′ (0) ≈ 1 2 h ( − 3 U 0 + 4 U 1 − U 2 ) , week 3 week 4 0 = u ′ (1) ≈ 1 2 h ( U N − 2 − 4 U N − 1 + 3 U N ) . week 6 week 8 So we have N + 1 equations. − 3 U 0 + 4 U 1 − U 2 = 0 , 1 h 2 ( − U 0 + 2 U 1 − U 2 ) + U 1 = F 1 := f ( x 1 ) , 1 h 2 ( − U 1 + 2 U 2 − U 3 ) + U 2 = F 2 := f ( x 2 ) , · · · , 1 h 2 ( − U N − 2 + 2 U N − 1 − U N ) + U N − 1 = F N − 1 := f ( x N − 1 ) , U N − 2 − 4 U N − 1 + 3 U N = 0 .

Math 610 Week 3 - sparse matrix and operations Section 700 - Recitation • A sparse matrix is a matrix with most of its elements zero. week 1 Generally we say a matrix is sparse when at least half if its week 2 elements are zeros. week 3 • Examples: the system matrix obtained from finite week 4 difference method (see week 2) and finite element method week 6 (see below). In particular, # non-zero elements = O ( N ). week 8 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 Week 3 - sparse matrix and operations Section 700 - Recitation • In order to store an element in a 2D matrix A , we only need to week 1 know its coordinations and value. week 2 • A natural sparse format is to use a (row id, col id, val) week 3 triplet, and store them in a sorted list first by row id and then week 4 by col id. Disadvantage: slow access. week 6 • Compressed Sparse Column (CSC) format. We create 3 arrays - week 8 one 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 Week 3 - sparse matrix and operations Section 700 - Recitation week 1 An example: compute x T A T with week 2 week 3     0 3 0 4 0 8 week 4 1 0 0 0 0 9 A T = week 6      , x =     0 0 0 5 0 10 week 8    2 0 0 6 7 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 Week 3 - sparse matrix and operations Section 700 - Recitation week 1 week 2 The i − th column in A T can be easily obtained from first checking week 3 col ptr ( i ) to col ptr ( i + 1) − 1 to obtain the positions, and get: week 4 • values val ( col ptr ( i )) , · · · , val ( col ptr ( i + 1) − 1), week 6 week 8 • 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 Week 3 - sparse matrix and operations Section 700 - Recitation • Intuitively, the matrix-vector operation Ax takes the rows of A , week 1 dot product with x . However, in our CSC format, it is expensive week 2 to extract rows. week 3 • To bypass the difficulty, we invoke the identity Ax = ( x T A T ) T . week 4 Now taking rows of A is equivalent to taking columns of A T , week 6 which is cheap in CSC format. week 8 • 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] x T = [8 9 10 11] So the 4-th column of result is 4 ∗ 8 + 5 ∗ 10 + 6 ∗ 11 = 148.

Math 610 Week 3 - sparse matrix and operations Section 700 - Recitation week 1 week 2 week 3 The algorithm looks as follows: week 4 week 6 function sol = sparse matrix multiplication(val, row id, week 8 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 Week 3 - sparse matrix and operations Section 700 - Recitation week 1 week 2 week 3 • A comparison of memory consumption between sparse and week 4 dense matrices. Recall a double precision floating-point week 6 value takes 8 bytes memory in Matlab. See an example of week 8 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 Week 3 - sparse matrix and operations Section 700 - Recitation week 1 week 2 week 3 week 4 Extra readings: week 6 • SuiteSparse, a suite of sparse matrix algorithms developed week 8 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 .