1
$\begingroup$

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["NDSolve`FEM`"];
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]

I obtain: enter image description here

which is exactly what I expected for $\omega$ very small. If I try this with the previous $\omega$ I get this:

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – rhermans Jul 16 '18 at 12:55
  • 1
    $\begingroup$ Definition of eqSyst ??? $\endgroup$ – Ulrich Neumann Jul 16 '18 at 13:11
  • $\begingroup$ This is mistake u[t,x0+S0]== psi0[x0+S0 + S0*Sin[t]]*Exp[-I*Q*t/omega]. $\endgroup$ – Alex Trounev Jul 16 '18 at 13:44
  • $\begingroup$ Sorry, Ulrich, I just added the definition of eqSyst. $\endgroup$ – Stefan Alexandru Ghinescu Jul 16 '18 at 14:47
  • $\begingroup$ Alex Trounev, what do you mean? What's wrong? $\endgroup$ – Stefan Alexandru Ghinescu Jul 16 '18 at 14:49
1
$\begingroup$

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

fig1

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Stefan Alexandru Ghinescu Jul 17 '18 at 6:33
  • $\begingroup$ I have modified the boundary condition for the second problem and now it does not give a warning for inconsistent boundary and initial conditions. $\endgroup$ – Stefan Alexandru Ghinescu Jul 17 '18 at 6:59
  • $\begingroup$ 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. $\endgroup$ – Alex Trounev Jul 17 '18 at 13:11
  • $\begingroup$ 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. $\endgroup$ – Stefan Alexandru Ghinescu Jul 17 '18 at 13:23
  • $\begingroup$ 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. $\endgroup$ – Alex Trounev Jul 17 '18 at 15:52

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.