5
$\begingroup$

I am trying to solve 3D Schrodinger equation (with mass equals to 1) and with no potential (I will add a potential later). My goal is to set boundary conditions on the outer boundary of a sphere and see that the result obtained by NDSolve is consistent with it.

The analytical solution for this equation is simply the spherical harmonics times the spherical Bessel function.

So first, I defined:

rCoor[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2]
θCoor[x_, y_, z_] := ArcTan[z, Sqrt[x^2 + y^2]]    
ϕCoor[x_, y_, z_] := ArcTan[x, y]
YlmCart[l_, m_, x_, y_, z_] := SphericalHarmonicY[l, m, θCoor[x, y, z], ϕCoor[x, y, z]]

Now I'm trying to solve the equation, with l=m=0 (the spherical harmonic paramters) and for energy parameter of Ef = 1:

Needs["NDSolve`FEM`"];
l = 0;
m = 0;
Ef = 1;
kf = Sqrt[2 Ef];
rmax = 5 kf^-1;
mesh = ToElementMesh[Ball[{0, 0, 0}, rmax]];

boundary[l_, m_, k_, x_, y_, z_] := 
 SphericalBesselJ[l, k rCoor[x, y, z]] YlmCart[l, m, x, y, z]
op = -(1/2) Laplacian[u[x, y, z], {x, y, z}] - Ef u[x, y, z];
sol = NDSolveValue[{op == 0,DirichletCondition[u[x, y, z] == 
    boundary[l, m, kf, x, y, z] , True]}, u, {x, y, z} ∈ mesh];

Now I'm plotting the result to compare with the analytical solution (checking only in the x direction, in which theta and phi are constant):

Plot[{Re[ sol[x, 0, 0]], 
  Re[boundary[l, m, kf, x, 0, 0]]}, {x, 0, rmax}]

enter image description here

as you can see, the result is just fine and perfectly aligned with the analytical result. Now, when I'm increasing the range of computation by changing rmax, the result start to deviate. This is the result for rmax = 30 kf^-1:

enter image description here

What is going on here? How can I solve this problem?

Improve this question
$\endgroup$
7
$\begingroup$

For rmax = 5 kf^-1, ToElementMesh discretizes the region of integration by default into about 7000 elements. For rmax = 30 kf^-1, it does the same. Thus, resolution is much reduced, and this is the source of the inaccuracy for the larger rmax. To increase resolution and display the mesh, use

mesh = ToElementMesh[Ball[{0, 0, 0}, rmax], MaxCellMeasure -> {"Length" -> 1.2}]
mesh["Wireframe"]

instead of the definition of mesh in the question. (It also is necessary to add Needs["NDSolveFEM"], which apparently was omitted from the question by oversight.) It yields a mesh of about 90,000 elements

enter image description here

which is as many as my computer can accommodate without becoming extremely slow. The solution then is

enter image description here

which is less accurate than desired but much better than the result depicted in the question. To obtain still better results, one should use some of the many options in ToElementMesh to generate a mesh with high resolution radially and low resolution azimuthally.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you, that was it! The only question is, how can I speed it up? I need it to work for much larger distances, even 1000 kf^-1. $\endgroup$ – Amit Abir Aug 8 '15 at 19:53
  • $\begingroup$ @AmitAbir To improve speed for very large rmax, it will be necessary to develop a mesh with high resolution in r and low resolution in angle. (I presume that you plan to work with low l and m). One way to do this would be to work in spherical rather than rectangular coordinates. In any case, it is difficult for me, or any reader, to provide suggestions without more detail. To obtain more responses, I suggest that you accept my answer, if it meets your purposes, and then ask a new question. Best wishes. $\endgroup$ – bbgodfrey Aug 8 '15 at 20:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.