# NDSolve Time-dependent Schrodinger Equation

I am currently trying to solve a time dependent Schrodinger equation (TDSE) for a rectangular potential barrier and a varying magnetic potential with a time dependent magnetic potential is oriented on the x-Axis and has a sinusoidal form.

The equation I'm trying to solve is: $$i\omega\partial_{\xi}u(x,\xi) = -\frac{\hbar^2 c^2}{2\mu c^2}\partial_{x,x}u(x, \xi) + V(x)u(x,\xi) + i \frac{\hbar c}{\mu c^2}q* A0*cos(\xi)\partial_{x}u(x,\xi)$$

$\xi = \omega t$.

The initial state of the particle is a stationary state with energy Q > 0 found by:

hbarc = 197.326;
mu = 4000.;
S0 = 3.0;
omega = 0.0001;
ec = 0.0854;
q = 2.*ec;
Q = 5.0;
V0 = 0.;
V1 = 20.;
a = 10.;
l = 5.;
x0 = 0.;
xN = a + l + 10.;
A0 = mu omega S0/(hbarc q);
k = Sqrt[2 mu Q]/hbarc;
potential[x_] = Piecewise[{{V1,a<=x<=a+l},{V0,x<a||x>a+l}};

eqSchUnpert={Q*u0[x]==-hbarc*hbarc/(2*mu)*D[u0[x],{x,2}]+potential[x]*u0[x]};
bc={u0[xN]==beta*Exp[I*k*xN],(D[u0[x],x]/.x->xN)==Exp[I*k*xN]*beta*I*k};
sol=ParametricNDSolve[Join[eqSchUnpert,bc],u0,{x, x0, xN},beta]


and, since it is not a bound state, any beta will do so I choose 10^-4:

psi0=u0[10.^-4]/.sol


Now, when adding the magnetic potential the code takes forever and doesn't bring any result:

eqSchPert = {I*omega*
D[u[t, x],t]==-hbarc*hbarc/(2 mu)D[u[t, x],{x, 2}]+
potential[x] u[t, x] + I*S0*omega*Cos[t]*D[u[t, x], x]};


with the initial and boundary contitions:

icPert = {u[0.,x]== psi0[x]};
bcPert = {u[t,x0+S0]== psi0[x0+S0 + S0*Sin[t]]*Exp[-I*Q*t/omega], u[t,xN-S0]== psi0[xN-S0 + S0*Sin[t]]*Exp[-I*Q*t/omega]};


NDSolve never finishes. If anyone knows what's wrong, please tell me. What I'm doing now is:

eqSyst=Join[eqSchPert, icPert, bcPert];

NDSolve[eqSyst,u,{t,0.,2.},{x,x0+S0,xN-S0},Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","DifferenceOrder"->"Pseudospectral","MinPoints"->300,"MaxPoints"-> 300}}]


EDIT 1

After searching a little more I gave FiniteElement a try and I think I'm getting close, but I'm not there yet. If I set $\omega$ = $10^{4}$ and solve the equation with:

omega = 10.^4;
icPert = {u[r, 0.] == psi0[r]};
bcPert = {u[k*(x0 + S0), t] ==
psi0[k (x0 + S0 + S0*Sin[t])]*Exp[-I*Q*t/omega],
u[k*(xN - S0), t] ==
psi0[k (xN - S0 + S0*Sin[t])]*Exp[-I*Q*t/omega]};
schEq = {I*omega*D[u[r, t], t] == -Q (D[u[r, t], {r, 2}]) + I*k*S0*omega*Cos[t] D[u[r, t], r] + potential[r]*u[r, t]};
eqSyst = Join[schEq, icPert, bcPert];


and I use a mesh made like this:

Needs["NDSolveFEM"];
mesh = ToElementMesh[Rectangle[{k (x0 + S0), 0.}, {k (xN - S0), tFinal}], "MaxBoundaryCellMeasure" -> 0.05, "MeshElementType" -> QuadElement];


the solution becomes:

solution = NDSolve[eqSyst, u, {r, t} \[Element] mesh];
u1 = u /. First[solution];
Plot3D[Abs[Evaluate[u1[r, t]]]^2, {r, k (x0 + S0), k (xN - S0)}, {t, 0., 2*\[Pi]}, PlotRange -> All, Mesh -> All]


which is exactly what I expected for $\omega$ very small. If I try this with the previous $\omega$ I get this: which is absurd. In both cases, I receive a warning message:

NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.


Please note that in the equations I passed to the dimensionless variable $x->\rho = k x$, with $k$ defined in the beggining of this post.

• Hi Stefan and welcome! To make the most of Mma.SE start by taking the tour now. It will help us to help you if you write an excellent question. Edit if improvable, show due diligence by giving a minimal working example of code you have tried yourself. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. – rhermans Jul 16 '18 at 12:55
• Definition of eqSyst ??? – Ulrich Neumann Jul 16 '18 at 13:11
• This is mistake u[t,x0+S0]== psi0[x0+S0 + S0*Sin[t]]*Exp[-I*Q*t/omega]. – Alex Trounev Jul 16 '18 at 13:44
• Sorry, Ulrich, I just added the definition of eqSyst. – Kvothe Jul 16 '18 at 14:47
• Alex Trounev, what do you mean? What's wrong? – Kvothe Jul 16 '18 at 14:49

Minimum working code, in which all equations have solutions. I made a replacement t/omega->t to reduce computation time.

hbarc = 197.326;
mu = 4000.;
S0 = 3.0;
omega = 0.0001;
ec = 0.0854;
q = 2.*ec;
Q = 5.0;
V0 = 0.;
V1 = 20.;
a = 10.;
l = 5.;
x0 = 0.;
xN = a + l + 10.;
A0 = mu omega S0/(hbarc q);
k = Sqrt[2 mu Q]/hbarc;
potential[x_] =
Piecewise[{{V1, a <= x <= a + l}, {V0, x < a || x > a + l}}];
eqSchUnpert = {Q*u0[x] == -hbarc*hbarc/(2*mu)*D[u0[x], {x, 2}] +
potential[x]*u0[x]};
bc = {u0[xN] == beta*Exp[I*k*xN], (D[u0[x], x] /. x -> xN) ==
Exp[I*k*xN]*beta*I*k};
psi0 = ParametricNDSolveValue[{eqSchUnpert, bc},
u0, {x, x0, xN}, {beta}];
Plot[Abs[psi0[10^-4][x]], {x, x0, xN}]
eqSchPert = {I*
D[u[t, x], t] == -hbarc*hbarc/(2 mu) D[u[t, x], {x, 2}] +
potential[x] u[t, x] +
I*S0*omega*Cos[t*omega]*D[u[t, x], x]}; icPert = {u[0., x] ==
psi0[10^-4][x]};
bcPert = {u[t, x0 + S0] == psi0[10^-4][x0 + S0],
u[t, xN - S0] == psi0[10^-4][xN - S0]};

sol1 = NDSolveValue[{eqSchPert, icPert, bcPert},
u[t, x], {t, 0., 2}, {x, x0 + S0, xN - S0},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 300, "MaxPoints" -> 1000,
"DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6];
Plot3D[Abs[sol1], {t, 0, 2}, {x, x0 + S0, xN - S0}, Mesh -> None,
ColorFunction -> Hue, AxesLabel -> {"t", "x"}] • Thank you very much! I know that NDSolve is able to find a solution for short time intervals. I want one the range $\omega$ t from 0 to 2 Pi. This is the reason I post this question. I don't know if this is possible. – Kvothe Jul 17 '18 at 6:33
• I have modified the boundary condition for the second problem and now it does not give a warning for inconsistent boundary and initial conditions. – Kvothe Jul 17 '18 at 6:59
• Thank you, Stefan. I made a calculation with these boundary conditions. To count for long times, we need to reduce the equation to a dimensionless form, determine the region $0\le x \le 2$, and lower the omega to $2\pi / 100$. In this case, there is a chance to count up to $t=100$. I launched such a task, but it has not yet been counted to the end, but it requires fast CPU and a memory of at least 10 GB. – Alex Trounev Jul 17 '18 at 13:11
• Of course we can scale the problem. But it must be scaled all together: the energy of the particle needs to be 4 orders of magnitude higher than $\omega$ and so on. You can choose any system of units as long as the quantities follow the ratio in my original post. – Kvothe Jul 17 '18 at 13:23
• Why do not you do it yourself? Now the task looks like a freak. It is necessary to give it a form convenient for calculation. So I could not calculate the task for t=100. The computer did not have enough memory. Although it seems that the task is simple. But the Mathematica is hanging on it. – Alex Trounev Jul 17 '18 at 15:52