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


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.


1 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

  • $\begingroup$ I have a dense matrix that I will invert, I just need to know how to fill the matrix given the difference equation. $\endgroup$
    – kotozna
    Commented May 25, 2020 at 20:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.