# 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. – Marius Ladegård Meyer Nov 18 '14 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. – AnotherShruggingPhysicist Nov 18 '14 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. – Marius Ladegård Meyer Nov 18 '14 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. – AnotherShruggingPhysicist Nov 18 '14 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. – Saran Nov 19 '14 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, 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[]], 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? – AnotherShruggingPhysicist Nov 21 '14 at 11:06