How to fill matrix entries for two-dimensional implicit finite-difference for the general case

Question

If I have derived a finite-difference formula for a 2D problem, for example something like:

$af_{i,j}+bf_{i-1,j}+cf_{i,j-1}+df_{i-1,j-1}=g_{i,j}$

where f is the unknown function on a grid and everything else is known, how do I fill the entries in the finite difference matrix $A_{i,j}^{i',j'}$ in the general case?

Can I construct some rules based on the above of the form "if $i>j-1$ then..."?

In every book on the subject they simply show the matrix and demonstrate it works. But I'm interested in the general case because I would like to derive some complicated finite difference formulae (much more complicated than, say, 2D diffusion). It might be better to use different notation to the one I used above.

Answer 1

It really depends on how the matrix will be used. In implicit schemes, you typically solve a system of the form $(I-\gamma A)x=b$, where $\gamma$ is some small number related to a time step. For anything much larger than 1-D and small 2-D problems, you have to think very hard about how to actually store the matrix and solve the problem. Two popular methods are

1. Sparse direct methods: MATLAB is particularly good at these and the backslash command will default to these if your matrix is stored as a sparse array. This requires one to actually build the matrix, although you do this by only specifying the nonzero values of the matrix to save space. In multiple dimensions, probably the easist way to do this is to use the Kronecker product to tensor multiple 1-D stencils together. For instance, if $D_1$ is a 1-D finite difference Laplacian and $D_2$ is the 2-D finite difference Laplacian, then one has the relation $$ D_2 = D_1\otimes I+I\otimes D_1. $$ This can be used to drastically simplify the process of building these matrices. One can also then add things like boundary conditions after the differential operaters have been tensored.

2. Krylov iterative methods: With Kyrlov methods, you compute the solution to a system without ever using the matrix itself, i.e., only using matrix-vector products $(I-\gamma A)v$. Since these matrices are so sparse, the matrix-vector product is relatively cheap and as a bonus, one never really needs to worry about actually forming the matrix if you don't want to. For example, you could write a function that takes in a vector $v$ of unknowns, forms it into a 2-D array (in 2 spatial dimensions), then loop over each index and apply the finite differen formula and reshape the result back into a vector and this would be perfectly fine for a Krylov method. Some popular Krylov methods are GMRES for general problems and Conjugate Gradient for spd problems