Hello,
I’m interested in implementing the laplacian operator \nabla^2=\partial_x^2+\partial_y^2, such that I can apply it to a matrix A_{ij} and obtain the corresponding matrix: B_{ij}=(\nabla^2A)_{i,j}. For example, the result of the action of the operator \nabla^2 on A, is given by a matrix B with elements:
B[i,j]=\frac{A[i+1,j]-2A[i,j]+A[i-1,j]}{dx^2}+\frac{A[i,j+1]-2A[i,j]+A[i,j-1]}{dy^2}
with periodic boudary conditions (of course this is just an approximation valid up to order dx^2+dy^2, and more accurate formulas can be used).
I was thinking of using DiffEqOperators.jl, but it seems that it is in the process of being deprecated. So does anyone have any suggestions on how I should implement it? Performance and efficiency are a big concern since I’m gonna use this in a SDE that I need to solve.
Thanks!
Don’t know how this fits on your efficiency concern as it (may) implies making a grid copy. Sorry for the weird syntax but this (grdmath) GMT.jl module is impossibly hard to convert into the nice keyword names syntax of the rest of the GMT.jl package.
using GMT
G = GMT.peaks();
Gc = gmt("grdmath ? CURV", G); # The Laplacian
You can always just write a loop in Julia. Loops are fast. For example, this post give example code for \nabla^2 in N dimensions with a loop: Seemingly unnecessary allocation within for loop - #8 by stevengj … writing code that only works in 2d is even easier.
For periodic boundary conditions, the naive approach is to compute i\pm 1 modulo the size in the loop body, but this will slow things down a lot. A very flexible approach is to use “ghost cells” (padding) as in the code linked above, which lets the loop body not care about the boundary conditions, which are updated in a separate (cheap) loop over the boundaries.
You may want to have a look at ParallelStencil.jl. Alternatively, I like defining my grid as a graph with Graphs.jl (see the grid function), for which you can easily derive a Laplacian matrix with laplacian_matrix, retrieved as a sparse array.
Thanks for the reply! Yes, using loops was initially my approach. However, I’m solving a system of stochastic differential equations using DifferentialEquations.jl .
The reason why I need the laplacian of a matrix is because my deterministic drift function F involves that operator. To be more specific, my equation looks something like
\begin{equation}
du=f(u,du,t)dt+gdW
\end{equation}
where u is a matrix, and F involves \nabla^2 u.
So I need to implement the laplacian operator inside the f function, and I think vectorized operations are preferred than loops in this case.
Loops can be just as fast as (and often faster than) vectorized operations, regardless of whether they are in your ODE function. An advantage of vectorized operations is that they can be sometimes easier for AD tools to handle, if you are using an implicit method that requires the Jacobian of your rhs.