# Finite difference Laplacian with five-point stencil
**URL:** https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014
**Category:** Numerics
**Tags:** diffeq
**Created:** [June 6, 2019, 4:23pm UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014 "2019-06-06T16:23:04Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@Vortico](https://discourse.julialang.org/u/Vortico)
#### Post date: [June 6, 2019, 4:23pm UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/1 "2019-06-06T16:23:04Z")
This is my first day learning Julia, and I’m not sure of an elegant way to compute the finite difference Laplacian on an NxM Float32 array.
\Delta f(x,y) \approx f(x-1,y) + f(x+1,y) + f(x,y-1) + f(x,y+1) - 4f(x,y)
`Base.diff` is implemented using `view(a, r1...) .- view(a, r0...)`, so I imagine this would look the same. However, after handling the boundaries, this solution seems messier than I’d think it would be.
The problem I’m testing is the 2D heat equation on a Cartesian grid using DifferentialEquations.jl.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 6, 2019, 5:20pm UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/2 "2019-06-06T17:20:10Z")
No need to dig through the standard library in the hopes of finding a built-in “vectorized” routine for this. A loop (in a function) will be clearer and faster.
```julia
function ∇²(f)
∇²f = zero(f) # initialize to zero since we don't touch the boundaries
for y = 2:size(f,2)-1, x = 2:size(f,1)-1 # compute ∇²f in interior of f
∇²f[x,y] = f[x-1,y] + f[x+1,y] + f[x,y-1] + f[x,y+1] - 4f[x,y]
end
return ∇²f
end
```
Note that I only looped over the “interior” of the array `f`. You could write separate loops to handle the boundaries. But a much easier solution in this kind of code is to use “ghost cells”: simply define extra array elements around the boundaries and assign them to whatever boundary condition you want (e.g. Dirichlet, periodic, …) before computing the Laplacian.
If you are performing additional calculations on the Laplacian (e.g. explicit timestepping of a heat equation), you can combine the above loop with your other calculations rather than allocating and returning a temporary `∇²f` array. There are also fancier tricks using [`Base.Cartesian`](https://docs.julialang.org/en/v1.1/devdocs/cartesian/) and `@nloops` to write a Laplacian implementation that works for any number of dimensions.
On the other hand, if you are performing sparse solves, e.g. solving Poisson’s equations, then you typically want the Laplacian operator as an explicit sparse matrix. On a Cartesian grid, this can be accomplished neatly with Kronecker products, as [shown in this notebook from one of my classes](https://nbviewer.jupyter.org/github/mitmath/18335/blob/master/notes/Nested-Dissection.ipynb).
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 6, 2019, 6:21pm UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/3 "2019-06-06T18:21:45Z")
You might want to take a look at this tutorial for optimizing these kinds of calculations: [http://juliadiffeq.org/DiffEqTutorials.jl/html/introduction/optimizing\_diffeq\_code.html](http://juliadiffeq.org/DiffEqTutorials.jl/html/introduction/optimizing_diffeq_code.html) . (And you can optimize it even more as well…)
---
### Author:  [@Vortico](https://discourse.julialang.org/u/Vortico)
#### Post date: [June 7, 2019, 1:24am UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/4 "2019-06-07T01:24:22Z")
Thanks! The sparse matrix method works great too, but I’m working with dense matrices for now.
For a 2-torus, would you recommend ghost cells in `u` or handling the indices in the explicit loop? Is there a way to create an array of size `(size(u,1)+2,size(u,2)+2)` and then create a “view” so that `u_ghost[0]` and `u_ghost[size(u,1)+1]` are valid indices?
---
### Author:  [@chakravala](https://discourse.julialang.org/u/chakravala)
#### Post date: [June 7, 2019, 1:06pm UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/5 "2019-06-07T13:06:49Z")
Using [Grassmann.jl](https://github.com/chakravala/Grassmann.jl) you can use automatic differentiation to compute the fully generalized Hodge De Rahm Laplacian and Dirac operators in any dimension with any metric. It’s a very simple automatic differentiation formula with Grassmann.jl, but I won’t explain the details until my presentation at JuliaCon 2019, which will be supplemented by a paper with detailed proofs.
Sorry if this doesn’t answer your question with finite difference, I’m just a bit excited to share what I’m working on and hope it will be useful once published.
With this new method, you won’t have to worry about what kind of stencil is used, as the algebra is all automatically taken care of.
---
### Author:  [@BeastyBlacksmith](https://discourse.julialang.org/u/BeastyBlacksmith)
#### Post date: [June 7, 2019, 2:18pm UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/6 "2019-06-07T14:18:33Z")
You could use an [Offset Array](https://github.com/JuliaArrays/OffsetArrays.jl) for this.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 7, 2019, 2:25pm UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/7 "2019-06-07T14:25:38Z")
> [@Vortico](#):
>
> For a 2-torus, would you recommend ghost cells in `u` or handling the indices in the explicit loop?
Ghost cells. See also: [Arrays with periodic boundaries - #4 by stevengj](https://discourse.julialang.org/t/arrays-with-periodic-boundaries/4015/4)
(You’ll need an explicit loop anyway to update the ghost cells on each iteration. But the whole approach is so much more flexible because it separates the boundary conditions from the stencil.)
---
### Author:  [@pxshen](https://discourse.julialang.org/u/pxshen)
#### Post date: [August 12, 2022, 12:10am UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/8 "2022-08-12T00:10:52Z")
Two libraries to do this:
1. [EquivariantOperators.jl](https://aced-differentiate.github.io/EquivariantOperators.jl/) scales the stencil coeffs wrt your finite difference cell dimensions, which can in general be noncartesian. Set `border = :circular` for periodic boundaries in the constructor.
```julia
using EquivariantOperators
cell = dx * [1 0; 0 1]
▽2 = Laplacian(cell)
▽2(your_array)
```
1. [Images.jl](https://juliaimages.org) does generic 1,2,1 stencil
```julia
using Images
imgl = imfilter(img, Kernel.Laplacian())
```
---
### Author:  [@SrivatsaPrasad](https://discourse.julialang.org/u/SrivatsaPrasad)
#### Post date: [February 21, 2024, 7:23pm UTC](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/9 "2024-02-21T19:23:31Z")
FYI, the link to Steven’s notebook with the sparse matrix implementation of the Laplacian is broken - seems like [this](https://nbviewer.org/github/mitmath/18S096/blob/409bf1c1cbc8ed0f70afeb0f885ddc382f5138be/lectures/other/Nested-Dissection.ipynb) is the new correct link.