Modern Physics (comments to: moloney@nextwork.rose-hulman.edu)  

MAPLE text and input lines 


PHYSICS RESOURCE PACKETS PROJECT            File: BARRIER.MWS
Rose-Hulman Institute of Technology         Authors: Dr. Kirtley & Perry Peters
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

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`);






This is a derivation of the compton wavelength formula.
It illustrates eliminating variables, dealing with analytic solutions.

It uses 'assume' to improve the chances of the derivation succeeding.

The 'solve' analytic solution of the two momentum equations gives the famous 'Root of ...' output.

To deal with this, one must use 'allvalues' to make sense of the output. This gives a solution for momentum, and angle phi. We select the momentum solution, and eliminate phi, making sure to select the positive momentum solution.

Download comptn2.mws  [Use 'Save File' and then open it in Maple]

PHYSICS RESOURCE PACKETS PROJECT            File: LEGENDRE.MWS Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4 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.MWS Rose-Hulman Institute of Technology         Authors: Perry Peters & Greg Williby http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4 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.MWS
Rose-Hulman Institute of Technology         Author: Perry Peters
http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4

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),r=0..scale*.75,theta=0..2*Pi,phi=0..Pi,
  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),r=0..3*scale,theta=0..2*Pi,phi=0..Pi,
  coords=spherical,grid=[20,20,30],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),r=0..2.3*scale,theta=0..2*Pi,phi=0..Pi,
  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: PSHELL.MWS Rose-Hulman Institute of Technology         Author: Perry Peters http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4 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.MWS Rose-Hulman Institute of Technology         Author: Dr. Mike Moloney http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4 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)
PHYSICS RESOURCE PACKETS PROJECT            File: SCHRO.MWS Rose-Hulman Institute of Technology         Author: Dr. Sudipa Kirtley http://www.rose-hulman.edu/~moloney         Software: Maple V Release 4 This is a trivial version of the solution to Schrodinger's equation of a particle moving freely in a one-dimensional box of length L. All units are in MKS system. restart; assume(m>0,L>0,hbar>0,E>0); _EnvAllSolutions:=true: Write down a simple wave equation, and find the solution to this differential equation. Notice that the term k^2 is equal to (2*m*E)/(hbar)^2, where E is the energy of the particle.  The wave function vanishes at the boundaries of the box, at x = 0 and at x = L. From this find the possible values of E, and the possible wave functions. k:=sqrt(2*m*E)/hbar; y:=rhs(dsolve({diff(f(x),x$2) + k^2*f(x) = 0,f(0)=0}, f(x))); eq1:=subs(x=L,y=0); E_levels:=solve(eq1, E); Note that the term _Z1 in the above output expression is an integer value included in the software; isolate that term. Also isolate the constant term in the y-expression, and plug it back in the wave function. Z:=op(1,op(3,E_levels)); ynew:=subs(op(1,lhs(eq1))=c,y); Now we need to find the value of the constant in front of the wave function from normalization conditions. arg:=simplify(subs(E=E_levels,ynew)); eq2:=int((arg)^2, x=0..L)=1; const:=solve(eq2,c); The exact wave function can be obtained now. y_exact:=subs(c=max(const),ynew); Insert here number of levels you want to examine. N:=5; The energy eigenvalues for the levels: for i from 1 to N do Energy.i:=subs(Z=i,E_levels); od; The corresponding wave functions for the different levels: for i from 1 to N do psi.i:=simplify(subs(E=Energy.i,y_exact)); od; Specify some constants such as the mass of the particle, the length of the box, and the value of hbar. Find the exact energy eigenvalues and the corresponding eigenfunctions, offset by the corresponding eigenvalue. values:=L=3E-10,hbar= 6.6260755e-34/(2*Pi),m=1E-30; assign(values[1]): for i from 1 to N do Energysub.i:=subs(values,Energy.i); y.i:= subs({values,E=Energysub.i},y_exact*hbar^2*Pi^2/(sqrt(2)*m*L^(3/2))+   Energysub.i); od; Plot the eigenfunctions associated with the corresponding eigenvalues. with(plots): g1:=plot({y.(1..N)},x=0..L,color=red): g2:=plot({Energysub.(1..N)},x=0..L,color=black,title=`Energy eigenvalue and the   corresponding wave function of a free particle in a box`): display({g1,g2});


Back to Modern Physics Resource Page