2d finite differences implementation for a Numerical Analysis lecture

Hi! I’m new to Julia world and would need some help in the transition from Matlab.
I am teaching a numerical analysis class where I use finite differences on a two dimensional (cartesian but not necessarily square) domain. The Matlab function delsq was very useful as I could then generate the grid by meshgrid, build the FD matrix from it, change boundary conditions to Neumann or Robin to a part of the boundary and also localise the source term, diffusion coeff etc…
I all want are simple teaching codes where I can illustrate all the basic concepts on finite differences in a genuinely 2d implementation (no kronecker of 1d!) - and not black box like in [https://github.com/SciML/DifferentialEquations.jl]
Has someone already tried that?

Maybe the codes in https://github.com/luraess/parallel-gpu-workshop-JuliaCon21/tree/main/scripts can be of use, e.g. here diffusion with a forward euler step https://github.com/luraess/parallel-gpu-workshop-JuliaCon21/blob/main/scripts/diffusion_2D_expl.jl.

3 Likes

Thanks a lot! This is very useful.

Check out this EquivariantOperators.jl tutorial and package website

See, for example, this post, which provides a simple N-dimensional implementation of a finite-difference method for a wave equation using a leapfrog explicit timestepping scheme. It would be pretty easy to change it to the diffusion equation, too.

Among other things, it illustrates computing the finite-difference Laplacian in N dimensions cleanly and efficiently using Julia’s Cartesian-indexing code, and the useful technique of “ghost cells” for boundary conditions (which allow you to separate the boundary-condition code from the finite-difference code).

(Whoops, sorry, I see that this is an old thread that was recently revived.)