Example for sparse matrix A: problem on n by n grid

Contents

You need to download an m-file

Download the following m-file and put it in the same directory with your other m-files:

Solve problem on n by n grid

We have 2500 unknowns $u_1,\ldots,u_{2500}$ arranged on a 50 by 50 grid:

For each unknown $u$ we have neighbors $u_W, u_E, u_S, u_N$ to the West, East, South, North and an equation

$$a_D u + a_W u_W + a_E u_E + a_S u_S + a_N u_N = \cdots$$

(if some neighbor doesn't exist the corresponding term is omitted)

For our example we choose aD=40, aW=-10.8, aE=-9.2, aS=-9.6, aN=-10.4, and the right side vector b where all entries are 1.

n = 50;                    % unknowns are on n by n grid
M = n^2;                   % number of unknowns
A = spgridmat(n,40,-10.8,-9.2,-9.6,-10.4); % M by M matrix in sparse format
b = ones(M,1);             % right hand side vector

u = A\b;                   % solve linear system
                           % (uses clever row and column permutations to minimize "fill-in")

contour(reshape(u,n,n)',30); axis equal; colorbar % plot solution as contour graph

Fill-In for standard pivoting strategy

Note the "fill-in" in the matrix U. The matrix A has only 12300 nonzero entries, but U has 125049 nonzero entries.

spy(A);
title('nonzero entries in A'); snapnow
[L,U,p] = lu(A,'vector');  % column pivoting, choosing candidate with largest absolute value
spy(U);
title('nonzero entries in U'); snapnow

Fill-In for optimized strategy with row and column pivoting

Now the nonzeros in U have a weird shape. But the number of nonzero entries is only 35913, compared to 125049 for the standard pivoting strategy.

[L,U,p,q] = lu(A,'vector');  % Note: this only works for A in sparse format
spy(U);
title('nonzero entries in U'); snapnow


Published with MATLAB® R2018a