Lecture 10: 29 September 2010 2 in an Lx×Ly box 2d finite-difference discretizations: discretized the 2d Laplacian operator -with M×N points, for Dirichlet (0) boundaries, so that u(mΔx,nΔy) is approximated by MN degrees of freedom um,n. Showed that simply applying the 1d center-difference rule along x and y in which the Laplacian at (m,n) 2" approximation for -5-point stencilresults in a (famous) "depends on u at (m,n) and the 4 nearest-neighbor points. In order to express this as a matrix A, however, we need to "flatten" the 2d grid of points um,n into a single column vector u with MN components. There are multiple ways to do this, but a standard and easy scheme is the "row-major" approach in which u is formed from the contiguous y-columns um,: concatenated in sequence. Given this, we then see how to form A by operating one y-column at a time, using the N×N identity matrix IN and the N×N discrete 1d -Laplacian KN (=DTD, -1 times the 1d laplacian we saw before). The ∂2∂y2 matrix Ay is simply a matrix with KN along the diagonal M times, which differentates each y-column block by block. The ∂2∂x2 matrix is a little tricker, but if we think of operating on whole columns then we see that it is just the KM matrix with the entries "multipled" by N×N identity matrices IN. In order to have a convenient way to express this, we use the Kronecker product notation A⊗B [kron(A,B) in Matlab], which multiplies the entries of A by the matrix B to create a matrix of matrices. In this notation, Ay=IM⊗KN/Δy2 and Ax=KM⊗IN/Δx2. Using this machinery, constructed A for M=N=10 in Matlab. Visualized the pattern of nonzero entries with spy. Solved for the eigenfunctions, and plotted a few; to convert a column vector u back into a 2d matrix, used reshape(u,N,M), and plotted in 3d with the surf command. The first few eigenfunctions can be seen to roughly match the sin(nxπx/Lx) sin(nyπx/Ly) functions we expect from separation of variables. However, M=N=10 is rather coarse, too coarse a discretization to have a really nice (or accurate) picture of the solutions. In order to increase M and N, however, we have a problem. If the problem has P=MN degrees of freedom, we need to store P2 numbers (8P2 bytes) just to store the matrix A, and even just solving Ax=b by Gaussian elimination takes about P3 arithmetic operations. Worked through a few numbers to see that even M=N=100 would have us waiting for 20 minutes and needing a GB of storage, and 3d grids (e.g. 100×100×100) seem completely out of reach. Saving grace, however, is sparsity: the matrix is mostly zero (and in fact the 5-point stencil A has at most 5P nonzero entries). This means that, first, you can store only the nonzero entries, greatly reducing storage. Second, it turns out there are ways to exploit the sparsity to solve Ax=b much more quickly, and there are also quick ways to find a few of the eigenvalues and eigenvectors. In Matlab, you exploit sparsity by using the sparse command and friends to create sparse matrices (in particular, speye creates a sparse identity matrix). Once you have a sparse matrix, Matlab automatically uses algorithms to exploit sparsity if you solve Ax=b by x=A\b and use the eigs function to find a few eigenvalues (instead of eig).Further reading: Section 3.5 of the Strang book on 2d finite differences, section 7.1 on sparsity.MIT OpenCourseWare http://ocw.mit.edu 18.303 Linear Partial Differential Equations: Analysis and Numerics Fall 2010 For information about citing these materials or our Terms of Use, visit:
View Full Document
Unlocking...