School of Computer Science
Computer Science COMP 199 (Winter term)
Excursions in Computing Science
%MATLABpak08cIV: SHOquarter.m, QHO.m
%
% QHO.m QHOampl.m
% function SHOquarter(E,k,m) THM 090309
% Plot x, xdot, Lagrangian, Hamiltonian vs t for phase space ellipse, "thruhalf"
% for the first quarter cycle.
% E total energy, k spring constant, m mass:
% L = m xdot^2/2 - k x^2/2
% H = m xdot^2/2 - k x^2/2 (not necessarily conserved nor total energy)
% x = sqrt(2E/k) f(t)
% xdot = sqrt{2E/m} g(t) constrained xdot = Slope_t x
% Thruhalf trajectory must start at time 0, end at time t2 = (pi/2) sqrt(m/k)
% because ellipse does. It also passes through (Ek/2,Em/2) (see below)
% A good test invocation is SHOquarter(1,2,8/pi^2): this gives E = 1, t2 = 1
function SHOquarter(E,k,m)
Em = sqrt(2*E/m)
Ek = sqrt(2*E/k)
km = sqrt(k/m)
t2 = pi/2/km
res = 50; % resolution of calculation
t = 0:t2/(res-1):t2;
% ellipse
xEllipse = Ek*cos(km*t);
xdotEllipse = -Em*sin(km*t);
Lellipse = m*xdotEllipse.^2/2 - k*xEllipse.^2/2;
Hellipse = m*xdotEllipse.^2/2 + k*xEllipse.^2/2; % This will be conserved
% thruhalf
quintSplineMat = [1,1,1,1;2,3,4,5;1/4,1/8,1/16,1/32;1,3/4,1/2,5/16];
% (from setting f(0) = 0 = f'(0), f(t2) = 1, f'(t2) = sqrt(k/m), f(t2/2) = 1/2,
% f'(t/2) = sqrt(k/m)/2 in f(t) = a(t/t2)^2 + b(t/t2)^3 + c(t/t2)^4 + d(t/t2)^5)
coeff = (quintSplineMat^-1)*[1;pi/2;1/2;pi/4];
a = coeff(1)
b = coeff(2)
c = coeff(3)
d = coeff(4)
xThruhalf = Ek*(1 - a*(t/t2).^2 - b*(t/t2).^3 - c*(t/t2).^4 - d*(t/t2).^5);
xdotThruhalf = -Em*(2/pi)*(2*a*(t/t2) + 3*b*(t/t2).^2 + 4*c*(t/t2).^3 + 5*d*(t/t2).^4);
Lthruhalf = m*xdotThruhalf.^2/2 - k*xThruhalf.^2/2;
Hthruhalf = m*xdotThruhalf.^2/2 + k*xThruhalf.^2/2; % This will not be conserved
% x--xdot ("phase") space
subplot(3,1,1)
plot(xEllipse,xdotEllipse,'b',xThruhalf,xdotThruhalf,'r')
title('Ellipse and "Thruhalf" fit in x--xdot space')
xlabel('x'), ylabel('xdot')
legend('ellipse','thruhalf','location','best')
% plot x, xdot
subplot(3,1,2)
plot(t,xEllipse,'b',t,xdotEllipse,'b-.',t,xThruhalf,'r',t,xdotThruhalf,'r-.')
title('x and xdot for ellipse and "thruhalf" trajectories in x-xdot space')
xlabel('t'), ylabel('x,xdot')
legend('x for ellipse','xdot for ellipse','x for thruhalf','xdot for thruhalf','Location','best')
% plot L, H
subplot(3,1,3)
plot(t,Lellipse,'b',t,Hellipse,'b-.',t,Lthruhalf,'r',t,Hthruhalf,'r-.')
title('L and H for ellipse and "thruhalf" trajectories in x-xdot space')
xlabel('t'), ylabel('L,H')
legend('Lagrangian for ellipse','Hamiltonian for ellipse','Lagrangian for thruhalf','Hamiltonian for thruhalf','Location','best')
% action: use the hand expansion of L on powers of t/t2
% T = 2,4a^2 3,12ab 4,16ac 5,20ad V = 0,1 2,-a 3,-b 4,-c 5,-d
% 4,9b^2 5,24bc 6,30bd 4,a^2 5,ab 6,ac 7,ad
% 6,16c^2 7,40cd 6,b^2 7,bc 8,bd
% *4/pi^2 8,25d^2 8,c^2 9,cd
% 10,d^2
c0 = -1; c2 = 16*(a/pi)^2 + a; c3 = 48*a*b/pi^2 + b;
c4 = (2/pi)^2*(16*a*c + 9*b^2) - a^2 + c;
c5 = (2/pi)^2*(20*a*d + 24*b*c) - a*b + d;
c6 = (2/pi)^2*(30*b*d + 16*c^2) - b^2 - a*c;
c7 = (2/pi)^2*40*c*d - b*c - a*d;
c8 = (2/pi)^2*25*d^2 - c^2 - b*d;
c9 = -c*d;
c10 = -d^2;
action = t2*(c0 + c2/3 + c3/4 + c4/5 + c5/6 + c6/7 + c7/8 + c8/9 + c9/10 + c10/11)
% = -41.4795 * t2
% function QHO THM 090316
% plot amplitudes, probabilities for levels 0--3of quantum harmonic oscillator
% (plots saved as ../QHO.fig. QHO.eps)
% Uses A = QHOampl(n,x,1)
function QHO
n1 = 0;
for n = 3:-1:0
% amplitude part
x = -2*(n+1):0.1*(n+1):2*(n+1);
A = QHOampl(n,x,1);
n1 = n1 + 1;
subplot(4,2,n1)
plot(x,A)
xlabel('x'), ylabel('amplitude')
title(['Q harmonic oscillator amplitude, n = ' num2str(n)])
% probability part using scatter
sizeU = 1000;
Up = rand(1,sizeU); % 1-by-sizeU array
P1 = A.^2;
P = P1/sum(sum(P1)); % probabilities should sum to 1
for k = 1:sizeU
Pindx = cumSumPos(Up(k),P);
qho(k) = x(Pindx);
end %for k
S = 5; % size
C = 'b'; % all blue
n1 = n1 + 1;
subplot(4,2,n1)
scatter(qho,rand(1,sizeU),S,C) % 1D -> 2D, randomly
xlabel('x') %, ylabel('probability')
title([num2str(sizeU) ' Q harmonic oscillators, n = ' num2str(n)])
end % for n
% function P = QHOampl(n,x.alpha) THM 090316
% find quantum harmonic oscillator wavefunction for level n, positions x
% energy (2n+1)/2 \hbar\omega alpha = m*omega/hbar: try alpha = 1
% Used by QHO.m
function P = QHOampl(n,x,alpha)
y = sqrt(alpha)*x;
P0 = sqrt(sqrt(alpha/pi))*exp(-y.^2/2);
switch n
case 0
P = P0;
case 1
P = P0*sqrt(2).*y;
case 2
P = P0.*(2*y.^2 - 1)/sqrt(2);
case 3
P = P0.*(2*y.^3 - 3*y)/sqrt(3);
otherwise disp(['QHOampl: unknown case n=' num2str(n)])
end