Introduction to the Maths Lab - Part 2

Introduction

In this session we will concentrate on the graphical capabilities of Mathematica; we will also show you how data can be imported into Mathematica from other applications and how it can be exported out of Mathematica.

1. Two Dimensional Plotting

A standard task in the physical laboratory is the plotting and analysis of two dimensional experimental data. In many cases this data will be held as a sequence of numbers in a file external to Mathematica (in fact, all computer-driven lab instruments can give you results directly in this form). The numbers can be stored in a variety of formats; for this exercise the file is simply a text file holding pairs of values of time and intensity separated by a space, each pair appearing on a new line. 

To see what is inside any file you simply type in two exclamation marks followed by the path to the file. The example from our Macintosh computer was:

! ! "Hard disk:Desktop Folder:expt1.dat"

The path will be different for you, and the way the path is expressed will look different on a Windows PC. So to find expt1.dat on your computer, first start Mathematica, open up a new notebook, type in two exclamation marks and then go to  the Input menu and select Get File Path ...  to locate the data file expt1.dat on your hard disk. When you click on Open the appropriate file path will appear on the input line. Press <shift-return> now and the contents of the file will appear.

To read the data into Mathematica we use the Import command. The following command opens up the disk file and reads the pairs of numbers into a list that we have called data. Try this command, using your own file path, and compare it to the contents you saw above.

data = Import["Hard disk:Desktop Folder:expt1.dat"]

Import is a sophisticated command with many options (see the help system for the details); you can tell it explicitly what kind of data format you have got, or, as here, you can allow it to make an assumption - when it sees the ".dat" in the filename it assumes that the data is going to be in a plain text "tabular" format.

We can plot the data list using the following command.

ListPlot[data]

The data is from a kinetic experiment measuring current in μA as a function of time in s, so we should label the plot appropriately. We also choose to join up the data points, colour the resulting line, put the plot into a frame, add a title and change the font.

\(plot = 
      ListPlot[data, Frame -> True, 
        FrameLabel -> {"t/s", "I/μA"}, 
RotateLabel -> False, 
        PlotLabel -> "Kinetic Data", 
        DefaultFont -> {"Helvetica", 9}, 

        PlotStyle -> RGBColor[1, 0, 0]]\)

It is anticipated that the data should follow first order kinetics; we can check this by plotting the (natural) logarithm of the current against the time. To do this we first need to load up an add-on package which allows all sorts of graphical manipulations

Needs[ "Graphics`Graphics`"]

and then plot the intensity data on logarithmic scales.

\(logplot = 
      LogListPlot[data, Frame -> True, 

        FrameLabel -> {"t/s", "I/μA"}, RotateLabel -> False, 
        PlotLabel -> "Kinetic Data", 

        DefaultFont -> {"Helvetica", 9}, 
        PlotStyle -> {RGBColor[1, 0, 0]}]\)

Exercises and Explorations 1

1: Logging the data

You can estimate the gradient and intercept from this plot and use these values in the model of the current as a function of the time, t.

current[t_] := (Subscript[I, 0]) (E^(-λ t))

Notice that you have to define the function via the "colon-equals" sign and each variable must be declared (on the left hand side) by being followed by an underscore. However, when you want to evaluate the function, for example at t = 10, you simply type current[10]; there is no underscore. Below we show how the function has been used to plot the model for rather poor guesses at the kinetic parameters.

\(logmodel = 
      LogPlot[current[t]\  /. {I\_0 -> \ 1.5, λ\  -> \  .3}, {t, 0, 
          4}, Frame -> True, 
FrameLabel -> {"t/s", "I/μA"}, 
        RotateLabel -> False, 
DefaultFont -> {"Helvetica", 9}, 
        PlotStyle -> {RGBColor[0, 0, 1]}]\)

Here we have introduced the idea of replacement rules when we asked for current[t] to be evaluated for {(Subscript[I, 0]) ->  1.5, λ  ->   .3}. (Using replacement rules means that the values for (Subscript[I, 0]) and λ are only "remembered" during the calculation; we haven't assigned permanent values to them.) Adjust the parameters until you find the best fit to the experimental data; checking your guesses by superimposing data and model:

Show[logplot, logmodel]

2: Exporting data from Mathematica

There may be occasions when you need to "export" calculations that you have done with Mathematica to other programs - for example Excel, CricketGraph, CurveFit or SigmaPlot. There are a number of ways of doing this, but the most straightforward is to generate a simple data file which such programs can read.

In this way, the model function that you worked out in question 1 can be exported as a table of values for t and current(t). Try the following command, substituting your calculated values for (Subscript[I, 0]) and λ, and an appropriate file path for your computer.

\(Export["Hard disk:Desktop Folder:mymodel.dat", 
Table[{t, 
          current[t] /. {I\_0 -> \ 1.5, λ\  -> \  .3}}, {t, 0.0, 4.0, 
          0.2}]
]\)

And to see the contents:

!! "Hard disk:Desktop Folder:mymodel.dat"

This model has a simple decaying graph. When a graph is more "wiggly" you'll probably need to increase the number of data points by setting the step value (0.2 above) appropriately.

If you are familiar with the plotting features of Excel, or  a general plotting program, try recreating with it the superimposed plot of the data points and model curve, based on the files "expt1.dat" and "mymodel.dat".

Import and Export are powerful commands, and Mathematica has many others as well for handling the input and output of data, text and graphics. A good place to start hunting for information about this is Section 2.11 of The Mathematica Book. (In the Help system, click on the button labelled "The Mathematica Book" and type "2.11" into the Go To box; the same information can be found on the web at http://documents.wolfram.com/v4/ .)

It would be a good idea at this point to save all your work from Section 1, and begin a fresh notebook for Section 2.

2. Surface Plots

A surface plot arises when we have a function z=f(x,y); that is, the height z of the surface is a function of position in the x-y plane.  A good chemical example of this arises when we consider the quantum mechanical wave function of a particle confined to a rectangular planar surface of side (Subscript[l, 1])in the x direction and (Subscript[l, 2])in the y direction. (This is a reasonable model of the wave function for a molecule adsorbed on to a crystal surface.) The wave function, ψ(x,y), for the particle is given by 

ψ(x,y)=(2/(Sqrt[(Subscript[l, 1]) (Subscript[l, 2])]))sin((((Subscript[n, 1]) πx)/(Subscript[l, 1])))sin((((Subscript[n, 2]) πy)/(Subscript[l, 2])))   (Subscript[n, 1])= 1, 2, ... and (Subscript[n, 2])= 1, 2, ...

where (Subscript[n, 1])and (Subscript[n, 2]) are the quantum numbers and the energy of each state is given by

E = \(TraditionalForm\`\(h\^2\/\(8  m\)\) \((n\_1\%2\/l\_1\%2 + 
            n\_2\%2\/l\_2\%2)\)\)

Ignoring the scaling constants, you can plot the three dimensional surfaces using the following code. We start by looking at the wavefunction for a particle confined to a square surface.

ψ[x_, y_] := 
      Sin[(((Subscript[n, 1]) π x)/(Subscript[l, 1]))] 
        Sin[(((Subscript[n, 2]) π y)/(Subscript[l, 2]))]

Plot3D[ψ[x, y]  /. {(Subscript[n, 1]) -> 1, (Subscript[n, 2]) -> 2, (Subscript[l, 1]) -> 1, 
          (Subscript[l, 2]) -> 1}, {x, 0, 1}, {y, 0, 1}]

As with two dimensional plots we can change the appearance of the surface plot in a large number of ways. We increase the resolution of the plot, label the axes, change the viewpoint, the length of the edges of the box, within which we plot ψ, and colour code the height of the function. The latter makes high points red and low points orange. The view point was set using the 3D ViewPoint Selector ... under the Input pull down.

\(Plot3D[ψ[x, y]\  /. {n\_1 -> 1, n\_2 -> 2, l\_1 -> 1, 
          l\_2 -> 1}, {x, 0, 1}, {y, 0, 1}, 
PlotPoints -> 30, 
      ColorFunction -> Hue, AxesLabel -> {"x", "y", "ψ"}, 

      ViewPoint -> {5.375, \ 2.227, \ 2.965}, 
      BoxRatios -> {1, 1, 1}]\)

Exercises and Explorations  2

1: Plotting (ψ^2)

The probability of finding the particle at a position (x, y) is equal to (ψ^2). Plot (ψ^2) for the quantum numbers used above. (N.B. you should use the normalisation constant given above.)

This would be a good time to try out the "RealTime3D" graphics feature that we introduced in Maths Lab 1. Note, however, that in RealTime3D mode, axes and axes labels (and other things) will not appear.

2: Contours

Plot the function from exercise 1 as a contour plot, making sure the axes are labelled.

3: Degeneracy

Show that the wave functions for (Subscript[n, 1])= 1, (Subscript[n, 2])= 2 and for (Subscript[n, 1])= 2, (Subscript[n, 2])= 1 are degenerate when (Subscript[l, 1])= (Subscript[l, 2]). That is, you should show that the energy for these two states is the same, but the wave functions are different. What is the relationship between these degenerate wavefunctions?

4: Accidental degeneracy

Obvious symmetries do not always exist between degenerate states. Show that when (Subscript[l, 2]) = 2  (Subscript[l, 1]) there is degeneracy between states (Subscript[n, 1 ])= 1, (Subscript[n, 2]) = 4 and  (Subscript[n, 1]) = 2, (Subscript[n, 2]) = 2, but that a simple symmetry operation does not connect the wavefunctions. If you have time, find another non-symmetrical pair of degenerate states for (Subscript[l, 2])= 2(Subscript[l, 1]). N.B. since (Subscript[l, 1])≠ (Subscript[l, 2])  now you will have to change the range over which you are plotting.


We recommend that you now save your work  from Section 2, quit and restart Mathematica, and begin a fresh  Notebook to work on Section 3.

3. Three Dimensional  Surface Representations

The surface plots of Section 2 are single valued functions of x and y, i.e. there is only one value of z for each position (x, y). Frequently we want to make three dimensional surface representations where the surface is a function of x, y and z as well. A good chemical example of this are the molecular orbitals of hydrogen.

In this case the spherical symmetry of the problem indicates that we will simplify the mathematical description if we use the spherical coordinates r, θ and ϕ. Here r is the radial distance from the origin, θ is the angle with respect to the z-axis and ϕ is the angle in the x-y plane with respect to the x-axis. (Hence x = r sinθ cosϕ, y = r sinθ sinϕ and z = r cosθ.) In these coordinates the hydrogenic wave function, Ψ, can be expressed as

Ψ(r,θ,ϕ) = R(r) Θ(θ) Φ(ϕ)

In other words the wave function is the product of separate functions of r, θ and ϕ. The angular parts, Θ(θ) Φ(ϕ), are functions called spherical harmonics and are what give us the hydrogenic molecular orbitals (henceforth MO). The MOs and their associated spherical harmonics are given below.

(OverBar[(UnderBar[({{MO, Spherical Harmonic}})])])
(UnderBar[({{s, Constant}, {(Subscript[p, z]), cosθ}, {(Subscript[p, y]), sinθ  sinϕ}, {(Subscript[p, x]), -}, {(Subscript[d, (z^2)]), 3 (cos^2) θ  -  1}, {(Subscript[d, yz]), cosθ sinθ  sinϕ}, {(Subscript[d, zx]), -}, {(Subscript[d, (x^2) - (y^2)]), \(\(\ \)\(\(sin\^2\) θ\ cos2ϕ\)\)}, {(Subscript[d, xy]), -}})])

We have created the three dimensional renditions of the spherical harmonics by loading  the following package

Needs["Graphics`ParametricPlot3D`"]

which contains code to plot surfaces in three dimensions using spherical coordinates. Below is the code to plot the (Subscript[p, z])orbital. The radial distance from the origin, r, is being plotted as a function of θ and ϕ. In this case according to the MO table r ∝ cosθ, however, because of the way in which SphericalPlot3D does its computations, it is necessary to make r ∝ |cosθ|. This is true for plotting all the spherical harmonics using this function. Opening up the Help file on ParametricPlot3D you may be able to work out why this is.

SphericalPlot3D[
      Abs[Cos[θ]], {θ, 0, π},  {ϕ, 
        0,   2 π}]

The code which follows labels things neatly, and increases the number of points being used to plot the surface to give it a smoother appearance.

\(SphericalPlot3D[
      Abs[Cos[θ]], {θ, 0, π}, \ {ϕ, 0, \ \ 2\ π}, 
      PlotRange -> All, 
	AxesLabel -> {"x", "y", "z"}, 
      PlotLabel\  -> \*""\!\(p\_z\)-orbital"", PlotPoints -> 25, 
	
      DefaultFont -> {"Helvetica", 8}, 
	
      AxesEdge -> {{\(-1\), \(-1\)}, {1, \(-1\)}, {\(-1\), \(-1\)}}]\)

Exercises and Explorations 3

Making molecular orbitals

Replicate the s, p and d orbitals shown above, making sure everything is correctly labelled. You will have to figure out how to produce (Subscript[p, x]), (Subscript[d, zx]), and (Subscript[d, xy]).