# Time dependent Schrödinger equation in 2D

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$$:

with the following colormap

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

• It looks like you try to reproduce some results from our paper onlinelibrary.wiley.com/doi/abs/10.1002/adts.201900011 ? – Alex Trounev Oct 8 '20 at 19:08
• Honestly not, I took it from sci-hub.st/10.1137/S1064827595295337, pg.19 @AlexTrounev – Vefhug Oct 8 '20 at 19:38
• I see that you took Example 7.3 extended on 2D. But your example has no stiffness. – Alex Trounev Oct 8 '20 at 22:08
• @AlexTrounev sorry. why does it have no stiffness? – Vefhug Oct 9 '20 at 23:34
• 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. – Alex Trounev Oct 10 '20 at 11:50

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", "|Ψ|"}]


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).

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}}

• 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 – Vefhug Oct 8 '20 at 9:41
• 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. – Henrik Schumacher Oct 8 '20 at 10:22
• @Vefhug As shown by Henrik, it's quite straightforward to solve the problem in Mathematica, why do you think it's too hard? – xzczd Oct 8 '20 at 10:23
• Better to remove the ∈ Ω in NMaximize, then there's no warning and the output is {1.44277, {x -> 0.127249, y -> 0.127515}}. – xzczd Oct 8 '20 at 12:40
• @HenrikSchumacher You have a typo in equation with I compare to paper, OP and xzczd. – Alex Trounev Oct 10 '20 at 11:32

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 ExperimentalNumericalFunction[…] as output in v8 and v9, which is obviously a (now fixed) bug. • Thanks @xzczd. Tensor product is a FD implementation? – Vefhug Oct 8 '20 at 11:46 • @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. ) – xzczd Oct 8 '20 at 12:30 • Sure, so it's the classical method of lines :) @xzczd – Vefhug Oct 8 '20 at 16:44 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] ` • @DanielHubner thanks. What is the maximum of$|\Psi|$? Just to see if it's identycal to mine – Vefhug Oct 8 '20 at 11:47 • Hi, I get a max of Abs[f[0,0]] = 1.37179 – Daniel Huber Oct 8 '20 at 11:51 • @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? – Vefhug Oct 8 '20 at 12:14
• 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? – Daniel Huber Oct 8 '20 at 12:38