I am trying to solve a hybrid Dirac/Schrodinger problem in 2D with a Coulomb potential. I am aware of the excellent post on the 2D Coulomb problem, but my problem lacks cylindrical symmetry so I can't turn the problem into a 1D radial problem. It will have to be 2D. Before tackling the real problem I thought to check my routine by solving the 2D Coulomb problem for either the Dirac or the Schrodinger case.

Schrodinger Case:

Coulomb2DSchrod={-D[u[x,y],x,x]-D[u[x,y],y,y] - (x^2+y^2)^-1 u[x,y]};
C2DS=Simplify[Coulomb2DSchrod/.{u->(fxn[ArcTan[#1], ArcTan[#2]]&)} /.

omega=Rectangle[{-Pi/2, -Pi/2)}, {Pi/2, Pi/2}];

NDEigenvalues[{C2DS, DirichletCondition[fxn[s,t]==0, Abs[s]==Pi/2||Abs[t]==Pi/2]}, 
  fxn, {s,t}∈ omega, 5]

Out:  {0.000606273, 0.00930986, 0.00930986, 0.0432714, 0.0521204}

which is of course, horribly wrong for the 2D problem. Increasing the resolution via:

Method -> {"SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" ->  
{"MaxCellMeasure" -> 0.0001}}}, "Eigensystem" -> {"Arnoldi", MaxIterations -> 40000}}

doesn't help at all. Mathematica complains about the pole at the origin, and gives zeros (basically) for all of the eigenvalues. In the posting I cite above (where the author changed to radial coordinates), Mathematica was clever enough to handle the pole, but here it can't.

Dirac Case (with mass to open gap):

Here I can get the code to run. However, the result of:

Coulomb2DDirac = 
{-I D[u2[x, y], x] - D[u2[x, y], y] - g/Sqrt[x^2 + y^2] u1[x, y] + mass u1[x, y], 
-I D[u1[x, y], x] +  D[u1[x, y], y] - g/Sqrt[x^2 + y^2] u2[x, y] - mass u2[x, y] };

and make a similar transformation to tangent coordinates, and exebute:

With[{max = 1, mass = 1, g = .2}, {ev, ef} = 
 NDEigensystem[{transEqn[Δ, δ, g],
  DirichletCondition[fxn1[s, t] == 0, 
   Abs[s] == Pi/2 || Abs[t] == Pi/2], 
  DirichletCondition[fxn2[s, t] == 0, 
 Abs[s] == Pi/2 || Abs[t] == Pi/2]}, {fxn1[s, t], 
fxn2[s, t]}, {s, t} ∈ omega, max, 
Method -> {"SpatialDiscretization" -> {"FiniteElement", 
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}, 
 "Eigensystem" -> {"Arnoldi", MaxIterations -> 40000}}]]

work at high resolution, the absolute value of the answer as a function of x is jagged:

Absolute value of solution plotted in x-direction

Both of the above problems are well posed and have finite answers. Is there a way to get Mathematica to handle either of these problems cleanly?

BTW, my apologies for poor formatting. It took me over an hour to enter this and it still looks pretty bad.

  • 1
    $\begingroup$ I've run into the sort of jaggedness seen in your second case when using FEM in the past. I suspect it's due to the behavior of InterpolatingFunction (the cusps are the vertices of your finite elements, I think); but I've never been quite clear on how to get around it. I'll be interested to see whether there are any good ways to get around this. $\endgroup$ Commented Aug 8, 2017 at 16:15
  • 1
    $\begingroup$ When I copy and past the code it does not run - so there is not much one can do. The first one has syntax errors and the second does not have all the code. $\endgroup$
    – user21
    Commented Aug 18, 2017 at 23:44