Solving the Helmholtz equation in polar coordinates

Question

I'm trying to (as a simple starter) find the solution to Helmholtz equation in polar coordinates -- I already know what it is and can derive it by hand, but just want to start here before asking Mathematica to do it in more complex coordinates. Sample derivation from Wikipedia.

My code:

pde2 = D[A[r, θ], {r, 2}] + D[A[r, θ], r]/r  + D[A[r, θ], {θ, 2}] / r^2 + k^2 A[r, θ] == 0

soln = DSolve[{pde2, A[a, θ] == 0}, A[r, θ], {r, θ}]

which I believe corresponds to:

$$ A_{rr} + \frac{A_{r}}{r} + \frac{A_{\theta\theta}}{r^2} + k^2 A = 0 $$ with boundary condition

$$ A(a, \theta) = 0 $$

One thing I'm wondering is that I'm not specifying $r \leq a$, but I can't work out how to do that.

Any help, including links to documentation I've missed or misread will be greatly appreciated.

Answer 1

Since your Dirichlet boundary conditions $A(a,\theta)=0$ is rotationally symmetric, you can solve the problem by separation of variables. Here is a way to do all the formal steps of this method in Mathematica. First I define only the left-hand side of the equation as an operator helmholtz, and then I introduce the separation ansatz to get a new form helmholtz2 on which the separation of variables can be performed.

helmholtz = 
  Function[A, 
   D[A, {r, 2}] + D[A, r]/r + D[A, {θ, 2}]/r^2 + k^2 A];

ansatz = ar[r] aθ[θ];

helmholtz2 = 
 Subtract @@ Simplify[helmholtz[ansatz]/ansatz == 0, r > 0]

$$\frac{\text{a$\theta $}''(\theta )}{\text{a$\theta $}(\theta )}+\frac{r \left(r \text{ar}''(r)+\text{ar}'(r)\right)}{\text{ar}(r)}+k^2 r^2$$

rSolution = 
 DSolve[Select[helmholtz2, FreeQ[#, θ] &] == C[1]^2, ar[r], r]

$$\left\{\left\{\text{ar}(r)\to c_2 J_{c_1}(k r)+c_3 Y_{c_1}(k r)\right\}\right\}$$

θSolution = 
 DSolve[Select[helmholtz2, FreeQ[#, r] &] == -C[1]^2, 
  aθ[θ], θ, GeneratedParameters -> B]

$$\left\{\left\{\text{a$\theta $}(\theta )\to B(2) \sin \left(c_1 \theta \right)+B(1) \cos \left(c_1 \theta \right)\right\}\right\}$$

generalSolution = ansatz /. Flatten[Join[rSolution, θSolution]]

$$\left(B(2) \sin \left(c_1 \theta \right)+B(1) \cos \left(c_1 \theta \right)\right) \left(c_2 J_{c_1}(k r)+c_3 Y_{c_1}(k r)\right)$$

Here the undetermined coefficients are named C[1] for the separation constant, C[2] and C[3] for the radial amplitudes, and B[1], B[2] for the angular function. The latter is displayed a little differently in the pasted output (TeXForm), but I need to do this because the solutions to the two separated equations must have separately named constants. This is what the GeneratedParameters is for.

In the separation step, I used the fact that helmholtz2 was correctly simplified by Mathematica to have only terms dependent on one variable at a time. Then I use Select to obtain the r dependent terms and set them equal to the separation constant C[1]^2, likewise with the angle-dependent term (equated to the negative of the same constant).

The rest is done by DSolve.

Here I confirm that the solution is in fact correct:

FullSimplify[helmholtz[generalSolution] == 0]

(* ==> True *)

Answer 2

When k is not specified, the equation that you have is an eigenvalue equation. Mathematica cannot solve EV equations directly, especially in 2D. You have to separate the variables yourself and obtain the two equations. Then, for each value of a you can use FindRoot to numerically find the EVs as no analytical formula exists for them.

In general, Mathematica will help you once you get to a 1D ordinary differential equation, but to get there you have to separate variables yourself. NDSolve can sometimes help, but I remember having trouble using it in more than 1D.