# Solving partial differential equation: boundary problem

Here we are solving a set of non-linear equations. And the parameters are

SetDirectory@NotebookDirectory[];
hbar = 1 ;
h = 2 Pi hbar;
num = 8 10^6;
n0 = 1;
q = 47 h;
c2n0 = 39 h;

{\[Omega]x, \[Omega]y, \[Omega]z} = 2 Pi {4.3, 8.8, 420};
rcx = 1;
rcy = 1/Sqrt[M \[Omega]y/hbar];
x0 = 20 rcx;
halfLx = x0;
M = 1;
c2 = c2n0/n0/(hbar \[Omega]x);
c0 = 27.6 c2;
numDiscrete = 200.;
\[Delta]x = 2 halfLx/(numDiscrete);
\[Delta]y = 2 halfLy/(numDiscrete);

datax = Table[x, {x, -halfLx, halfLx, \[Delta]x}];
datay = Table[y, {y, -halfLy, halfLy, \[Delta]y}];



And the initial state for \[Psi]0 is

datapsi0 = ToExpression@Import["datapsi0.csv"];

nx = Length@datapsi0;
ny = Length@First@datapsi0;
data\[Psi]Interpolation =
Catenate[
ParallelTable[{{datax[[ix]], datay[[iy]]},
datapsi0[[ix, iy]] // Chop}, {ix, 1, nx}, {iy, 1, ny}]];
initial\[Psi]1[x_, y_] := Interpolation[data\[Psi]Interpolation][x, y]


The code for solving equations is

Sx = 1/Sqrt[2] ( {
{0, 1, 0},
{1, 0, 1},
{0, 1, 0}
} );
Sy = I/Sqrt[2] ( {
{0, -1, 0},
{1, 0, -1},
{0, 1, 0}
} );
Sz = ( {
{1, 0, 0},
{0, 0, 0},
{0, 0, -1}
} );
vecF = {Sx, Sy, Sz};
vec\[Psi][t_, x_,
y_] := {\[Psi]1[t, x, y], \[Psi]0[t, x, y], \[Psi]minus1[t, x,
y]};
vecFbar[t_, x_,
y_] := {(Conjugate[vec\[Psi][t, x, y]] . Sx .
vec\[Psi][t, x, y]), (Conjugate[vec\[Psi][t, x, y]] . Sy .
vec\[Psi][t, x, y]), (Conjugate[vec\[Psi][t, x, y]] . Sz .
vec\[Psi][t, x, y])};

Vtrap[x_, y_] :=
1/2 M (x^2 + (\[Omega]y/\[Omega]x)^2 y^2); (*normlized potential*)
\[CapitalOmega] = 2 Pi 150;(*Rabi frequency*)
\[Omega]0 = 2 Pi 291 10^3;(*Larmor frequency*)
\[Omega] = \[Omega]0;(*Frequency of RF *)
\[CapitalDelta] = \[Omega] - \[Omega]0;
\[Delta]0 = 1000 2 Pi;
f = 60;
\[Delta][t_] := \[Delta]0 Sin[2 Pi f t/\[Omega]x];
Hs[t_] := (hbar (\[CapitalDelta] - \[Delta][t]) Sz + q Sz . Sz -
hbar \[CapitalOmega] Sx)/(hbar \[Omega]x); (*the factor 1/(hbar \
\[Omega]x) is for dimensionless*)
eqs = {I hbar D[\[Psi]1[t, x, y],
t] == -(hbar^2/(2 M)) Laplacian[\[Psi]1[t, x, y], {x, y}] +
Vtrap[x, y] \[Psi]1[t, x, y] + ((Hs[t] . vec\[Psi][t, x, y])[[
1]]) + c0 (Abs[\[Psi]1[t, x, y]]^2 + Abs[\[Psi]0[t, x, y]]^2 +
Abs[\[Psi]minus1[t, x, y]]^2) \[Psi]1[t, x, y] +
c2 ((vecFbar[t, x, y] . vecF) . vec\[Psi][t, x, y])[[1]],
I hbar D[\[Psi]0[t, x, y],
t] == -(hbar^2/(2 M)) Laplacian[\[Psi]0[t, x, y], {x, y}] +
Vtrap[x, y] \[Psi]0[t, x, y] + ((Hs[t] . vec\[Psi][t, x, y])[[
2]]) + c0 (Abs[\[Psi]1[t, x, y]]^2 + Abs[\[Psi]0[t, x, y]]^2 +
Abs[\[Psi]minus1[t, x, y]]^2) \[Psi]0[t, x, y] +
c2 ((vecFbar[t, x, y] . vecF) . vec\[Psi][t, x, y])[[2]],
I hbar D[\[Psi]minus1[t, x, y],
t] == -(hbar^2/(2 M)) Laplacian[\[Psi]minus1[t, x, y], {x, y}] +
Vtrap[x, y] \[Psi]minus1[t, x,
y] + ((Hs[t] . vec\[Psi][t, x, y])[[3]]) +
c0 (Abs[\[Psi]1[t, x, y]]^2 + Abs[\[Psi]0[t, x, y]]^2 +
Abs[\[Psi]minus1[t, x, y]]^2) \[Psi]minus1[t, x, y] +
c2 ((vecFbar[t, x, y] . vecF) . vec\[Psi][t, x, y])[[3]]};
implicitRegion = Rectangle[{-halfLx, -halfLy}, {halfLx, halfLy}];
bc = {\[Psi]1[0, x, y] == 0, \[Psi]0[0, x, y] ==
initial\[Psi]1[x, y], \[Psi]minus1[0, x, y] == 0,
DirichletCondition[{\[Psi]1[t, x, y] == 0, \[Psi]0[t, x, y] ==
0, \[Psi]minus1[t, x, y] == 0}, True]};

solv = NDSolve[{eqs, bc}, {\[Psi]1, \[Psi]0, \[Psi]minus1}, {t, 0,
2}, {x, y} \[Element] implicitRegion];


The results show

Table[Plot3D[
Evaluate[ Abs[\[Psi]0[t, x, y] /. solv[[1]]]^2], {x, -halfLx,
halfLx}, {y, -halfLy, halfLy}, PlotLabel -> t, PlotRange -> All,
PlotPoints -> 100], {t, 0, 2, .5}]


My question is

The boundary condition DirichletCondition[{\[Psi]1[t, x, y] == 0, \[Psi]0[t, x, y] == 0, \[Psi]minus1[t, x, y] == 0}, True] requires \[Psi]0==0 at the edge, why does the wave function \[Psi]0 grows in the four corners? This phenomenon is unreasonable physcally. And what should we solve this problem？

Note: More mathematical details can be found in this paper.

Decreasing x0 from 20 rcx to 10 rcx and run the code again, we find the fearure of growing in the corners disappers

Table[Plot3D[
Evaluate[ Abs[\[Psi]0[t, x, y] /. solv[[1]]]^2], {x, -halfLx,
halfLx}, {y, -halfLy, halfLy}, PlotLabel -> t, PlotRange -> All,
PlotPoints -> 100], {t, 0, 10, 2}]


And what is the reason here?

@Roland Franzius, I use the PeriodicBoundaryCondition

bc = {\[Psi]1[0, x, y] == 0, \[Psi]0[0, x, y] ==
initial\[Psi]1[x, y], \[Psi]minus1[0, x, y] ==
0, \[Psi]1[t, -halfLx, y] == \[Psi]1[t, halfLx, y], \[Psi]1[t,
x, -halfLy] == \[Psi]1[t, x, halfLy], \[Psi]0[t, -halfLx,
y] == \[Psi]0[t, halfLx, y], \[Psi]0[t, x, -halfLy] == \[Psi]0[t,
x, halfLy], \[Psi]minus1[t, -halfLx, y] == \[Psi]minus1[t,
halfLx, y], \[Psi]minus1[t, x, -halfLy] == \[Psi]minus1[t, x,
halfLy](*,DirichletCondition[{\[Psi]1[t,x,y]==0,\[Psi]0[t,x,y]==
0,\[Psi]minus1[t,x,y]==0},True]*)};



The results show

solv = NDSolve[{eqs, bc}, {\[Psi]1, \[Psi]0, \[Psi]minus1}, {t, 0,
8}, {x, -halfLx, halfLx}, {y, -halfLy, halfLy}(*,{x,
y}\[Element]implicitRegion*)];
Table[DensityPlot[
Evaluate[ Abs[\[Psi]0[t, x, y] /. solv[[1]]]^2], {x, -halfLx,
halfLx}, {y, -halfLy, halfLy}, PlotLabel -> t, PlotRange -> All,
PlotPoints -> 100], {t, 0, 8, 2}]


Although the wave function \[Psi]0 doesn‘t grow in the corners, but the obvious expanding feature is not desired.

So where is the problem?

• Perhaps, this should be edited into your question or asked as a separate question. Commented Mar 26, 2023 at 18:42

You model a system of three fields, each obeying Schödinger equation with a osillator potential and interacting via the third order interaction term (as far as I see)

i d_t psi_k == - Lap psi_k + V psi_k + WW(psi_1,psi_2,psi_3)


The one-particle linear potential term Vtrap*psi is a confining for all energies, so the psi_k move in the pit insofar as the total energy conserved.

But even with energy conservation for all fields, I dont't see on inspection, if the scattering is producing higher frequency components, that naturally blur out a spatial concentrated initial field.

This fact is behind the spatial wave-packet spreading of a Gaussin wave packet, that is a Fourier composition of wave components of all energies in (0, infinity)

Comment on code: If a function like Hs depends on variable t,x you should define it as Hs[t_,x_] in order to enable the NDSolve algorithm to read and transfers a transparent map to its code. Otherwise the first argument of NDSolve should be enclosed in Evaluate[] in order to make dependenecies transparent.

• Your answer probably would be more appropriate as a comment. Commented Mar 27, 2023 at 1:39