Mathematica QM Animation of delta potential

[Poirot]

I'm trying to use the Animate function on Mathematica to show a gaussian wave packet passing through a delta potential. I'm quite new to Mathematica and this is by far the hardest thing I've had to do so please bear with me.

I effectively want to solve the integral:
$\phi_k(x) = \left\{ \begin{array}{ll} e^{ikx} + R_k e^{-ikx} & \quad x < 0 \\ T_k e^{ikx} & \quad x > 0 \end{array} \right.$
I have found the reflection and transmission coefficients:
$R_k =\frac{1}{\frac{\hbar^2ik}{sm} -1} \\ T_k =\frac{1}{1-\frac{sm}{\hbar^2ik}}$
where s is the strength of the delta
$V(x) = s\delta(x)$

I then need to plug my phi into the integral
$\int_{-\infty}^{\infty} \phi_k(x) e^{-(k-k_0)^2/\alpha^2} e^{-i\hbar k^2 t/2m} dk$

I've been told the integrals are exactly solvable but not to worry about that just yet so I've been plugging in the limits k-> -100 and 100 in the mean time. I am able to solve the first part of the x<0 by hand but the thing that's causing me the most issues is the transmission and reflection coefficients as the denominator goes to infinity when k=1 (I have set all my constants to 1 in the mean time for simplicity except k0 the initial momentum which I've toggled with to see what happens). Every time I put my piecewise function into Mathematica to integrate it spits out the integral unsolved. I am not really sure how I'm supposed to animate this as I can only solve 1 part of the integral.

So really the first issue I'm having is with the integral, I am *hoping* once I can do that it should all fall into place but I've been stuck for a long time now.

Thanks for any help.

 

[phyzguy]

Have you tried using NIntegrate instead of Integrate? NIntegrate does numerical integration and should have no trouble integrating a piecewise function.

 

[Poirot]

Have you tried using NIntegrate instead of Integrate? NIntegrate does numerical integration and should have no trouble integrating a piecewise function.

I have tried NIntegrate and I get an error "
"The integrand \[Psi][x,t,k] has evaluated to non-numerical values \
for all sampling points in the region with boundaries \
{{-\[Infinity],0.}}"
Which is think is due to the integral containing x and t in the exponents (which aren't to be integrated as I need to plot these after).

 

[phyzguy]

You need to define the NIntegrate function using := instead of = so that it won't try to evaluate it until it needs the value to plot, so then x and t will be defined as numerical values. Also, you can't NIntegrate to infinity, you have to give it a numerical value. Try posting your code (enclosed in Code tags - see below or look at the BBCode guides link below.) so we can try to help.

x=y

 

[Poirot]

You need to define the NIntegrate function using := instead of = so that it won't try to evaluate it until it needs the value to plot, so then x and t will be defined as numerical values. Also, you can't NIntegrate to infinity, you have to give it a numerical value. Try posting your code (enclosed in Code tags - see below or look at the BBCode guides link below.) so we can try to help.

x=y

Ahh I think I see what you mean, I've just given it a go and had a similar error message

En = \[HBar]^2*k^2/2*m;
Rk = 1/(((\[HBar]^2*I*k)/s*m) - 1);
Tk = 1/(1 - (s*m/(\[HBar]^2*I*k)));
R = 1/((2*\[HBar]^2*Energy)/(s^2*m) + 1);
T = 1/(((s^2*m)/(2*\[HBar]^2*Energy)) + 1);
\[HBar] = 1;
m = 1;
s = 1;
\[Alpha] = 1;
k0 = 2;
\[Phi]k[x_, k_] := 
  Piecewise[{{Exp[I*k*x] + Rk*Exp[-I*k*x], x < 0}, {Tk*Exp[I*k*x], 
     x > 0}}];
\[Phi]k1[x_, k_] := Exp[I*k*x] + Rk*Exp[-I*k*x];
\[Phi]k2[x_, k_] := Tk* Exp[I*k*x] ;
f[k_] := Exp[-(k - k0)^2/\[Alpha]^2];
\[Psi][x_, t_, k_] := \[Phi]k[x, k]*Exp[-\[HBar]*k^2*I*t/2*m]*f[k]

g := NIntegrate[\[Psi][x, t, k], {k, -100, 100}]
Animate[Plot[Abs[g]^2, {x, -20, 20}], {t, 0, 20}]

There's quite a lot of preamble as I've been playing around for quite a long time and a fair amount isn't used in this attempt.

I can't immediately see what the problem is but I'm certain there is one!

Thanks for your help by the way, much appreciated!

 

[phyzguy]

I don't have access to Mathematica right now, but don't you need to also define Rk and Tk as functions of k? i.e. Rk[k_] := ... and Tk[k_] := ...

 

[Poirot]

I don't have access to Mathematica right now, but don't you need to also define Rk and Tk as functions of k? i.e. Rk[k_] := ... and Tk[k_] := ...

Ahh you're right thank you.

Just ran and didn't work however :/.

 

[phyzguy]

I also see that you need to define g[x_, t_], not just g. After you've done this, try evaluating some value of g, such as g[5.0, 7.0]. This will test whether you have that much working before you plug it into Animate.

 

[Poirot]

I also see that you need to define g[x_, t_], not just g. After you've done this, try evaluating some value of g, such as g[5.0, 7.0]. This will test whether you have that much working before you plug it into Animate.

Ahh brilliant, the check worked, can take values except x=0 which is to be expected. Although when I implement the Animate function, it comes up saying $Aborted which I think may be due to the integral blowing up?

 

[phyzguy]

I think the $Aborted just means it is taking too long. Try doing some simple plots first, like Plot[g[x,7],{x,-10,10,1}], and try to get a feeling for how long your whole animation is going to take. Then, there are options in Animate to tell it to do the evaluations and save them, only displaying the animations after they are all done.

 

[Poirot]

I think the $Aborted just means it is taking too long. Try doing some simple plots first, like Plot[g[x,7],{x,-10,10,1}], and try to get a feeling for how long your whole animation is going to take. Then, there are options in Animate to tell it to do the evaluations and save them, only displaying the animations after they are all done.

Thank you so much, this is the most progress I've made in days!

It looks like t=7 is towards to end of the animation, and I need to do a bit of tweaking as I think the wave is propagating the wrong way but that should be an easy fix.
So can I get Animate to complete the calculations first and then plot when it's finished?

 

[Poirot]

I think the $Aborted just means it is taking too long. Try doing some simple plots first, like Plot[g[x,7],{x,-10,10,1}], and try to get a feeling for how long your whole animation is going to take. Then, there are options in Animate to tell it to do the evaluations and save them, only displaying the animations after they are all done.

I'm also getting an error message every time I run the plot saying

"NIntegrate failed to converge to prescribed accuracy after 9 \
recursive bisections in k near {k} = {9.76553}. NIntegrate obtained \
1.93516*10^-11+1.15758*10^-11\ I and 1.3159348010134315`*^-15 for the \
integral and error estimates."

and for some values I get a bonus error message

"NIntegrate failed to converge to prescribed accuracy after 9 \
recursive bisections in k near {k} = {10.0143}. NIntegrate obtained \
5.24754*10^-16+4.59702*10^-17\ I and 8.477444514838961`*^-15 for the \
integral and error estimates. "

Not sure if I need to worry as they do pop the graphs out after running for a while eventually.[/Code]

 

[phyzguy]

We're getting to the limit of my knowledge, and I can't currently access my copy of Mathematica. Probably NIntegrate is keeping way more precision than you need for your animation. There are parameters AccuracyGoal and PrecisionGoal that you can set to reduce the required accuracy and speed things up. I'm not sure how to tell Animate to keep running even if it takes a long time, but I'm pretty sure there is a parameter that controls this. Try Googling these things and see what you find. Also, Mathematica has an online documentation site. For example, look at these:

http://reference.wolfram.com/language/ref/AccuracyGoal.html

http://reference.wolfram.com/language/ref/Animate.html

Also, there is a ListAnimate function where you can keep all of the frames and then use ListAnimate to animate them.