1
$\begingroup$

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: enter image description here

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?

$\endgroup$
4
  • $\begingroup$ What is your Mathematica version? v12.2 can't reproduce your results $\endgroup$ Sep 25, 2023 at 15:03
  • $\begingroup$ MMA 13.3 reproduces this result! $\endgroup$ Sep 25, 2023 at 17:10
  • $\begingroup$ 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. $\endgroup$ Sep 25, 2023 at 21:36
  • $\begingroup$ Hi, using Mathematica 13.3 on Linux. $\endgroup$
    – kacper
    Sep 26, 2023 at 7:22

1 Answer 1

2
$\begingroup$

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

Needs["NDSolve`FEM`"]
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 *)
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Sep 25, 2023 at 21:44
  • 1
    $\begingroup$ 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? $\endgroup$
    – kacper
    Sep 26, 2023 at 7:42
  • $\begingroup$ 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. $\endgroup$
    – user21
    Sep 26, 2023 at 9:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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