Modern Physics (comments
to: moloney@nextwork.rose-hulman.edu)
MAPLE text and input lines
PHYSICS RESOURCE PACKETS PROJECT File: BARRIER.MS
Rose-Hulman Institute of Technology Authors: Dr. Kirtley & Perry Peters
http://www.rose-hulman.edu/~moloney Software: Maple V Release 3
The idea is to create a wave packet of several waves and watch the packet
develop with time. The individual waves are not just plane waves, but rather
each one represents a plane wave interacting with a barrier. The full wave
(coming in, reflected, inside barrier, going out the other side) is all
represented in one function using heaviside (step) functions.
Solution:
restart;
assume(A>0,B>0,F>0,k>0,m>0,hbar>0);
assume(x,real);
assume(t,real);
Constant outside the barrier
k:=sqrt(2*m*Eng)/hbar;
Constant inside the barrier
j:=sqrt(2*m*(V-Eng))/hbar;
Let us define the waves in the three regions: one to the left of the well, one
in the center of the barrier, and the third to the right of the barrier. Let us
call them psi1, psi2, and psi3 respectively.
psi1:=A*exp(I*k*x)+B*exp(-I*k*x);
psi2:=C*exp(j*x)+D*exp(-j*x);
psi3:=F*exp(I*k*x);
Let us now define the entire wave function in all three regions in terms of a
step function.
psi:=psi1*Heaviside(a-x)+psi2*(Heaviside(x-a)-Heaviside(x-b))+psi3
*Heaviside(x-b);
Let us now define the following boundary conditions: if we assume that the well
extends from x=a to x=b, then the proper wave functions as well their
derivatives are continuous at these two boundaries.
eqn1:=subs(x=a,psi1=psi2);
eqn2:=subs(x=a,diff(psi1,x)=diff(psi2,x));
eqn3:=subs(x=b,psi2=psi3);
eqn4:=subs(x=b,diff(psi2,x))=subs(x=b,diff(psi3,x));
eqnsset:={eqn1,eqn2,eqn3,eqn4};
Now to find the coefficients B, C, D, and F in terms of A
ceff1:=solve(eqnsset,{B,C,D,F});
Values to later be substituted in for the variables.
values:=(a=0,b=1*10^(-9),A=1,V=1.4*10^(-19),hbar=1.055*10^(-34),m=9.1*10^(-31),
v=6*10^5);
Substitute "values" into B, C, D, and F
ceff:=simplify(subs(values,ceff1));
Substitute B, C, D, and F into the psi.
psifin:=subs(ceff,psi);
Gaussian wave packet for an electron traveling with a speed of 5 x 10^5 m/s.
The packet is made up of N fourier components around the central k value (ko).
This is in effect a discrete fourier transform, so that one will have aliasing
of the wave function at twice the 'sampling frequency'. The group velocity is
twice the phase velocity for nonrelativistic particles, so one should be able to
see the envelope overtaking individual peaks as the packet moves along. Any
numerical integration will have some graininess to it, and so some aliasing is
to be expected in any case.
ko := m*v/hbar;
delta_k := ko/10;
c1:=hbar/(2*m);
c2:=1/(2*delta_k)^2;
kRange:=4*delta_k;
Number of fourier components
N:=4;
Find the energies for each of the components
for i from 1 to N do
k.i:=subs(values,ko-kRange/2 + (i-1)*(kRange/N));
E.i:=subs(values,1/2*k.i^2*hbar^2/m);
od;
i:='i';
Equation for the wave packet in terms of t, E, and k.
packet:=subs(values,exp(-I*c1*(kg)^2*t)*psifin);
Substitute in the component energies and k values and add up the fourier
components.
s:=convert({seq(subs({kg=k.i,Eng=E.i},packet),i=1..N)},`+`);
Probability of Psi
psistarpsi:=evalf(s*conjugate(s));
The real part if Psi
realpsival:=Re(evalf(s));
The imaginary part of Psi
impsival:=Im(evalf(s));
Relative energy levels of barrier and fourier components.
scalefactor:=5*10^19;
levels:=seq(scalefactor*E.i,i=1..N);
barrier:=subs(values,scalefactor*V*(Heaviside(x-a)-Heaviside(x-b)));
The rectangular pulse is the relative energy level of the potential barrier.
The four horizontal lines are the relative energies of the fourier components.
plot({subs(t=0,psistarpsi),barrier,levels},x=-1.5*10^(-8)..1.5*10^(-8),
numpoints=300,title=`|psi|^2`);
plot({subs(t=0,realpsival),levels,barrier},x=-1.5*10^(-8)..1.5*10^(-8),
numpoints=300,title=`Re[psi]`);
plot({subs(t=0,impsival),barrier,levels},x=-1.5*10^(-8)..1.5*10^(-8),
numpoints=300,title=`Im[psi]`);
Load the plots package.
with(plots):
Animation showing the motion of the wave packet.
animate({psistarpsi,barrier,levels},x=-1.5*10^(-8)..1.5*10^(-8),t=0..2*10^(-14),
numpoints=300,frames=20,color=red,title=`|Psi|^2`);
Zoom in on the previous plot to view the motion of the wave packet at the
barrier.
animate({psistarpsi,barrier,levels},x=-5*10^(-9)..5*10^(-9),t=0..2*10^(-14),
numpoints=300,frames=20,color=red,title=`|Psi|^2`);
PHYSICS RESOURCE PACKETS PROJECT File: LEGENDRE.MS
Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney Software: Maple V Release 3
Objective: To derive the first eight Legendre polynomials in two different ways:
a) Using Maple's built in Legendre function
b) By using the formula given: 1/(2^n*(n!))*(d/dx)^n[(x^2 - 1)^n]
Solution:
Load the orthopoly package
with(orthopoly);
Part A
Generate the first eight Legendre polynomials
p1 := P(1,x);
p2 := P(2,x);
p3 := P(3,x);
p4 := P(4,x);
p5 := P(5,x);
p6 := P(6,x);
p7 := P(7,x);
p8 := P(8,x);
Part B
By using a simple FOR loop, the formula given at the top can be evaluated for
several values of n.
for n from 1 to 8 do
P.n:= simplify(1/(2^n*(n!))*diff((x^2 - 1)^n, x$n));
od;
Plot the first eight Legendre polynomials
plot({P.(1..8)},x=-1.5..1.5,y=-1.5..1.5,title=`Plot of first eight Legendre
Polynomials`);
PHYSICS RESOURCE PACKETS PROJECT File: LGHTBULB.MS
Rose-Hulman Institute of Technology Authors: Perry Peters & Greg Williby
http://www.rose-hulman.edu/~moloney Software: Maple V Release 3
A light bulb is modeled as a blackbody object.
a) Find the energy distribution for the light bulb using Planck's radiation
formula for T = 2000K to T = 7000K.
b) Find the frequency and wavelength of maximum energy output of the light
bulb if T = 2800K
c) What percentage of the light bulb's total energy is emitted in the visible
spectrum (350nm - 750nm)?
Solution:
Part A
Clear all variables.
restart;
Load the plots package
with(plots):
Planck's radiation formula in eV
u:=6.2415E18*8*Pi*h*nu^3/(c^3*exp(h*nu/(k*T))-1);
define constants
h:=6.626*10^(-34);
k:=1.381*10^(-23);
c:=2.998*10^8;
Display the function is decimal form.
evalf(u);
reset variable
T:='T':
Generate functions of u for different values of T
for q from 2 to 7 do
u.q:=subs(T=q*1000,u);
od;
Plot these function on one graph.
plot({u.(2..8)},nu=0..20*10^14,title=`Energy Distribution Plots for
T=2000K-7000K`);
plot3d(u,nu=0..20*10^14,T=0..8000,title=`3D Graph of Energy Distribution`);
Show animation of energy distribution when temperature changes.
T:='T':
animate(u,nu=0..20*10^14,T=0..8000,color=red);
Part B
Find the frequency of maximum energy at a given temperature.
T:=2800;
up:=diff(u,nu):
evalf(up);
Maximum output is where derivative of u is 0.
Graph u' with a scale factor to make it more visible.
plot(up*10^15,nu=0..10^15,title=`Graph of u'(nu) at T=2800`);
Find the frequency where the derivative is 0.
frequency:=fsolve(up=0,nu=5E13..2E14);
Find the corresponding wavelength.
wavelength:=c/frequency;
Convert to nm
wavelength_nm:=10^9*wavelength;
Part C
Find the frequencies corresponding to 750 nm and 350 nm.
f1:=c/750E-9;
f2:=c/350E-9;
Energy in the visible spectrum
Energy_visible:=evalf(int(u,nu=f1..f2));
Total energy
Percentage of energy in the visible spectrum
Energy_visible/Energy_tot*100;
PHYSICS RESOURCE PACKETS PROJECT File: ORBITAL.MS
Rose-Hulman Institute of Technology Author: Perry Peters
http://www.rose-hulman.edu/~moloney Software: Maple V Release 3
Plots 's', 'p', & 'd' orbital shells using the implicitplot3d command.
restart;
with(plots):
assume(r,real);
assume(theta,real);
assume(phi,real);
assume(a0,real);
used in radius range
scale:=10^(-95/10);
100
Enter in the wave function.
psi:=exp(-r/a0)/(sqrt(Pi)*a0^(3/2));
psi2:=psi*conjugate(psi);
Switch theta and phi for the coordinate system that Maple uses.
psi2sw:=subs({theta=phi,phi=theta},psi2);
Substitute in the value for a0
psi2sub:=subs(a0=5282/1000*10^(-11),psi2sw);
This plots a surface of uniform probability. Change "10^27" below to a
different value to see a surface with a different probability.
implicitplot3d(psi2sub=10^(27),theta=0..2*Pi,phi=0..Pi,r=0..scale*.75,
coords=spherical,grid=[15,15,15],title=`n=1, l=0, m=0`);
210
psi:=r*exp(-r/(2*a0))*cos(theta)/(4*sqrt(2*Pi)*a0^(3/2)*a0);
psi2:=psi*conjugate(psi);
Switch theta and phi for the coordinate system that Maple uses.
psi2sw:=subs({theta=phi,phi=theta},psi2);
Substitute in the value for a0
psi2sub:=subs(a0=5282/1000*10^(-11),psi2sw);
implicitplot3d(psi2sub=10^(26),theta=0..2*Pi,phi=0..Pi,r=0..3*scale,
coords=spherical,grid=[20,20,20],title=`n=2, l=1, m=0`);
320
psi:=r^2*exp(-r/(3*a0))*(3*cos(theta)^2-1)/(81*sqrt(6*Pi)*a0^(3/2)*a0^2);
psi2:=psi*conjugate(psi);
Switch theta and phi for the coordinate system that Maple uses.
psi2sw:=subs({theta=phi,phi=theta},psi2);
Substitute in the value for a0
psi2sub:=subs(a0=5282/1000*10^(-11),psi2sw);
implicitplot3d(psi2sub=10^(27),theta=0..2*Pi,phi=0..Pi,r=0..2.3*scale,
coords=spherical,grid=[15,15,15],title=`n=3, l=2, m=0`);
Find the probability that an electron exists around an atom.
3D spherical integration
Int(r^2*R^2,r=0..infinity)*Int(Theta^2*sin(theta),theta=0..Pi)
*Int(Phi^2,phi=0..2*Pi);
a0:=5282/1000*10^(-11);
n=1 l=0 m=0
Phi:=1/sqrt(2*Pi);
Theta:=1/sqrt(2);
R:=2*exp(-r/a0)/a0^(3/2);
ans:=int(r^2*R^2,r=0..infinity)*int(Theta^2*sin(theta),theta=0..Pi)
*int(Phi^2,phi=0..2*Pi);
PHYSICS RESOURCE PACKETS PROJECT File: PLANCK.MS
Rose-Hulman Institute of Technology Author: Dr. Paul Mason
http://www.rose-hulman.edu/~moloney Software: Maple V Release 3
The Planck function has two wavelengths at which the emitted power is exactly
half of the peak power. Determine a temperature dependence for each of these
wavelengths (similar to Wien's displacement law) and find the fraction of the
total power which is emitted between these two wavelengths.
Solution:
restart;
Planck's Law is given by:
Wlambda := (2*Pi*h*c^2)/(lambda^5)/(exp((c*h)/(k*(lambda*T))) - 1);
Where the following constants are defined:
c := 2.99792458e8;
h := 6.6260755e-34;
k := 1.380658e-23;
The first step is to find the peak wavelength. This is simply a matter of
finding the wavelength at which the derivative of the Planck function is zero.
diff1:= diff(Wlambda,lambda)=0;
Maxlambda:= solve(diff1,lambda);
Substituting the peak wavelength back into the Planck function allows us to find
the peak power emitted. It turns out that the peak power obeys its own version
of Wien's displacement law:
PeakPower := evalf(subs(lambda=Maxlambda,Wlambda));
Now, to find the wavelengths at which the power is only half of this peak value
we need to solve:
eqn2 := evalf(Wlambda-PeakPower/2)=0;
Unfortunately, the solve command can't handle this formula because of the fifth
power. To get around this limitation, we use fsolve and specify the
temperature and wavelength range to look for a solution. First, specify a
temperature.
T := 300;
Next, look at the peak wavelength to get a clue as to the appropriate wavelength
ranges to look.
Maxlambda;
Now, try the fsolve command.
l1 := fsolve(eqn2,lambda,.8e-7...8e-6);
This is the shorter of the two wavelengths. Try again for the longer
wavelength.
l2 := fsolve(eqn2,lambda,.8e-6...8e-5);
The next step is to find the power emitted between the two wavelengths.
mid := evalf(int(Wlambda,lambda=l1..l2));
And from the long wavelength to infinity.
long := evalf(int(Wlambda,lambda=l2..infinity));
Finally, we find the total power emitted (but we need a substitution before this
integration will work).
temp := subs(lambda = x / T, Wlambda);
SB := evalf(int(temp/T, x=0..infinity));
midfrac := mid / SB;
So about 63% of the total energy is emitted between the two half-power
wavelengths.
longfrac := long / SB;
33% is emitted at long wavelengths. That only leaves
1-midfrac-longfrac;
T := 'T';
about 4% at short wavelengths. These fractions are independent of temperature.
This means that the two half-power wavelengths should obey Wien displacement
law type formulae.
shortl := l1/300*T;
longl := l2/300*T;
PHYSICS RESOURCE PACKETS PROJECT File: PSHELL.MS
Rose-Hulman Institute of Technology Author: Perry Peters
http://www.rose-hulman.edu/~moloney Software: Maple V Release 3
Generate a density plot of a p orbital shell.
restart;
with(plots):
a0:=5.282E-11;
Wave function in spherical coordinates.
psi:=evalf(1/(4*sqrt(2*Pi)*a0^(3/2))*r*exp(-r/(2*a0))*cos(theta));
Wave function in rectangular coordinates.
psixy:=subs({r=sqrt(x^2+y^2),cos(theta)=y/sqrt(x^2+y^2)},psi);
Plot range.
Range:=9.3;
densityplot(psixy^2,x=-10^(-Range)..10^(-Range),y=-10^(-Range)..10^(-Range),
grid=[80,80]);
PHYSICS RESOURCE PACKETS PROJECT File: RAYPLANCK.MS
Rose-Hulman Institute of Technology Author: Dr. Mike Moloney
http://www.rose-hulman.edu/~moloney Software: Maple V Release 3
Rayleigh-Jeans and planck radiation formulas.
Plot the spectral energy density of a perfect blackbody over a large range of
frequencies, at a particular temperature, T. To do this, assume that the inside
walls of the body are perfect reflectors for a series of standing
electromagnetic waves. The number of independent standing waves in the
frequency range f and df turns out to be (8 pi f^2df)/c^3. Multiply this factor
by kT, the average energy of each standing wave, where k is Boltzmann's
constant, to obtain the spectral energy density u(f)df. This is Rayleigh-Jeans
function for blackbody radiation.
At high frequencies Rayleigh-Jeans theory breaks down (ultraviolet catastrophe).
This led to the discovery of the quantum picture of radiation by Planck. In the
new theory, the spectral density is given by
u(f)df = (8pi h/c^3)(f^3 df /(exp(hf/kT)-1)).
a) Go to a new, dimensionless variable x = hf/kT and obtain the Rayleigh-Jeans
formula in terms of x.
b) Put the Planck formula in terms of x.
c) Do a series expansion of the Planck formula in x. Show that when only the
first term in x is retained in the expression exp(x)-1 the expressions in a)
and b) are the same. [This requires expanding the planck formula to order
x^3.]
d) Plot Planck and Rayleigh-Jeans spectral densities and find the region in x
where the two agree. Is this region a high-frequency or a low-frequency
region (hf>>kT or hf<< kT) ?
Solution:
restart;with(plots):
Rayleigh:=kT*8*Pi*f^2*df/c^3;
Planck:=8*Pi*h*f^3*df/(c^3*(exp(h*f/kT)-1));
Raysub:=subs(f=x*kT/h,df=dx*kT/h,Rayleigh);
Plnksub:=subs(f=x*kT/h,df=dx*kT/h,Planck);
ex1:=series(Plnksub,x,4);
n:=op(1,ex1); # see what the first operand is
n*x^2/Raysub;
Comparing the integrals comes down to x^2 dx vs x^3dx/(exp(x)-1).
rayplot:=plot(x^2,x=0..2,color=red):
plnkplot:=plot(x^3/(exp(x)-1),x=0..2,color=blue):
display({rayplot,plnkplot},title=`Rayleigh and Planck vs x`);
The plots agree really well for small x ( hf << kT)
Back
to Modern Physics Resource Page
