Poisson equation solver

General Usage
,
1

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?

2

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
3

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).

4

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

1 Like
5

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…

6

The examples might help here, there’s several examples of PDEs with different spaces: GitHub - JuliaApproximation/ApproxFunExamples: Examples for using ApproFun.jl

2 Likes
7

@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?

1 Like
8

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.

3 Likes
9

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?

10

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
11

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

12

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
13

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