A two-dimensional finite-difference mesh in Matlab

Here are a few ways to generate the matrix for the system of linear equations in Homework 1 in Matlab.

The graph of the matrix is a k-by-k square mesh, so there are n = k^2 unknowns u(i,j), each related to its four neighbors by an equation of the form

 u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4*u(i,j) = 0 

(at the boundary of the mesh, some neighbors are missing and the right-hand side is nonzero). 

The k^2-by-k^2 coefficient matrix is block tridiagonal, with a k-by-k block corresponding to each "line" of unknowns u(:,j). The k-by-k blocks themselves are either tridiagonal (the blocks on the block diagonal) or identity matrices (the blocks above and below the block diagonal). 

The slick way to generate this matrix is as a Kronecker product. The Kronecker product of matrices X and Y has a block the size of Y for each element of X; it contains every product of an element of X and an element of Y. Say "help kron" to Matlab for an example. 

Our coefficient matrix is the sum of two Kronecker products: 

e = ones(k,1); 
a = spdiags([-e 2*e -e] , [-1 0 1], k, k); 
I = speye (k, k); 
A = kron(a,I) + kron(I,a); 

Here a is a k-by-k tridiagonal matrix, I is a (sparse) k-by-k identity, and A is our k^2-by-k^2 block tridiagonal matrix of coefficients. The code for kron ("type kron") is a compact and slightly tricky example of how to build sparse matrices in Matlab! 

You can also create A with "blockdiags", which is a routine to construct block diagonal matrices in the meshpart toolbox from http://www.cerfacs.fr/algor/Softs/MESHPART/index.html. 

a = blockdiags ([-1 4 -1], -1:1, k, k); 
I = speye (k, k); 
A = blockdiags ([-I a -I], -1:1, k, k); 

The meshpart toolbox also has a routine "grid5" that creates this whole matrix, together with the xy coordinates of the mesh points. So, you could say 

k = 20; 
[A,xy] = grid5(k); 
spy(A); 
gplotg(A,xy); 

to see the matrix and its graph.