The exact solution for the driven oscillator which we’ll compare to
\operatorname{In}[1]:=\operatorname{Clear}[ m , \omega, \hbar, \Omega, A ] ;
\ln [2]=f\left[t_{-}\right]=A \operatorname{Sin}[\Omega t] ;
\operatorname{In}[3]=\operatorname{xc}\left[ t _{-}\right]=\text {Simplify }[\omega \text { Integrate }[ f [t p] \sin [\omega(t-t p)],\{t p, 0, t\}]] ;
\operatorname{In}[4]=\operatorname{psi}\left[n_{-}, x_{-}\right]:=\left(\frac{m \omega}{\pi \hbar}\right)^{1 / 4} \frac{1}{\sqrt{2^{n} n !}} \text { HermiteH }\left[n, \sqrt{\frac{m \omega}{\hbar}} x\right] \operatorname{Exp}\left[-\frac{1}{2}\left(\sqrt{\frac{m \omega}{\hbar}} x\right)^{2}\right] ;
\ln [5]:=\Psi\left[n_{-}, x_{-}, t_{-}\right]=
\text { Simplify }\left[\text { psi } [ n , x – x c [ t ] ] \operatorname { E x p } \left[\frac { I } { \hbar } \left(-\left(n+\frac{1}{2}\right) \hbar \omega t+m x c^{\prime}[t]\left(x-\frac{x c[t]}{2}\right)+\right.\right.\right.
\left.\left.\frac{m \omega^{2}}{2} \text { Integrate }[f[t p] \times c[t p],\{t p, 0, t\}]\right]\right] ;
\text { Numerics } \Omega=\omega / 5 .
Physical constants
\operatorname{In}[6]= m =1 ; \omega=1 ; \hbar=1 ; \Omega=.2 ; \quad A =1 ;
Numerical constants
\ln [7]=N x=100 ; \text { xmin }=-10 . ; \quad \operatorname{xmax}=10 . ; dx =\frac{x \max -x m i n}{N x+1} ;
Nt =1000 ; d t =\frac{2 \pi / \Omega}{ Nt } ;
Construct H.
The Identity Matrix
\ln [8]=\text { Openone }=\text { SparseArray }\left[\left\{i_{-}, i_{-}\right\} \rightarrow 1 .,\{N x, N x\}\right] ;
The second derivative operator
\operatorname{In}[9]=D 2=\frac{1.0}{d x^{2}} \text { SparseArray }\left[\left\{\left\{i_{-}, i_{-}\right\} \rightarrow-2 .,\left\{i_{-}, j_{-}\right\} / ; A b s[i-j]=1 \rightarrow 1 .\right\},\{N x, N x\}\right] ;
The position operator
\ln [10]=X=\text { SparseArray }\left[\left\{i_{-}, i_{-}\right\} \rightarrow x m i n+d x i,\{N x, N x\}\right] ;
The Hamiltonian
I n[11]:=H\left[t_{-}\right]=-\frac{\hbar^{2}}{2 m} D 2+\frac{1}{2} m \omega^{2} X \cdot X-m \omega^{2} X f[t] ;
Construct the time-evolution operator
The numerator in Cayley’s form
\operatorname{In}[12]=\text { Uplus }\left[t_{-}\right]=\text {Openone }+I H\left[t+\frac{d t}{2}\right] \frac{d t}{2} ;
The denominator in Cayley’s form
\ln [13]==\text { Uminus }\left[t_{-}\right]=\text {Openone }-I H\left[t+\frac{d t}{2}\right] \frac{d t}{2} ;
Solve the t-dep Schrödinger equation
\ln [14]==\operatorname{Psi0}=\text { Eigenvectors }[H[\theta]][[N x]] / \operatorname{Sqrt}[d x] ;
\operatorname{In}[15]:=\operatorname{Psi}=\operatorname{Psi 0} ;
Data =
Table [
(* Step the wave function forward in time *)
\text { Psi }=\text { Linear Solve[Uplus }[(i-1) d t], \text { Uminus }[(i-1) d t] . \text { Psi }] ;
If [
(* Only Plot 100 frames regardless of number of time steps *)
\operatorname{Mod}[i, N t / 100]=0 ,
(* The exact solution *)
g\left[x_{-}\right]=\Phi[0, x, i d t] ;
(* The instantaneous ground state *)
h\left[x_{-}\right]=p s i[0, x-f[i d t]] ;
(* Plotting commands *)
Show [
\{ Plot [\{\operatorname{Re}[g[x]], \quad \operatorname{Im}[g[x]], \quad A b s[g[x]], \quad A b s[h[x]]\} ,
\{x,-3,3\}, \text { PlotRange } \rightarrow\{-1,1\}, \text { AxesLabel } \rightarrow\left\{" x ", "|\Psi|^{\prime \prime}\right\} ,
\text { Plotstyle } \rightarrow \text { \{Blue, Red, Black, \{Black, Dashed\}\}] } ,
\text { ListPlot[ }\{R e[\text { Psi], Im[Psi], Abs [Psi] }\}, \text { PlotRange } \rightarrow\{-1,1\} ,
\text { DataRange } \rightarrow\{ xmin + dx , xmax – dx \} ,
\text { Plotstyle } \rightarrow \text { PointSize[Medium]] } ,
\}
],
Nothing
],
\{\dot{i}, 1, Nt \}
];
\text { In[17]:= ListAnimate [Data] } .
(b)\Omega=5 \omega:
\text { Numerics } \Omega=5 \omega .
Physical constants
\ln [18]:= m =1 ; \omega=1 ; \quad \hbar=1 ; \Omega=5 ; \quad A =1 ;
Numerical constants
\ln [19]=N x=100 ; \text { xmin }=-10 . ; \text { xmax }=10 . ; d x=\frac{x \max -\operatorname{xmin}}{N x+1} ;
Nt =1000 ; d t =\frac{2 \pi / \Omega}{ Nt } ;
Construct H
The Identity Matrix
\operatorname{In}[20]=\text { Openone }=\text { SparseArray }\left[\left\{i_{-}, i_{-}\right\} \rightarrow 1 .,\{N x, N X\}\right] ;
The second derivative operator
\operatorname{In}[21]= D 2=\frac{1.0}{ dx ^{2}} \text { SparseArray }\left[\left\{\left\{ i _{-}, i_{-}\right\} \rightarrow-2 .,\left\{i_{-}, j_{-}\right\} / ; \text {Abs }[i-j]=1 \rightarrow 1 .\right\},\{ Nx , N X \}\right] ;
The position operator
\operatorname{In}[22]==X=\text { SparseArray }\left[\left\{i_{-}, i_{-}\right\} \rightarrow \operatorname{xmin}+d x i,\{N x, N x\}\right] ;
The Hamiltonian
\text { In }[23]=H\left[t_{-}\right]=-\frac{\hbar^{2}}{2 m} D 2+\frac{1}{2} m \omega^{2} X \cdot x-m \omega^{2} X f[t] ;
Construct the time-evolution operator
The numerator in Cayley’s form
\operatorname{In}[24]=\text { Uplus }\left[t_{-}\right]=\text {Openone }+I H\left[t+\frac{d t}{2}\right] \frac{d t}{2} ;
The denominator in Cayley’s form
\operatorname{In}[25]=\text { Uminus }\left[t_{-}\right]=\text {Openone }-I H\left[t+\frac{d t}{2}\right] \frac{d t}{2} ;
Solve the t-dep Schrödinger equation
\operatorname{In}[26]=\operatorname{Psi} 0=\text { Eigenvectors }[H[\theta]][[N x]] / \operatorname{Sqrt}[d x] ;
\operatorname{In}[27]:=\operatorname{Psi}= Psi 0 ;
Data =
Table [
(* Step the wave function forward in time *)
\text { Psi }=\text { Linear Solve[Uplus }[(i-1) d t], \text { Uminus }[(i-1) d t] . \text { Psi] } ;
If [
(* Only Plot 100 frames regardless of number of time steps *)
\operatorname{Mod}[i, N t / 100]=0 ,
(* The exact solution *)
g\left[x_{-}\right]=\Phi[0, x, i d t] ;
(* The instantaneous ground state *)
h\left[x_{-}\right]=\operatorname{psi}[0, x-f[i d t]] ;
(* Plotting commands *)
Show[
\{\operatorname{Plot}[\{\operatorname{Re}[g[x]], \operatorname{Im}[g[x]], \quad \operatorname{Abs}[g[x]], \quad A b s[h[x]]\} ,
\{x,-3,3\}, \text { PlotRange } \rightarrow\{-1,1\}, \text { AxesLabel } \rightarrow\{" x ", "|\Psi| "\} ,
\text { PlotStyle } \rightarrow \text { \{Blue, Red, Black, \{Black, Dashed\} \}] } ,
\text { ListPlot [\{Re[Psi], Im[Psi], Abs[Psi] \}, PlotRange } \rightarrow\{-1,1\} ,
\text { DataRange } \rightarrow\{ xmin + dx , xmax – dx \} ,
\text { Plotstyle } \rightarrow \text { PointSize[Medium]] } .
\}
],
Nothing
],
\{i, 1, \text { Nt }\}
];
\text { In[29]: ListAnimate [Data] } .
\text { Numerics } \Omega=6 / 5 \omega .
Physical constants
\operatorname{In}[30]:= m =1 ; \omega=1 ; \hbar=1 ; \Omega=6 / 5 ; A =1 ;
Numerical constants
\ln [3]:=N x=100 ; \operatorname{xmin}=-10 . j \operatorname{xmax}=10 . j dx =\frac{\operatorname{xmax}- xmin }{ Nx +1} ;
Nt =1000 ; dt =\frac{2 \pi / \Omega}{ Nt } ;
Construct H
The Identity Matrix
\ln [32]=0 \text { penone }=\text { SparseArray }\left[\left\{i_{-}, i_{-}\right\} \rightarrow 1 .,\{N x, N x\}\right] ;
The second derivative operator
\ln [33]: D 2=\frac{1.0}{ d x^{2}} \text { SparseArray }\left[\left\{\left\{i_{-}, i_{-}\right\} \rightarrow-2 .,\left\{i_{-}, j_{-}\right\} / ; \text {Abs }[i-j]=1 \rightarrow 1 .\right\},\{ Nx , Nx \}\right] ;
The position operator
\operatorname{In}[34]:=X=\text { SparseArray }\left[\left\{i_{-}, i_{-}\right\} \rightarrow \operatorname{xmin}+d x i,\{N x, N x\}\right] ;
The Hamiltonian
\operatorname{In}[35]=H\left[t_{-}\right]=-\frac{\hbar^{2}}{2 m} D 2+\frac{1}{2} m \omega^{2} X \cdot X-m \omega^{2} \times f[t] ;
Construct the time-evolution operator
The numerator in Cayley’s form
\ln [36]:=\text { Uplus }\left[t_{-}\right]=\text {OpenOne }+I H\left[t+\frac{d t}{2}\right] \frac{d t}{2} ;
The denominator in Cayley’s form
\operatorname{In}[37]=\text { Uminus }\left[t_{-}\right]=\text {Openone }-I H\left[t+\frac{d t}{2}\right] \frac{d t}{2} ;
Solve the t-dep Schrödinger equation
\ln [38]:=\operatorname{Psi0}=\text { Eigenvectors }[H[\theta]][[N x]] / \operatorname{Sqrt}[d x] ;
\operatorname{In}[39]:=\operatorname{Psi}=\operatorname{Psi} 0 ;
Data =
Table [
(* Step the wave function forward in time *)
\text { Psi }=\text { LinearSolve[Uplus }[(i-1) d t], \text { Uminus }[(i-1) d t] . \text { Psi] } ;
If [
(* Only Plot 100 frames regardless of number of time steps *)
\operatorname{Mod}[i, Nt / 100]=0 ;
(* The exact solution *)
g\left[x_{-}\right]=Ψ[0, x, i d t] ,
(* The instantaneous ground state *)
h\left[x_{-}\right]=p s i[0, x-f[i d t]] ;
(* Plotting commands *)
Show[
\{P \operatorname{lot}[\{\operatorname{Re}[g[x]], \quad \operatorname{Im}[g[x]], \quad A b s[g[x]], \quad A b s[h[x]]\} ,
\{ x ,-3,3\}, \text { PlotRange } \rightarrow\{-1,1\}, \text { AxesLabel } \rightarrow\left\{" x^{\prime \prime}, "|\Psi| "\right\} ,
\text { Plotstyle } \rightarrow \text { \{Blue, Red, Black, \{Black, Dashed\}\}] } ,
\text { ListPlot[\{Re[Psi], Im[Psi], Abs[Psi] }\}, \text { PlotRange } \rightarrow\{-1,1\} ,
\text { DataRange } \rightarrow\{x \min +d x, x \max -d x\} ,
\text { PlotStyle } \rightarrow \text { PointSize[Medium]] } .
\}
],
Nothing
],
\{i, 1, Nt \}
];
ListAnimate[Data]
