# Non-separable partial differential equation in polar coordinates

I'm trying to solve the Schroedinger equation in 2D for a system interacting via a dipole potential. This means, in effect, I'm trying to solve the nonlinear PDE

$$-\frac{1}{r} \frac{\partial}{\partial r}\left( r \frac{\partial \psi}{\partial r} \right) - \frac{1}{r^2} \frac{\partial^2 \psi}{\partial \phi^2} + \frac{d^2}{r^3}\left(A - B \cos 2\phi \right)\psi(r,\phi)=E\, \psi(r,\phi)$$

with A, B, d, E constants.

I'm having several difficulties getting Mathematica to understand what I want to do, so hopefully someone can help point me in the right direction. Firstly I would like to constrain $\phi$ to be an angular variable, i.e. $\psi(r,\phi)=\psi(r,\phi+2\pi)$: but setting this as one of the boundary conditions in NDSolve is giving lots of errors. I would also like to be able to set at least two small-radius boundary conditions, $\psi(c,\phi)=1$ and $\partial \psi / \partial r |_{r=c}=1$ where c is some small number: when I attempt this Mathematica complains I have overdetermined the system.

Does anyone have any suggestions as to how to proceed?

ClearAll["Global*"]; c =1/1000; d = 1; energy = 0.5; th = \[Pi]/2;
t = NDSolve[{-(1/r) D[f[r, phi], {r, 1}] - D[f[r, phi], {r, 2}]
- 1/r^2 D[f[r, phi], {phi, 2}] + d^2/r^3 (LegendreP[2, Cos[th]] -  3/2 Sin[th]^2 Cos[2 phi]) f[r, phi]
== energy f[r, phi],
(f[r, phi] /. r -> c) == 1,
(D[f[r, phi], {r, 1}] /. r -> c) ==1,
(f[r,phi]) == (f[r, phi + 2 \[Pi]])}, f, {r, c, 10}, {phi, 0, 2 \[Pi]}];
u[r_, phi_] := Evaluate[f[r, phi] /. t];
norm = NIntegrate[u[r, phi], {r, c, 10}, {phi, 0, 2 \[Pi]}];
ParametricPlot3D[{r Cos[phi],r Sin[phi], (u[r, phi])/norm}, {r, c,10}, {phi, 0, 2 \[Pi]}, AxesLabel -> {x, y}, ImageSize -> Large]


Edit: Periodicity seems to drop out of just setting continuity at the boundaries; however, I'm now running into a problem where I get different qualitatively results if I use different MaxStepSizes in NDSolve (compare MaxStepSize=1 and 0.1). So NDSolve is still not happy.

• If I replace (f[r,phi]) == (f[r, phi + 2 \[Pi]]) with (f[r,0]) == (f[r, 2 \[Pi]]), I get an interpolating function, but with a large error warning message. Nov 18, 2014 at 7:49
• @MariusLadegårdMeyer I have also tried this, but I'm not entirely convinced by it: it will give continuity of $\psi$ at the boundary, but won't force it to be a function periodic over all $\phi$, which is what the real solution should be. Nov 18, 2014 at 9:07
• Not sure either, but I plotted the solution for many values of r as function of phi and they are all periodic. However, the solution blows up at phi = \[Pi], so that's the reason for the error message. Nov 18, 2014 at 10:18
• Hmm, you're right, it does look periodic, so perhaps that isn't too much of a problem. The blowing up, on the other hand, is still not appreciated ... although off the top of my head I'm not sure what can be done to deal with that. Nov 18, 2014 at 11:27
• It seems that NDSolve is sometimes a bit picky when it comes to elliptic equations. Just a quick query: are you sure there is actually a bound state solution for these parameters? I've run Eigensystem on the discretised Schrödinger operator and can't seem to find any sensible solution. Nov 19, 2014 at 11:59

I suppose the issue here is that stationary-state Schrödinger should really be treated as an eigenvalue problem. The value E = 0.5 might or might not correspond to a solution.

A way to do this is to use finite differences to discretize the LHS operator into a matrix then diagonalize it. You could do it in polar coordinates as written, using periodic interpolation for φ derivatives. Or you could do it in Cartesian coordinates, which I reckon is rather simpler here:

d = 1;
A = LegendreP[2, Cos[th]];
B = 3/2 Sin[th]^2;
th = Pi/2;

eps = 1/255;
grid = Table[i, {i, -20, 20, 40*eps}];

identityMatrix = SparseArray[{i_, i_} :> 1, Length[grid]];
d2 = NDSolveFiniteDifferenceDerivative[Derivative[2], grid, DifferenceOrder -> 4]["DifferentiationMatrix"];

(* tensor product grid *)
d2x = KroneckerProduct[d2, identityMatrix];
d2y = KroneckerProduct[identityMatrix, d2];
gridT = Flatten[Outer[List, grid, grid], 1];
matrixT[f_] := SparseArray[{i_, i_} :> f@@gridT[[i]], Length[gridT]]

schrodingerMatrix = - (d2x + d2y) + matrixT[Function[{x, y}, (d^2 ((A - B) x^2 + (A + B) y^2))/(x^2 + y^2)^(5/2)]];

(* Get the 200 smallest eigenvalues *)
eigs = Eigensystem[N@schrodingerMatrix, -200];

(* Pick out only real, positive eigenvalues: everything else is spurious *)
sols = {};
For[i = 1, i <= Length[eigs[[1]]], i = i + 1,
If[Abs@Im@eigs[[1, i]] < 10^-8 && Re@eigs[[1, i]] > 0,
sols = Append[sols, {eigs[[1, i]], eigs[[2, i]]}];
];
];
Clear[eigs]


This gives approximate eigenvalues

sols[[;;,1]] = {0.759627, 0.713689, 0.651773, 0.623434, 0.580288, 0.496437, \
0.473723, 0.446067, 0.291805, 0.221632, 0.200286, 0.180741, 0.173824, \
0.158172, 0.105491, 0.079869, 0.0166472, 0.0017899, 0.000675304, \
0.000439906, 0.000117715, 0.0000151839, 4.48709*10^-6}


I personally wouldn't place a great deal of trust in these values as I haven't done proper checks. The corresponding eigenfunctions are probably, qualitatively OK, but again don't take my words for it.

The E ≈ 0 state:

The E ≈ 0.496 state:

• Thanks for your answer Saran. I have my suspicions about the eigenfunctions there, even qualitatively: the system has a reflection symmetry about the x-axis, which isn't present in those eigenfunctions. With regards to the equation as an eigenvalue equation you are of course correct; however unless I'm very mistaken it should be a continuous spectrum of eigenvalues (corresponding sort of to the incident energy of the particle scattering off the potential). Presumably discretizing the LHS is what has made this into a discrete spectrum? Nov 21, 2014 at 11:06