6
$\begingroup$

I have the following Schrödinger equation in $2D$:

\begin{cases} \partial_t \Psi(x,t) = V(x,t) \Psi(x,t) \quad x \in [-10,10]^2\\ \Psi(x,0)=\exp( \frac{1}{2} (-x^2+y^2)) \end{cases}

where the potential $V(x,t)=\mathbb{i} \Bigl( \frac{1}{2} \Delta - (x^2+y^2) - \sin^2(t) (x+y) \Bigr)$ with homogeneous Dirichlet boundary conditions. I need the solution at time $T=1$.

Using second order finite differences, I obtain the following plot, plotting $|U|$ at $T=1$:enter image description here

with the following colormap

enter image description here

I'd like to use Mathematica to check my results and to try what comes out by changing some parameters, but I don't know how to solve it properly. Could someone show the plot of the surface I should obtain with Mathematica, and, if possible, the right code-snippet?

EDIT:

I had a different initial data, now my plot seems to agree with the on of Henrik

enter image description here

$\endgroup$
7
  • $\begingroup$ It looks like you try to reproduce some results from our paper onlinelibrary.wiley.com/doi/abs/10.1002/adts.201900011 ? $\endgroup$ Oct 8 '20 at 19:08
  • $\begingroup$ Honestly not, I took it from sci-hub.st/10.1137/S1064827595295337, pg.19 @AlexTrounev $\endgroup$
    – Vefhug
    Oct 8 '20 at 19:38
  • $\begingroup$ I see that you took Example 7.3 extended on 2D. But your example has no stiffness. $\endgroup$ Oct 8 '20 at 22:08
  • $\begingroup$ @AlexTrounev sorry. why does it have no stiffness? $\endgroup$
    – Vefhug
    Oct 9 '20 at 23:34
  • $\begingroup$ If we put like in the paper κ = 10 and µ = 100 then we get at t=1 numerical solution with several pikes and message from the system. In this case method of lines is preferable. $\endgroup$ Oct 10 '20 at 11:50
7
$\begingroup$

Something like the following should do. It employ the finite element method.

Ω =   DiscretizeRegion[Rectangle[{-10, -10}, {10, 10}], MaxCellMeasure -> (1 -> 0.5)];
sol = NDSolveValue[
   {
    D[Ψ[x, y, t], t] == I/2 Laplacian[Ψ[x, y, t], {x, y}] - I ((x^2 + y^2) + (x + y) Sin[t]^2) Ψ[x, y, t], 
    DirichletCondition[Ψ[x, y, t] == 0, True],
    Ψ[x, y, 0] == Exp[-1/2 (x^2 + y^2)]
    },
   Ψ,
   {t, 0, 1},
   {x, y} ∈ Ω
   ];
Plot3D[Abs[sol[x, y, 1]], {x, y} ∈ Ω, PlotRange -> All, AxesLabel -> {"x", "y", "|Ψ|"}]

enter image description here

Looks a bit different from OP's solution, but that could be to a copying error... Anyways, this shows roughly how the PDE can be solved.

For further details (in particular on how to increase the accuracy of the solution), please refer to the documentation (https://reference.wolfram.com/language/FEMDocumentation/tutorial/FiniteElementOverview.html).

enter image description here

Finding the maximum:

NMaximize[{Abs[sol[x, y, 1]], -10 <= x <= 10, -10 <= y <= 10}, {x, y}]

{1.38754, {x -> -0.0632606, y -> -0.0637582}}
$\endgroup$
13
  • $\begingroup$ Many thanks, I'm trying to plot it but it doesn't show anything. Does Mathemtica plots something to you? If so, could you please add it, so that I can accept your answer? @HenrikSchumacher $\endgroup$
    – Vefhug
    Oct 8 '20 at 9:41
  • $\begingroup$ Hm. What is you version of Mathematica? I used NDSolveValue instead of NDSolve because it returns a function an not a rule. But NDSolveValue was only added quite recently... Also there has quite some progress in Mathematica's capabilities, so the version number might be critical. $\endgroup$ Oct 8 '20 at 10:22
  • $\begingroup$ @Vefhug As shown by Henrik, it's quite straightforward to solve the problem in Mathematica, why do you think it's too hard? $\endgroup$
    – xzczd
    Oct 8 '20 at 10:23
  • 1
    $\begingroup$ Better to remove the ∈ Ω in NMaximize, then there's no warning and the output is {1.44277, {x -> 0.127249, y -> 0.127515}}. $\endgroup$
    – xzczd
    Oct 8 '20 at 12:40
  • 1
    $\begingroup$ @HenrikSchumacher You have a typo in equation with I compare to paper, OP and xzczd. $\endgroup$ Oct 10 '20 at 11:32
6
$\begingroup$

FiniteElement isn't necessary for this problem. The old good TensorProductGrid handles the problem quite well:

system = With[{Ψ = Ψ[x, y, t]}, 
          {D[Ψ, t] == I (Laplacian[Ψ, {x, y}]/2 - ((x^2 + y^2) + Sin[t]^2 (x + y)) Ψ),
           Ψ == 0 /. {{x -> -10}, {x -> 10}, {y -> -10}, {y -> 10}},
           Ψ == Exp[-1/2 (x^2 + y^2)] /. t -> 0}];

sol = NDSolveValue[system, Ψ, {t, 0, 1}, {x, -10, 10}, {y, -10, 10}];

Plot3D[Abs@sol[x, y, 1], {x, -10, 10}, {y, -10, 10}, PlotRange -> All, PlotPoints -> 50]

NMaximize[Abs[sol[x, y, 1]], {x, y}]   
(* {1.4014, {x -> -0.0593488, y -> -0.0593488}} *)

Test passes in v12.1.1.


Futher tests show v9.0.1 and v8.0.4 have difficulty in solving the system with defaullt setting, so this turns out to be another example indicating NDSolve is improved silently these years. Nevertheless, with the magic of Pseudospectral, we can still solve the problem in v8 and v9:

If[$VersionNumber < 9, Laplacian = D[#, x, x] + D[#, y, y] &;
  NDSolveValue = #2 /. First@NDSolve[##] &];

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

system = With[{Ψ = Ψ[x, y, t]}, 
          {D[Ψ, t] == I (Laplacian[Ψ, {x, y}]/2 - ((x^2 + y^2) + Sin[t]^2 (x + y)) Ψ),
           Ψ == 0 /. {{x -> -10}, {x -> 10}, {y -> -10}, {y -> 10}},
           Ψ == Exp[-1/2 (x^2 + y^2)] /. t -> 0}];

sol = NDSolveValue[system, Ψ, {t, 0, 1}, {x, -10, 10}, {y, -10, 10}, 
    Method -> mol[55]]; // AbsoluteTiming
(* v8.0.4: {178.4673377, Null} *)
(* v9.0.1: {40.305892, Null} *)

FindMaximum[Abs@sol[x, y, 1], {x, y}]
(* v8.0.4: {1.38975, {x -> -0.0438577, y -> -0.0438577}} *)
(* v9.0.1: lstol warning, {1.38918, {x -> -0.0439239, y -> -0.043924}} *)

NMaximize isn't used to find the maximum because it spits out a Experimental`NumericalFunction[…] as output in v8 and v9, which is obviously a (now fixed) bug.

$\endgroup$
3
  • $\begingroup$ Thanks @xzczd. Tensor product is a FD implementation? $\endgroup$
    – Vefhug
    Oct 8 '20 at 11:46
  • $\begingroup$ @Vefhug FD Yes. (But notice TensorProductGrid is used only for spatial dimension i.e. x and y direction in this problem. It's ODE solver that's used in t direction. ) $\endgroup$
    – xzczd
    Oct 8 '20 at 12:30
  • $\begingroup$ Sure, so it's the classical method of lines :) @xzczd $\endgroup$
    – Vefhug
    Oct 8 '20 at 16:44
3
$\begingroup$

You can simply solve this equation using NDSolve.

Note, I rewrote your equation a bit more towards standard form.

V[x_, y_, t_] := (x^2 + y^2 +  Sin[t]^2 (x + y));
eq = {I  Derivative[0, 0, 1][f][x, y, 
      t] == -Laplacian[f[x, y, t], {x, y}]/2 + V[x, y, t] f[x, y, t], 
   f[x, y, 0] == Exp[-1/2 (x^2 + y^2)], 
   DirichletCondition[f[x, y, t] == 0, True]};
sol = NDSolve[eq, f, {x, -10, 10}, {y, -10, 10}, {t, 0, 1}]

fu[x_, y_] = Abs@f[x, y, 1] /. sol;
Plot3D[fu[x, y], {x, -10, 10}, {y, -10, 10}, PlotRange -> All]

enter image description here

$\endgroup$
4
  • $\begingroup$ @DanielHubner thanks. What is the maximum of $|\Psi|$? Just to see if it's identycal to mine $\endgroup$
    – Vefhug
    Oct 8 '20 at 11:47
  • $\begingroup$ Hi, I get a max of Abs[f[0,0]] = 1.37179 $\endgroup$ Oct 8 '20 at 11:51
  • $\begingroup$ @DanielHubner I obtain the max as $(0,0)$ and it is $1.124$. Is it normal to have so large errors in two dimensional plots? $\endgroup$
    – Vefhug
    Oct 8 '20 at 12:14
  • $\begingroup$ Certainly no. But, I gave you the absolute value of f, not f.Conjugate[f] what is more standard. Could it be that your value comes from Psi^2? $\endgroup$ Oct 8 '20 at 12:38

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.