# Solving 2D poisson equation in polar coordinates (finite differences)

**URL:** <https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280>\
**Category:** Numerics\
**Tags:** finitediff, differentialequation, approxfun\
**Created:** [January 9, 2022, 12:47pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280 "2022-01-09T12:47:58Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Aran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aran/32/32584_2.png) [@Aran](https://discourse.julialang.org/u/Aran)\
**Post date:** [January 9, 2022, 12:47pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/1 "2022-01-09T12:47:58Z")

</div>

Hello all,  
I am looking to solve the following poisson equation  
\nabla^2 \psi = -\omega \quad \text{and} \quad \nabla^2 = \frac{\partial ^2}{{\partial \rho^2 }} + \frac{1}{\rho }\frac{\partial}{\partial \rho} + \frac{1}{{\rho ^2 }}\frac{\partial ^2 }{{\partial \phi ^2 }}.

Where \omega is known over a rectangular domain of grid points of N by M .  
I am using finite differences to represent the derivatives of \psi

I am curious whether there are any pre-written Poisson solvers, I have seen a related post below, but didn’t know how to apply any of the recommendations. I have experimented a little with sparse matrices and made a solver on my own, however I would much rather adapt something that is pre-written and that has been verified to be working correctly and efficiently.  
Maybe something from ApproxFun.jl or DifferentialEquations.jl?

> [@How to tackle 2D Poisson equation?](https://discourse.julialang.org/t/how-to-tackle-2d-poisson-equation/48185):
>
> Can anyone point me in the right direction for solving a 2D Poisson equation in a circular region? I’m a little overwhelmed by the number of different Julia packages which a google search returns, and it can be hard to work out what’s current, which packages are abandoned or superseded by others, etc. Ideally I’m looking for something which is Julia all the way through, rather than a wrapper for a third party application. I don’t know much about numerical methods for differential equations, h…

Any help is very much appreciated.  
Many thanks

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 9, 2022, 2:56pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/2 "2022-01-09T14:56:22Z")

</div>

> [@Aran](#):
>
> I am curious whether there are any pre-written Poisson solvers

[https://gridap.github.io/Tutorials/dev/pages/t001\_poisson/](https://gridap.github.io/Tutorials/dev/pages/t001_poisson/)

---

<div class="post-metadata">

**Author:** ![Aran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aran/32/32584_2.png) [@Aran](https://discourse.julialang.org/u/Aran)\
**Post date:** [January 9, 2022, 4:13pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/3 "2022-01-09T16:13:12Z")

</div>

Many thanks @stevengj , I did have a read through the gridap.jl documentation for Poisson solvers, but it was quite tricky for me to understand especially since I haven’t used Finite Elements.

Do you happen to know of any matrix packages (for finite differences) that have already been written for solving Poisson equations?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [January 9, 2022, 4:22pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/4 "2022-01-09T16:22:46Z")

</div>

Handling boundary conditions is the tricky part in finite differences. That is easy with finite elements (such as gridap).

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 9, 2022, 10:06pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/5 "2022-01-09T22:06:49Z")

</div>

[https://github.com/SciML/MethodOfLines.jl](https://github.com/SciML/MethodOfLines.jl)

MethodOfLines.jl should handle this automatically. It’s a little early in its development so just open an issue if you run into any issues, but there are test problems with Poisson’s equation already:

[https://github.com/SciML/MethodOfLines.jl/blob/master/test/pde\_systems/MOL\_NonlinearProblem.jl](https://github.com/SciML/MethodOfLines.jl/blob/master/test/pde_systems/MOL_NonlinearProblem.jl)

so in theory it should be fine.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 10, 2022, 2:51pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/6 "2022-01-10T14:51:49Z")

</div>

It is possible in MultivariateOrthogonalPolynomials.jl, and Poisson was one of the examples in an older version that was built on top of ApproxFun.jl: [https://github.com/JuliaApproximation/MultivariateOrthogonalPolynomials.jl/blob/master/examples/Disk%20PDEs.ipynb](https://github.com/JuliaApproximation/MultivariateOrthogonalPolynomials.jl/blob/master/examples/Disk%20PDEs.ipynb)

However, the code has gone “stale” so will be hard to get working, and assumes one can sample ω anywhere, not just a rectangular grid.

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [January 18, 2022, 6:32pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/7 "2022-01-18T18:32:38Z")

</div>

Hi Aran !  
I implemented the discretization of a 2D poisson equation in polar coordinates with finite differences as an example for a paper on a new Krylov method specialized for nonsymmetric linear systems.  
The code is available [here](https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/main/test/get_div_grad.jl#L140-L176) with an [example](https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/main/test/test_utils.jl#L205-L211).  
 ![poisson](https://global.discourse-cdn.com/julialang/original/3X/0/5/051c875280ab5a068c0d0b25ac927645c6416fda.png)

I based my implementation on the following article: [https://onlinelibrary.wiley.com/doi/10.1002/num.1](https://onlinelibrary.wiley.com/doi/10.1002/num.1)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 18, 2022, 8:17pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/8 "2022-01-18T20:17:26Z")

</div>

> [@amontoison](#):
>
> I implemented the discretization of a 2D poisson equation in polar coordinates with finite differences as an example for a paper on a new Krylov method specialized for nonsymmetric linear systems.

If you don’t get a symmetric matrix with Poisson’s equation, you’re doing it wrong, because -\nabla^2 (with appropriate boundary conditions, e.g. Dirichlet, Neumann, or periodic) is a symmetric linear operator (and is also positive definite or semidefinite). (This is still true even for non-constant coefficients, e.g. if you have -a\nabla \cdot (b\nabla) for scalar functions a(x), b(x) \> 0, as long as you use an inner product weighted by 1/a.) This symmetry (+ definiteness) has rather fundamental qualitative consequences for the physics of problems with \nabla^2, especially diffusion and wave problems, so you really want to preserve it in your discretization.

In particular, in cylindrical coordinates, \nabla^2 is symmetric with a weighted inner product (the r \, dr Jacobian factor of polar integrals), and so you should be able to get a symmetric matrix with a diagonal change of basis if you discretized correctly.

(Galerkin FEM discretizations automatically preserve the symmetry because the corresponding weak form is a symmetric bilinear form.)

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [January 18, 2022, 8:59pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/9 "2022-01-18T20:59:50Z")

</div>

Thanks for your explanations @stevengj !

The resulting linear systems is unsymmetric in my case because the second term of the following PDE is discretized with centered differences (u\_{i+1,j} - u\_{i-1,j} / (2 \Delta r).  
 ![polar](https://global.discourse-cdn.com/julialang/original/3X/c/6/c64a08c6b99d6f667895cfecc2277ee7f4c694ed.png)  
Is it also possible to keep the symmetry with finite differences discretizations ?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 18, 2022, 9:28pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/10 "2022-01-18T21:28:49Z")

</div>

> [@amontoison](#):
>
> The resulting linear systems is unsymmetric in my case because the second term of the following PDE is discretized with centered differences

To make it easier to understand, you’ll want to write the first two terms together in a more symmetrical form:

\frac{\partial^2 u}{\partial r^2} + \frac{1}{r} \frac{\partial u}{\partial r} = \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u}{\partial r} \right)

which is now clearly a symmetric operator under the weighted inner product \langle u, v \rangle = \int u v \, r dr with appropriate boundary conditions. (Understanding the symmetry of the infinite-dimensional operator is really critical to discretizing it properly.)

Now, suppose you discretize the \partial / \partial r first-derivative operator by an (n+1)\times n matrix D which takes in n values \vec{u} on a grid G, plus two points on the boundaries (I’m assuming Dirichlet boundaries for simplicity) and returns the centered differences at n+1 points G' (halfway between the original grid points including the boundaries). After a little algebra, the discretized second-derivative operator above is then the n\times n matrix operation:

\left. \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u}{\partial r} \right) \right|\_G \approx -R^{-1} D^T R' D \vec{u}

where R is a diagonal n \times n matrix of r values on G, R' is a diagonal (n+1) \times (n+1) matrix of r values on the grid G', and -D^T turns out to be the center-difference operator on G' that gives you the approximate derivative on G. But -R^{-1} D^T R' D is clearly similar to a symmetric (and positive-definite!) matrix under a diagonal change of basis \vec{v} = R^{1/2}\vec{u}, which corresponds exaclty to the r Jacobian factor in the inner product for the continuous operator. Equivalently, it’s symmetric under a weighted inner product \langle \vec{u}, \vec{v} \rangle = \vec{u}^T R \vec{v}. (Any Krylov method, such as conjugate gradient, is easily adapted to support modified inner products.)

The above works if your domain is an annulus with Dirichlet boundaries, but if your domain includes the origin, you need to be a little more careful in setting up finite differences because of the coordinate singularity. You typically want to impose Neumann and not Dirichlet boundaries at the origin, and you may want to choose your discretization so that the grid G does not include the origin (so that you don’t divide by zero), or otherwise you can work out the derivative rule at the origin by a limiting procedure.

(This can be combined with a periodic center-difference discretization in θ, which is trivial symmetrical, [using a Kronecker product](https://github.com/mitmath/18303/blob/fall16/lecture-10.pdf), and the Kronecker product of symmetric matrices is symmetric.)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 18, 2022, 10:56pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/11 "2022-01-18T22:56:15Z")

</div>

> [@stevengj](#):
>
> If you don’t get a symmetric matrix with Poisson’s equation, you’re doing it wrong

This also means that Poisson is probably a poor test case for non-symmetric iterative methods — even if you discretize it badly and get a non-symmetric matrix, it is _close_ to being similar to a symmetric matrix (because it is _converging_ to a symmetric operator as you refine the discretization). So its eigenvalues will be close to the real line and its eigenvectors will be close to orthogonal (under the r-weighted inner product), which often has a big impact on Krylov methods.

You may want to test with something that is much farther from Hermitian, e.g. a Helmholtz equation with PML absorbing boundaries. See [this example Julia code](https://nbviewer.org/urls/dl.dropbox.com/s/s7x9kojyioib8ba/Helmholtz2d.ipynb) for such an absorbing-Helmholtz finite-difference solver in 2d via Kronecker products.

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [January 19, 2022, 12:17am UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/12 "2022-01-19T00:17:49Z")

</div>

> I have experimented a little with sparse matrices and made a solver on my own, however, I would much rather adapt something that is pre-written and that has been verified to be working correctly and efficiently

Your problem seems simple enough that actually it could be easier to learn how to verify and optimize your own code than to try to adapt an existing one (which you would still need to verify in order to make sure you don’t have any wrong signs, aren’t inputting data in an incorrect order, etc).

To verify your code, you need to set up a “test” against a known exact solution. If you don’t know any exact solution, you can “manufacture” one by yourself as follows:

1. Pick any known analytical function that satisfies the boundary conditions for your problem
2. Apply the Laplacian operator to it
3. Set that as -\omega
4. Apply your solver with the right hand side you just created.

As for performance, as long as you are using a sparse matrix solver you should be doing fine.

Of course, there exist advanced methods in the literature like spectral solvers, or FEM, but I believe it is also ok to use second-order finite differences as well… as long as you are using moderate values of N,M, or don’t need to call this as a part of an expensive fluid solver, etc.

Regards!

---

<div class="post-metadata">

**Author:** ![Aran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aran/32/32584_2.png) [@Aran](https://discourse.julialang.org/u/Aran)\
**Post date:** [January 20, 2022, 8:32am UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/13 "2022-01-20T08:32:37Z")

</div>

Hi @amontoison, this is absolutely great! Thank you so much for code sharing. This really will really help going forward. I have written a similar sort of sparse matrix structure using code based on [this](https://github.com/Emadmasroor/SimpleNavierStokes.jl/blob/9ded56a19908f5805d67c24c9baa585740d96112/src/functions.jl#L45-L66) cartesian version. Once I have tested version up and running I will post back here. Have you tried the methodoflines approach recommended by Chris in this thread?

---

<div class="post-metadata">

**Author:** ![Aran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aran/32/32584_2.png) [@Aran](https://discourse.julialang.org/u/Aran)\
**Post date:** [January 20, 2022, 8:41am UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/14 "2022-01-20T08:41:44Z")

</div>

This Kronecker product operator is absolutely amazing! Thank you so much for the notes link. Is this a commonly used operator? I have never seen it before and feel like a complete novice now that I have! I feel that my lecturer on tensor products never quite bridged the gap between theory and application.  
Once again, thank you for your kind reply.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [January 20, 2022, 1:10pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/15 "2022-01-20T13:10:23Z")

</div>

In case of separable variables like the Poisson equation on a Cartesian grid with constant (or separable coefficient) you may consider direct solvers based on tensorial decomposition of the Laplacian operator like [https://github.com/triscale-innov/LidJul.jl/blob/master/src/poisson2D\_TT.jl](https://github.com/triscale-innov/LidJul.jl/blob/master/src/poisson2D_TT.jl)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 20, 2022, 1:10pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/16 "2022-01-20T13:10:28Z")

</div>

> [@Aran](#):
>
> This Kronecker product operator is absolutely amazing! Thank you so much for the notes link. Is this a commonly used operator?

Yes, it’s a [standard](https://en.wikipedia.org/wiki/Kronecker_product) way to combine linear operators acting along different dimensions of an n-dimensional grid and write multilinear tensor-product operators as ordinary matrices acting on [flattened (vectorized) column vectors](https://en.wikipedia.org/wiki/Vectorization_(mathematics)), and is incredibly useful for expressing finite-difference equations on tensor-product grids.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 20, 2022, 4:27pm UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/17 "2022-01-20T16:27:25Z")

</div>

> [@LaurentPlagne](#):
>
> In case of separable variables like the Poisson equation on a Cartesian grid with constant (or separable coefficient) you may consider direct solvers based on tensorial decomposition of the Laplacian operator

If you write things in the correct form, you don’t need to limit yourself to separable coefficients to use tensor products. But then you can’t use tensor-product direct solvers (analogous to exploiting separable solutions analytically) and instead need to use generic sparse-direct solvers (which are quite fast in 2d) or more complicated methods (e.g. iterative methods, alternating direction preconditioners, etc).

In particular, you write down a sparse gradient operator G via Kronecker products of 1d finite-difference operators on your Cartesian/tensor-product mesh, then an operator like a \nabla \cdot (b\nabla) with arbitrary variable coefficients a(\vec{x}) and b(\vec{x}) can be written as a matrix -A G^T B G where A and B are diagonal matrices. As a bonus, writing the discretization in this form guarantees that your matrix is symmetric negative-definite (under a diagonal change of basis \sqrt{A}), like the original continuous operator.

See e.g. [this example notebook](https://github.com/mitmath/18303/blob/fall16/min-max-examples.ipynb) (which is for a pre-1.0 version of Julia … an updated version of the code is [here](https://github.com/mitmath/18303/blob/spring19/lecture_notes/minmax.ipynb)) demonstrating a Kronecker-product construction of -\nabla \cdot (c\nabla) for arbitrary non-constant coefficients c.

---

<div class="post-metadata">

**Author:** ![BJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bjr/32/36613_2.png) [@BJR](https://discourse.julialang.org/u/BJR)\
**Post date:** [May 26, 2022, 5:47am UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/18 "2022-05-26T05:47:30Z")

</div>

I am trying to implement this equation with the finite difference method as described. But after reading the explanation of stevemgj, I am confused. I am not a physicist. can you please clarify if this method is correct and use be used as a test case?

---

<div class="post-metadata">

**Author:** ![BJR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bjr/32/36613_2.png) [@BJR](https://discourse.julialang.org/u/BJR)\
**Post date:** [May 31, 2022, 4:55am UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/19 "2022-05-31T04:55:18Z")

</div>

Dear [amontoison](https://discourse.julialang.org/u/amontoison),

In this Julia implementation, what are the values of functions, f, and g that appear in the differential equation? I don’t see them in the code? could please clarify? Further, the code only returns A, b, but does not have any technique to solve the Au= b equation for u?

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [June 7, 2022, 1:47am UTC](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280/20 "2022-06-07T01:47:03Z")

</div>

Sorry for the delayed response @BJR,

`f` and `g` are defined [here](https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/main/test/test_utils.jl#L235-L241).  
I only return `A` and `b` because I use this linear system to test the different methods in Krylov.jl.  
For instance, BiLQ can solve this linear system.

With these functions `f` and `g`, we know that exact solution of the PDE.

 ![poisson](https://global.discourse-cdn.com/julialang/original/3X/1/8/182f63ce11d9426e87d213dce6a89586c7675d0e.png)

[Next page](https://discourse.julialang.org/t/solving-2d-poisson-equation-in-polar-coordinates-finite-differences/74280.md?page=2)