# NDSolveValue result contradicts initial condition

I'm trying to solve the Schroedinger equation in a box with two holes defined by:

l = 10.0;
th = 0.1;

rg =
RegionDifference[
RegionDifference[
RegionDifference[
Rectangle[{-l , -l} , {l , l}] ,
Rectangle[{-th , -0.5} , {th , 0.5}]
],
Rectangle[{-th , 0.6} , {th , l + 0.9}]
],
Rectangle[{-th , -0.6} , {th , -l-0.9}]
];


the region looks like this:

The numerical solution is obtained using:

schEq = I D[psi[x , y , t] , {t , 1}] +
1/2 Laplacian[psi[x , y , t] , {x , y}] == 0;
ic = psi[x, y , 0] == 10.0 Exp[-2 ((x + 5.0)^2 + y^2)];

sol = NDSolveValue[
{
schEq,
ic,
DirichletCondition[psi[x , y , t] == 0 , True]
},
psi,
{x , y} \[Element] rg , {t , 0 , 0.1}];


This executes without problems, however:

sol[-5.0 , 0.0 , 0.0]


evaluates to 8.80568 instead of 10. This seems to contradict the initial conditions. Why is this value different from 10?

• What is your Mathematica version? v12.2 can't reproduce your results Sep 25, 2023 at 15:03
• MMA 13.3 reproduces this result! Sep 25, 2023 at 17:10
• Note: on MM 13.1 (MacOS, M1 chip) the region creation fails with the code given. However, rg = RegionDifference[Rectangle[{-l, -l}, {l, l}], RegionUnion[Rectangle[{-th, -0.5}, {th, 0.5}], Rectangle[{-th, 0.6}, {th, l + 0.9}], Rectangle[{-th, -l - 0.9}, {th, -0.6}]]]; gives the same region and works fine. Sep 25, 2023 at 21:36
• Hi, using Mathematica 13.3 on Linux. Sep 26, 2023 at 7:22

Let's have a look at the mesh that is automatically generated:

Needs["NDSolveFEM"]
mesh = ToElementMesh[rg]


Now, we evaluate the ic at the coordinates of that default mesh:

data = Function[{x, y}, 10. E^(-2 ((5. + x)^2 + y^2))] @@@
mesh["Coordinates"];


Looking at the MinMax we get:

MinMax[data]
(* {5.11195*10^-282, 9.22845} *)


So, we see that the larges value is not even 10. For the interpolation of the values you will always need several values close by - those are even smaller such that we do not get to the value of 10. That's just how interpolation works.

One thing you could do is add a mesh coordinate that {-5,0}.

mesh = ToElementMesh[rg, "IncludePoints" -> {{-5, 0}}];
data = Function[{x, y}, 10. E^(-2 ((5. + x)^2 + y^2))] @@@
mesh["Coordinates"];
MinMax[data]
(* {5.11195*10^-282, 10.} *)


But what you really want is a finer then the default mesh:

mesh = ToElementMesh[rg, "MaxCellMeasure" -> 0.05];
sol = NDSolveValue[{schEq, ic,
DirichletCondition[psi[x, y, t] == 0, True]},
psi, {x, y} \[Element] mesh, {t, 0, 0.1}];
sol[-5.0, 0.0, 0.0]
(* 9.99076 - 4.42003*10^-18 I *)

• It's also instructive to superimpose the contour for which the initial wave function is equal to 9 on top of the mesh: Show[MeshRegion[mesh], ContourPlot[ic[[2]] == 9, {x, -10, 10}, {y, -10, 10}, ContourStyle -> Red]] On my machine there's maybe one mesh point inside this contour. This also shows why the finer mesh helps resolve the issue. Sep 25, 2023 at 21:44
• Thank you everyone! This was very instructive. I'm thinking of diving a little deeper into PDEs in Mathematica. Is the electronic documentation a good resource? Sep 26, 2023 at 7:42
• Good question: but my opinion might be slightly biased - as I wrote substantial parts of the FEM/PDEModeling stuff ;-) but seriously, if find things that are not clear, let me know and I'll see if I can explain them better. Sep 26, 2023 at 9:26