Problem #4 - Wave Packet in a One-Dimensional Box

Define the initial condition for Psi

Clear[Sigma,k0,Psi];
Sigma[a_]:= a/8;
k0[a_]:= 10/a;
Psi0[x_,a_]:=(1/(2Pi*Sigma[a]^2))^(1/4) * 
	Exp[-((x-a/2)^2)/(2*Sigma[a]^2) + I*k0[a]*x];
Psi0[x,a]

Checking the normalization condition

NIntegrate[Abs[Psi0[x,1]]^2,{x,-Infinity,Infinity}]

We note that the initial condition is not normalized. We can normalize
it by multiplying the expression for Psi[] by 2^(1/4)

Clear[Sigma,k0,Psi0];
Sigma[a_]:= a/8;
k0[a_]:= 10/a;
Psi0[x_,a_]:=2^(1/4)*(1/(2Pi*Sigma[a]^2))^(1/4) * 
	Exp[-((x-a/2)^2)/(2*Sigma[a]^2) + I*k0[a]*x];
Psi0[x,a]

NIntegrate[Abs[Psi0[x,1]]^2,{x,-Infinity,Infinity}]

Here is a sample plot of Psi to check our work

Plot[Re[Psi0[x,1]],{x,0,1}]

Define the stationary solution Phi

Clear[Phi];
Phi[x_,a_,n_]:= Sqrt[2/a]*Sin[n*Pi*x/a];
Phi[x,a,n]

Define the initial co-efficients Cn0

Clear[Cn0];
Cn0[n_,x_,a_]:=
	NIntegrate[Psi0[x,a]*Phi[x,a,n],
		{x,-Infinity,Infinity}];

We note that we can choose any value for a, since
the result for Cn0[] does not change with a; a is
merely a scaling factor for the length of the box.

Print[Cn0[1,x,1]];   (* a = 1 *)
Print[Cn0[1,x,2]];   (* a = 2 *)

Calculation of the first 10 co-efficients

Clear[CnList];
CnList = Table[Cn0[n,x,1],{n,1,10}];
Print["The following are the Cn(0)'s: "];
Do[Print["n = ", n, ": ", CnList[[n]]] ,{n,1,10}]

Simulation of the wave in time

Clear[En,Psi];
h = 6.63*10^-34;
m = 9.11*10^-31;
a = 1;
En[n_,a_]:=(h^2*Pi^2*n^2)/(2*m*a^2);
Psi[x_,t_]:=Sum[Exp[-I*En[n,a]*t/h]*Phi[x,a,n]*
	CnList[[n]],{n,1,10,1}];
tStep = .02*m*a^2/h;
Do[Plot[Re[Psi[x,t]],{x,0,a},PlotRange -> {-2,2}],
	{t,0,10*tStep,tStep}]

tStep

In this setup, x0=a/2, or the middle of the box. Decreasing x0 shifts the graph to the left.
k0 is the wave number of the exponential in x. Decreasing k means reducing the number
of waves per unit x, and therefore having a longer wavelength.

Clear[Sigma,k0,Psi0];
Sigma[a_]:= a/8;
k0[a_]:= 10/a;
Psi0[x_,a_]:=2^(1/4)*(1/(2Pi*Sigma[a]^2))^(1/4) * 
	Exp[-((x-a/4)^2)/(2*Sigma[a]^2) + I*k0[a]*x];
Plot[Re[Psi0[x,1]],{x,0,1}] (* smaller x0 *)

(* The result is centered at x = 0.25 *)

Clear[Sigma,k0,Psi0];
Sigma[a_]:= a/8;
k0[a_]:= 1/a;
Psi0[x_,a_]:=2^(1/4)*(1/(2Pi*Sigma[a]^2))^(1/4) * 
	Exp[-((x-a/2)^2)/(2*Sigma[a]^2) + I*k0[a]*x];
Plot[Re[Psi0[x,1]],{x,0,1}] (* smaller k0 *)

(* The result only has 1 wave *)