# Poisson equation solver

I would like to solve the poisson equation on an annulus (the problem is also rotation symmetric).
Is there a recommended package for this? Ideally I would like something that does not need manual installation of dependencies and can do the job in < 50 lines of code say. I do remember, that there was some basic PDE support in `DiffEq`, but it seems that it is gone?

Have you looked at Approxfun.jl? If you can assume the problem is radially symmetric you just have an ODE BVP, which they show how to solve in the readme for several examples.

2 Likes

Thanks thats a clever idea! I guess I could also use DiffEq for the ODE then. I am still interested in a PDE solution (I want to tweak the problem later and then it is not rotation symmetric anymore).

You can also use the space `Fourier() * Chebyshev(a..b)` for an annulus

1 Like

That is good to know. I figured there was an easy way to use ApproxFun.jl to do this, but havenâ€™t really learned the package well-enough yet to know what it is myselfâ€¦

The examples might help here, thereâ€™s several examples of PDEs with different spaces: https://github.com/JuliaApproximation/ApproxFunExamples

1 Like

@dlfivefifty Can you help me with the boundary conditions? I just want Dirichlet = 0. condition along the boundary of the annulus. So I need no condition in the â€śFourier directionâ€ť and the Dirichlet condition in the â€śChebyshev directionâ€ť, right?. But I am not sure to code this?

For zero Dirichlet the easiest is to incorporate it into the basis, this seems to work well:

``````using ApproxFun

S = JacobiWeight(1,1,Jacobi(1,1,1..2)) * Fourier()
r = Fun((r,Î¸) -> r, Chebyshev(1..2) * Fourier())
f = Fun((r,Î¸) -> exp(r*cos(Î¸)), Jacobi(1,1,1..2) * Fourier())
Dr = Derivative([1,0]) : S
DÎ¸ = Derivative([0,1]) : S;
r2L = r^2 * Dr^2 + r * Dr + DÎ¸^2
u = \(r2L, r^2 * f; tolerance=1E-8)
``````

Though itâ€™s not particularly fast.

2 Likes

Thanks! I am trying to understand this a bit better. I am not very familiar with Jacobi polynomials, is it correct that `JacobiWeight(1,1,Jacobi(1,1,1..2))` is a basis that consists of functions vanishing at the boundary points? And `operator : S` restricts the domain of the operator to a function space `S`?

Pretty much. `JacobiWeight(b, a, space)` on `-1..1` weights the `space` by `(1+x)^b*(1-x)^a`, so if `a = b = 1` then we get vanishing at the boundary. If the interval is other than `-1..1`, everything is translated to `-1..1` first. `Jacobi(1,1,1..2)` is then the space, it could probably be `Chebyshev(1..2)`, but Jacobi polynomials have some special relationships that mean the operators are a bit sparser.

But given how slow the calculation is, this not clearly being leveraged properly. Probably because the operators are falling back to slow `getindex` instead of fast operator building routines. Thereâ€™s a long term goal to make this â€śhigh performanceâ€ť via LazyArrays.jl + InfiniteArrays.jl + BlockBandedMatrices.jl, but thatâ€™s probably a year away from coming together. (This is a really â€śeasyâ€ť problem as there are no corners, so if the generation of the operators was fast enough, it should feel instantaneous.)

1 Like

Weighting a space by `w` means replacing the measure `mu` by `w^(-2) * mu` and each basis element `f_i` by `w*f_i`?

Thereâ€™s no notion of â€śmeasureâ€ť/inner product in ApproxFun, a space is usually just a span of a basis. The second part is correct: the basis is weighted as `w*f_i`.

1 Like

Thanks again for all the detailed answers! And of course for creating `ApproxFun`. It amazes me how easy to use and powerful it is at the same time!

1 Like