VKeld[kappa_][d_,chi_][x_,y_] := - (1/(4*chi)) * (StruveH[0,Sqrt[x^2 + y^2 + d^2]/(2*Pi*chi/kappa)] - BesselY[0,Sqrt[x^2 + y^2 + d^2]/(2*Pi*chi/kappa)])
Module[
{
nmax = 6, mx = 0.15, my = 1.2, pot = VKeld, kappa = 4.89,
chi = QuantityMagnitude[Quantity[0.41,"Nanometers"],"BohrRadius"],
d = 0 (* Want to also solve for spatially separated electrons and holes, but focusing on the direct case for now *),
eps = 49 * 10^-6 (* This is the smallest value of MaxCellMeasure I can use before my computing cluster crashes due to memory limitations *),
maxi = 10^12,
vn = None (* Since I need to perform the reverse coordinate transform back to cartesian coordinates, no point in having Mathematica try to normalize the transformed eigenfunctions *),
tds,
tdstrans,
shift=10,
stranseps=N[(Pi/2)-$MachineEpsilon],
evs,efs
},
tds = (-(1/2)*((1/mx)*D[f[x, y], {x, 2}, {y, 0}] + (1/my)* D[f[x, y], {x, 0}, {y, 2}]) + pot[kappa][d, chi][x, y]*f[x, y] + shift*f[x, y]);
tdstrans = Simplify[tds /. {f -> (u[ArcTan[#1], ArcTan[#2]] &)} /. {x -> (Tan[\[Xi]]), y -> (Tan[\[Psi]])}, {(Pi/2) > \[Xi] > -(Pi/2), (Pi/2) > \[Psi] > -(Pi/2)}];
{evs,efs}=NDEigensystem[
{
tdstrans,
DirichletCondition[u[\[Xi], \[Psi]] == 0,Abs[\[Xi]] == (stranseps) || Abs[\[Psi]] == (stranseps)]
},
u[\[Xi], \[Psi]],
{\[Xi], \[Psi]} \[Element] Rectangle[{-stranseps, -stranseps}, {stranseps, stranseps}],
nmax,
Method->{"SpatialDiscretization"->{"FiniteElement",{"MeshOptions"->{"MaxCellMeasure"->eps}}},"Eigensystem"->{"Arnoldi","MaxIterations"->maxi},"VectorNormalization"->vn}
];
{
evs-shift,
Head/@efs
}
]
EDIT 4/16:
I realize that my post is quite long and rather open-ended regarding the help I'm looking for, so I'll ask more specific questions here...
- One possible cause of the inconsistencies with
NIntegratecould be theInterpolatingFunctions returned byNDEigensystem. I've read elsewhere thatNDEigensystemcannot perform calculations with arbitrary precision, that it is instead locked to$MachinePrecision.
1a) Are there other options to NDEigensystem that affect the InterpolatingFunctions?
1b) In general, are there other options for NDEigensystem that might yield a better overall result besides changing MaxCellMeasure?
1c) Is there a way to perform post-processing on the InterpolatingFunctions that would somehow make their behavior more predictable, or the results of their integration more reliably accurate?
Perhaps the double-coordinate-transform between Cartesian and tangential coordinates causes trouble with the
InterpolatingFunction. Yesterday I tried working in 2D Polar coordinates so that I'd only have to transform $r \to \tan(\psi)$, $\psi \in [0,\infty); \theta \in [0, 2*\pi]$. However, implementing aPeriodicBoundaryConditionseems to consume a tremendous amount of memory - I have 32 GB RAM and can run the Cartesian-transformed system (whose domain isRectangle[{-Pi/2,-Pi/2},{Pi/2,Pi/2}], for a total area of $\pi^2$) with aMaxCellMeasure = 9 x 10^-6, but for the case of Polar-transformed coordinates, (with domainRectangle[{0,0},{Pi/2,2*Pi}], again an area of $\pi^2$), I run out of RAM at onlyMaxCellMeasure = 5 x 10^-4, a severe drop in accuracy. If I remove thePeriodicBoundaryCondition, I can again runNDEigensystemwithMaxCellMeasure = 9 x 10^-6without memory issues. Is there something fundamental aboutPeriodicBoundaryCondition` that causes it to consume significantly more memory? I can provide the code for how I'm implementing the Polar-transformed system, either in this question or in a separate question, if necessary.For the Cartesian-transformed system, the only boundary condition that I'm imposing is the
DirichletConditionthat the function must go to zero at the boundaries. Since this is a second-order differential equation, I realize that I should have a second boundary condition. Imposing Von Neumann BC, that is, specifying that the first derivative of the function should also go to zero at the boundaries, would seem to make the most sense. I'll note, however, that I did not do this previously when I was only solving the 1D Schodinger equation for isotropic electron/hole masses, yet my results seemed to be fine in that case (e.g. physically accurate, well-converged, etc.). I checked the documentation forNeumannValuebut the implementation of that function seemed to be very different thanDirichletCondition(and the examples given in the docs only described the implementation inDSolve,NDSolve, etc., which takes PDEs in a different form thanNDEigensystem). Can I implement Neumann BC using something likeDirichletCondition[D[f[\[Xi],\[Psi]],\[Xi]]==0,\[Xi]==Pi/2](and a similar statement for the\[Psi]coordinate)? Also, I realize that the derivative will be non-trivial due to the coordinate transformation, but I'm just sketching out the code for now. If that's the correct implementation, can I expect that the eigenfunctions will be more well-behaved?I'll admit that I'm not familiar with the advantages/disadvantages of the specific integration strategies for the
"GlobalAdaptive"method ofNIntegrate. As I mentioned, though, I tried otherMethods such asAdaptiveMonteCarloand did not find that it made my integration results more consistent. Is it possible, though, that different"GlobalAdaptive"strategies would be particularly well-suited to integrating coordinate-transformedInterpolatingFunctions in a large (that is, large enough to be "infinite") Cartesian domain?
