Mathematica Notebook to Solve Problem 1.17
             Physics 176    Prof. Greenside    1/12/2011

Overview

The following is a Mathematica notebook that illustrates, in the context of solving Problem 1.17 on page 9 of the Schroeder text, some of the Mathematica commands that you will find useful in Physics 176 such as plotting a table of data, plotting functions of interest, and fitting a function to model some data and then comparing the fitted function with the data. This notebook will also give you some guidance of how to organize Mathematica calculations in a format suitable for a homework assignment.

When using Mathematica for a homework assignment, I recommend that you write out algebraic statements or derivations that are more than a line long on paper (I recommend using pen, not pencil, since ink does not smear) rather than try to use Mathematica for derivations, but do feel free to use Mathematica for simple calculations like those in the first homework assignment and of course for generating various plots. You should then print out and attach your Mathematica notebook to the back of the written part of your homework assignment so the TA and I can see how you obtained your answers (and to help you find a mistake if you got an incorrect answer). It is possible to typeset your homework assignments rather beautifully with Mathematica (in fact, people have typeset entire books with Mathematica) but it is also time consuming to do so and generally it will be faster to just write things out on paper.

Mathematica has two software components, a notebook interface (which you are using right now) in which you type and execute commands and view results, and a kernel which receives commands from the notebook, executes them, and then returns output to be displayed by the notebook. Most of the time, people use the notebook interface to run commands interactively, but it is possible to submit commands to the kernel directly via a text file, for example for a long calculation that might last several days. (There won't be any of those in 176!)

The Mathematica notebook is organized into vertical hierarchical units called "cells". To create a cell, move the cursor anywhere into a notebook window (typically at the bottom of an existing cell or between cells) and observe when the cursor changes from a vertical bar to a horizontal line segment, at which point you can click to indicate that you want to create a cell at that location. Once you have indicated a cell location by a mouse click, choose what kind of information you want to fill the cell with by selecting the menus Format/Style (i.e., select the menu Format, then select the sub menu Style) and then choose a particular cell style. It is also possible to split a particular cell into two pieces at the location of the blinking cursor (select the menus Cell/Divide Cell). This allows you to introduce a new cell at that location, e.g., insert some Mathematica code to execute in the middle of a text cell.

For homework assignments, you will likely use just four kinds of cells: Section, Subsection, Text, and Code. A Section cell with a single line of text is a good way to give a title to the entire notebook, say "Assignment 2 for Physics 176" with your name and date. (The first cell at the top of this notebook that says "Mathematica Notebook to Solve..." is a Section cell.)  A Subsection cell is a good way to provide a short title or one-line description for a sizeable unit such as the solution to a problem or part of a problem, say part (a) of problem 1.17; the word "Overview" above is in a Subsection cell. A Text cell is a good way to write descriptive text before and after some Mathematica code. (This cell which you are currently reading is a Text cell).  A Code cell is a good way to write and display your Mathematica code that will get executed. You will see Code cells further below in this notebook because they are shaded, gray for the default notebook and other colors for other styles.

Each cell is identified by one or more blue vertical brackets on the right side of the notebook. The innermost (left most) of the possibly several brackets denotes the vertical range of the current cell. If you double-click carefully on the second innermost bracket (second from the most left), you will collapse the cell to a single line, leaving a descriptive title. The bracket then becomes much shorter and has a downward arrow at one end to indicate that the cell is collapsed. Double-clicking on the collapsed bracket expands the contents of the cell again to make it fully visible. (Try collapsing this cell and expanding it again.) Collapsing and expanding cells is a convenient way to manage information in a long Mathematica notebook like this one.

You can delete a cell by clicking on the innermost bracket bounding that cell and by then hitting the Delete key. It is often useful to delete an output cell (one created by Mathematica after you executed some command) if you want to just show the text and executable code or if you execute some commands just to explore some results and then want to clean up afterwards. For example, after running the commands in this notebook to make sure that they worked correctly, I deleted the output cells so that you could see how the output is produced as you execute various commands. To delete all output cells, select the menu option Cell/Delete All Output.

Mathematica is a complicated piece of software and it can take many months to become familiar with its capabilities and foibles. For 176, you will only need to use some rather elementary commands and syntax but the more you learn, the greater the variety of scientific and mathematical problems you can solve. I encourage you to play with Mathematica on your own and to use it in other courses or for your own interests. To make the quickest progress, plagiarize enthusiastically by grabbing segments of Mathematica code from this notebook, from the Mathematica documentation, from the Internet, and from your own notebooks.

To continue, collapse this cell (double click on the second left most bracket on the right) to get it out of the way. Then expand the next cell (double click on the bracket that has a downward arrow head) to read it.

Typesetting in Mathematica (can skip on first read but useful)

There are two good tutorials on the website http://www.wolfram.com/broadcast/screencasts/, "Using the Mathematica Basic Math Assistant Palette" and "Using the Mathematica Writing Assistant Palette" that can help you get started in terms of typesetting notebooks and mathematical expressions. In a nutshell, bring up the Writing Assistant Palette by selecting the menus Palettes/Writing Assistant and then experiment by clicking on various buttons. However, the following comments should provide you with enough knowledge to produce rather nicely typeset notebooks.

Mathematical symbols involving English letters like B and t  are typeset by simply typing the usual letters in an italic font. To switch to an italic font, simply type "control-i" (hold the control key down on your keyboard while striking the letter i), type one or more letters, then type control-i a second time to leave italic mode. Or you can just select a letter or phrase with the mouse and type "control-i" to make the entire selection italic. You should try these two approaches right now with a simple example, for example try typesetting Newton's famous equation F=ma. 

Note: to play with Mathematica without losing your place in this notebook or without altering this notebook, open up a second Mathematica notebook window by using your mouse to go to the File menu then select New and then select Notebook. You can use this new notebook window as a convenient scratchpad.

Symbols involving non-English characters like the Greek lower case pi "π" or the multiplication symbol "×" can be inserted by clicking on the symbol as its appears in one of the palettes or by using a keyboard shortcut, which is revealed to you if you hold your mouse cursor over a symbol in the palette for about a second (the shortcut appears in the lower left corner of the palette). Thus you can type the Greek letter π by typing the three successive characters "Esc" "p" "Esc" (where Esc means hitting the Escape key on your keyboard), and you can type the multiplication symbol × by typing the three successive characters "Esc", "*" (that is an asterisk, on the number 8 key) ,"Esc", and you can type a plus-or-minus symbol ± by typing the four symbols "Esc", "+", "-", "Esc", and so on. Try typing these three symbols (π, ×, and ±) now for yourself.

Expressions involving powers and ratios like 1.01 × (10^5) N/(m^2) can also be typeset  by either using the mouse with a palette or by using keyboard shortcuts. (The latter are definitely faster but take longer to learn). To use the mouse, select the Mathematica menu "Palettes/Other/Basic Typesetting". This pops up a Basic Typesetting window that shows various mathematical symbols like the mathematical times symbol × (which you can find displayed under the entry "Operators" in the "Basic Typesetting" window, click on the right-pointing triangle to reveal the list of operators) that you can then select by just clicking on the symbol. A keyboard shortcut for that symbol again displays in the lower left corner of the "Basic Typesetting" window when you place your mouse cursor on the symbol. (A quick test: use your mouse and the palette to find the keyboard shortcut for the integration symbol ∫  and verify that the shortcut works.)

The Basic Typesetting menu also provides an easy way to typeset combinations of symbols such as powers and ratios by selecting what is called a template, an example of which is the Superscript template (■^□) which is the first entry on the 8th line of the "Basic Typesetting" window. For example, to typeset the expression (((α+β))^γ), I pressed the keys "(", "Esc", "a", "Esc", "+", "Esc", "b", "Esc", and ")" to produce the expression (α+β). I then selected that expression with the mouse (dragged the cursor through the expression to highlight it), then clicked on the Superscript template (■^□) in the palette, then typed the characters "Esc", "g", "Esc", "control-space" to put a gamma in the exponent. The character "control-space" is needed to jump out of a given template. Try this right now, to make sure you can also typeset this expression.

Note that you can type the tab  key to jump from one slot to the next in these templates. You should also type "control-space" (hold the control key down while typing space) to jump out of a template to the next point in your sentence.

Mathematica also allows one to typeset powers and ratios easily using keyboard commands. For example, if you type a symbol, then type the "raise-to-a-power" command "control-^" (hold down the control key while typing the carat symbol ^), a template will pop up with a slot for the exponent which you can then type in and then type control-space to exit the template. The command "control-/" will give you a fraction template, and other typesetting commands can be found by looking up "typesetting" in the Help/Documentation Center window. 

As a quick exercise, see if you can type set the following expression using the keyboard commands control-^ (raise to power), control-/ (form ratio), control-space (leave a mathematical template), and Esc p Esc for pi. 
                                                                         (((a^2) + (b^2))/((c^2) - (d^2))) - ((((π/n)))^(1/2)) .
You will have to select the entire numerator (a^2)+(b^2) with your mouse before typing control-/ to place that expression in the numerator of a ratio.  Don't forget to type control-space to exit a given template like a ratio so that you can type in other parts of the expression.   Hint: typeset ((((π/n)))^(1/2)) by typing these keys in succession: "(", "Esc", "p", "Esc", "control-/", "n",  "control-space", ")", "control-^", 1", "/",  "2", "control-space".

To continue, collapse this cell (double-click on the second left-most bracket on the right) and expand the next cell (double click on the short downward-pointing arrow on the right).

Solution to Problem 1.17(a)

1.17(a):  The goal of Problem 1.17 is to relate the van der Waals equation (a useful empirical generalization of the ideal gas law PV =nRT) for nitrogen gas  to a systematic extension of the ideal gas law called the virial expansion which holds broadly for low-density gases. The problem is not important for the course but does provide a good opportunity to show you how to use Mathematica. (Don't worry for now what is the van der Waals equation or what is a virial expansion.)

Let us first answer part (a) of Problem 1.17: "For each temperature in the table, compute the second term in the virial equation B(T) /(V /n)."  With only six rows of data in the table on page 9, it is straightforward to just calculate the necessary six numbers using a hand calculator, but we will do it the longer way using Mathematica  since what you learn will allow the efficient visualization and analysis of data and mathematical expressions many times later in the course.

We observe that the second virial coefficient B(T) and the volume per moleV /n have the same physical units so the ratio B(T) /(V /n) is a dimensionless number that can be compared directly to the number 1, to see if the ratio is a small correction or not to the ideal gas law. We can obtain the value of V /n without calculating the volume V or the number of moles n by rewriting the ideal gas law PV = nRT  in the formV /n = RT / P, where the temperature values are given in the table and we are told in the problem that the pressure P = atmospheric pressure = 1.01 × (10^5) N/(m^2). The universal gas constant R has the value 8.31 J/(mol·K). (I got these values for P and R from the inside front cover of the Schroeder book, and we do not need more than three significant digits for our calculations since the B values in the table have only three significant digits.) Since B has units of (cm^3) / mol, we need to convert its values to SI units of (m^3)/mol when calculating the desired ratio, i.e. we need to multiply each table entry for B by  ((((10^(-2)) m/cm))^3)   =(10^(-6 ))(m^3)/(cm^3).

Let's start by looking visually at the table of temperature T  versus second virial coefficient B on page 9; it is always a good idea to get a sense of what the data look like. We begin by storing all the rows of the table in a Mathematica variable which I call "tableTB". Variable names can be arbitrarily long but, until you have more experience, you should avoid choosing variable names that begin with a capital letter since Mathematica intrinsic functions all start with a capital letter. You should also not use an underscore in a variable name---don't write table_T_B as one would in C++ and in some other computer languages---because underscore has a special meaning in Mathematica for matching symbolic expressions. For those who know some programming: you don't need to allocate memory for variables nor do you have to declare the type of a variable (say character or double), Mathematica figures this out for you automatically by context.

One of the fundamental ways to organize information in Mathematica is by something called a list, which is similar to a mathematical  set: it is simply a collection of objects separated by commas between two braces {...}. The objects in a list can be anything including numbers, symbolic expressions, graphs, sounds, digital images, and other lists. 

A rectangular table of numbers is represented in Mathematica as a list of lists. Here we will make each row of the table a list of two elements, e.g., the first row will become the list {100., -160.}, the second row will be the list {200, 35.0}, and so on. (For scientific or engineering applications, always insert a decimal point for a number since Mathematica will otherwise treat the number as an exact integer and might use arbitrarily many digits when manipulating an integer or ratio of integers, leading to time consuming calculations.) We then create a list of the row lists to get our table and store the information in the variable tableTB using the equal signs "=" which is an operator in Mathematica that means "store the stuff on the right side in the variable on the left side". Here is the corresponding Mathematica code in a cell of type "Code":

tableTB = {                (* version 1: table of temp versus B in SI units *)
  { 100, -160.0 (10^(-6)) }, 
  { 200,  -35.0 (10^(-6)) },    (* note the conversion factors from (cm^3)/mol to (m^3)/mol *)
  { 300,   -4.2 (10^(-6)) }, 
  { 400,    9.0 (10^(-6)) },
  { 500,   16.9 (10^(-6)) },
  { 600,   21.3 (10^(-6)) }
}

I created the nicely shaded region of Mathematica code by clicking in the cell and selecting the menus Format/Style/Code.  

The information between the symbols (* and *) are comments that are  ignored by Mathematica and that provide a convenient way to explain what is going on in the Mathematica code. I have also multiplied each entry for B by the SI conversion factor (10^(-6)). I typeset this number by typing the six characters "1", "0", "control^", "-", "6", and "control-space". Try this yourself to see how it works.

In Mathematica , a space between two expressions implies multiplication of the expressions, just as is the case for symbolic expressions like (x+y) (u-v)  so you don't need to insert an explicit multiplication symbol like *  or × . Unfortunately, this assumption is also a common way of making a mistake in Mathematica since a symbol like ab to the eye looks mathematically like the symbol a times the symbol b while Mathematica actually treats this as a variable name called "ab", which is completely different.

Mathematica doesn't care how I type the numbers so long as the syntax is correct. For example, I could have typed instead:

tableTB={{100,-160.0 (10^(-6))},{200,-35.0 (10^(-6))},{300,-4.2 (10^(-6))},{400,9.0 (10^(-6))},   (* version 2 *)
		{500,16.9 (10^(-6))},{600,21.3 (10^(-6))}}

with no spaces and obtained the exact same result. Although this is more compact, its compactness is a bad thing since this second form is hard to read: it is more difficult to identify which rows I have typed (or even how many), whether I have  typed the numbers correctly, and whether I have matching braces.  In version 1, I have made the code more readable by indenting each row list, by using one line per row, by aligning the decimal points for the B entries, by adding extra spaces after opening braces and before closing braces, and by adding some comments to explain what is going on.

Please  take the time to use extra space and extra lines to format your Mathematica code so that it is easily readable by yourself and by others. This will save you many headaches when trying to find  errors that inevitably arise when using Mathematica. 

You can execute the above cell of code (either one) by clicking anywhere in the cell with your mouse and then by typing "Shift-Enter", i.e. while holding down the Shift key, hit the Enter key (also called the Return key). Please do this now to see what happens. If this is the first time you are executing code in this notebook, a warning will pop up saying "Do you want to automatically evaluate all the initialization cells in the notebook "problem-1.17-mma-tutorial.nb?" Here click "No", otherwise all the other pieces of code in this notebook will be executed, preventing you from running the code yourself for the first time.

Mathematica  does several things when you execute this code: it evaluates all expressions like the product of the conversion factors with the B values, it interprets the equal sign = to mean "store the evaluated list of lists in the variable tableTB"  and does so, and Mathematica automatically saves the output produced in an output variable Out[n] where n is the nth output expression created, and Mathematica displays the resulting output. The variables Out[n] provide a second way (besides a variable name) to retrieve and use earlier information. (Mathematica also uses a percent sign "%" as a symbolic way to refer to the most recent output, so you can use that output quickly for the next calculation.)

Note that Mathematica's default output is not user friendly, it produces a dense and difficult to read list. You can use the command TableForm as an easy way to get more readable output, for example execute the following cell:

TableForm[ tableTB ]    (* display variable tableTB in more friendly form *)

Note how Mathematica switches from fixed point notation (-0.00016) to floating point notation (9.`*^-6) when the magnitude of a number becomes too small or large. This default behavior can be changed to whatever you want, e.g., print all numbers in floating point form, but I won't take the time to explain how in this notebook.
 
The Mathematica command Dimensions provides a quick way to get a sense of what kind of data a particular variable contains, without looking at all the information associated with the variable (which might be many pages long). Thus execute the following cell:

Dimensions[ tableTB ]   (* show dimnensions of list TableTB *)

You should get the list {6, 2} which indicates that the variable tableTB is a list of six objects, each of which itself has dimension 2.  What do you think you will get for the value of Dimensions[ { a, {b,c}, {d,e, f} } ]? Try it and see if you can understand Mathematica's answer. (Suggestion: you can open up a second notebook window by selecting the menu option File/New/Notebook and type and execute various expressions in that second window so you don't mess up the contents of this window.)

Once you have a list like tableTB, you can access the list's elements by using a double bracket notation [[-]]. For example, the expression tableTB[[1]] with a single index 1 is the first element of the list, which means the list {100. , -160. (10^(-6))} corresponding to the first row. The expression "tableTB[[4,2]]" refers to the second component of the fourth member of the list, i.e., the value of B in the 4th row.  Execute each of the following cells to make sure you understand this notion:

tableTB[[3]]    (* third element of the list tableTB, i.e., the 3rd row of the table *)

tableTB[[3,1]]    (* first element, i.e., the temperature, for the third row*)

tableTB[[3,2]]    (* second element , i.e., the B value, for the the 3rd row of the table *)

{ tableTB[[2,2]], tableTB[[3,2]], tableTB[[4,2]] }   (* list of B values for rows 2, 3 and 4 of the table *)

Using the wrong number of brackets is a common and frustrating error in Mathematica, e.g. try executing "tableTB[3,1]" to see what happens. Mathematica uses single brackets for arguments of a function, e.g., Sin[π] denotes sin(π), the value of the sine function at the number π.  Mathematica tries to  interpret the expression tableTB[3,1] as the function tableTB evaluated with two arguments 3 and 1 (what we would write mathematically as f(3,1)) but can't make any progress since a list of lists is not a function. Mathematica therefore returns the expression without change, which can be confusing for beginners (and painful if the expression turns out to be many pages long).

Let us now plot the list of data using the Mathematica command ListLinePlot. Do this by executing the next cell:

g1 = ListLinePlot[ tableTB ]

This plot teaches us some science:   we see that the second virial coefficient B increases monotonically over the temperature range 100 to 600 K and starts to saturate (flatten out) at higher temperatures. Does it make sense scientifically that a correction to the ideal gas law might increase with temperature?

The Mathematica function ListLinePlot is smart enough to interpret each list of two numbers as a (x,y) coordinate pair, to draw line segments between successive points, to automatically determine the horizontal and vertical range of numbers, and to choose tick marks appropriately. All  details of the plot can be changed but I'll explain how to do that later in the semester (or you can explore this on your own, using the online Mathematica documentation).  ListLinePlot is widely used to get a quick visual sense of discrete data sets.

You can learn more about ListLinePlot by selecting the menus Help/Documentation Center and then type "ListLinePlot" into the Search box of the Documentation Center window that pops up. You can then read about the syntax and see various examples that illustrate various features. You can also execute, or edit and execute, any of the examples that appear in the Mathematica documentation, a useful feature. Executing a command like ?ListLinePlot in a notebook gives a brief summary of how to use the function, try this and see what happens.

Note that I saved the plot in a variable g1.  (The variable name is arbitrary, the first letter of g is just a hint to me that the variable contains a graphics object.) In Mathematica, anything can be saved in a variable and later recalled or manipulated. For example, to display the plot again later on, type "Show[ g1]".

It is easy to add labels to this plot, which would then make it self contained and so more suitable for handing in as part of a homework assignment. Just execute the following cell and mimic the syntax as needed for future cases:

g2 = ListLinePlot[                            
       tableTB ,                                (* data to plot *)
       PlotLabel -> "Table of B versus T" ,      (* add a label to the top of the plot *)
	   Frame -> True ,                           (* add a frame around the entire plot *)
       FrameLabel -> {"Temperature T (in K)", "Second Virial Coeff B (m^3)/mol"}   (* add labels to x and y axes *)
     ]

I typed the short right arrow → by typing successively the two characters "-" and ">" (minus sign followed by greater-than sign), which Mathematica automatically replaces with the little arrow. The three new lines  in the above code after tableTB are called "options" that represent optional parameters that modify the behavior of the function. You can learn what options are available for any function by typing "Options[function name]", e.g. execute "Options[ListLinePlot]".  You can use your typesetting tricks to include mathematical symbols and expressions in the plot labels, as well as ratios of units like (m^3)/mol.

Note again how I wrote the Mathematica code to be readable: I indented the entries relative to ListLinePlot, I used one line for each option,  I used spaces between operators like the arrow so the operator is clearly visible, I added comments to explain what is going on in a particular line, and I put the last bracket "]" on a separate line so I can easily see that it matches up with the opening bracket [ of ListPlot. Compare version 1 above to this equivalent but ugly alternative form:

g3 = ListLinePlot[tableTB,PlotLabel->"Table of B versus T",Frame->True,
     FrameLabel->{"Temperature T (K)","Second Virial Coeff B (m^3)/mol"}]

As you can guess, it is possible to change every detail of a plot, e.g., the size and fonts of the labels, the number and location of tick marks, the aspect ratio of the plotting box, and so on until you have a journal-quality figure. For 176, please label your plots but don't worry about making them beautiful. Spend your time instead thinking about the science contained in the plots.

OK, OK, let's finally answer part (a), namely compute the second term in the virial expansion B(T) / (V /n) to see how small it is compared to 1, which in turn tells us how big the correction is to the ideal gas law. We can do this by using the Mathematica command Table, an important and frequently used command that generates a new list, either from an existing list or from mathematical expressions.  (Please look up Table in the Mathematica Documentation Center to learn more.)

The following command calculates the ratio B/(V /n) = B/(RT /P) = B P /(RT) for each choice of T and B and saves the value in a list of pairs { T, BP/(RT) } so we can use ListLinePlot to visualize the resulting numbers. To understand the Mathematica code below, recall that the temperature T  in the ith row is tableTB[[i,1]], and the value of B in the ith row is tableTB[[i,2]].  Please execute the following cell and think about what the code means:

R = 8.31  ;          (* universal gas constant R in SI units *)
P = 1.01 (10^5) ;       (* atmospheric pressure in SI units N/(m^2) or pascals *)
Bvalues = Table[ 
  {                                          (* a list of two numbers *)
    tableTB[[i,1]] ,                            (* temperature T *)
    tableTB[[i,2]] P / ( R tableTB[[i,1]] )     (* calculate the expression BP/(RT) *)
  } ,    
  {i, Length[tableTB]}         (* iterator that varies i from 1 to the number of rows *)
]

Table here creates a list of two objects consisting of the ith temperature T and the corresponding ith value of B /(V /n) = BP /(RT). Table then takes the six lists and puts them in a single list that I call Bvalues, which you should think of as a 6×2 matrix or table.

The range of indices i is controlled by the expression {i, Length[tableTB]}, which is called an iterator in Mathematica and that is required by the Table function. The simplest iterator has the form {i,5} which tells the variable i to vary from 1 to 5.  A more complex iterator can start with other  values and use a specified spacing, e.g., {i, -4, 4, 2} will iterate over the successive indices -4, -2, 0, 2, 4 in steps of 2. 

(Note: for those of you familiar with Matlab, the iterator {i,5} would correspond to the Matlab notation 1:5, the iterator {i,-4,4,2} would be the notation -4:2:4.)

To see how iterators work execute the following code:

Table[ i , {i, 10, -6, -2} ]  (* create list of iterator index i *)

and make sure you understand why you got the answer. Then change the various values and re-execute the code to see what happens for the iterators {i,-6,10},  {i,10,-6},   {i,-6,10,2},  {i,10,-6,-2},  {i,0,1,0.2},  and  {i,-5π,7π,2π}. You can also experiment with generating other kinds of tables, e.g., type Table[ {i,  i^2, i^3}, {i,-2,2} ] or try a table of tables, Table[ Table[ i j , {i,  1,3}], {j,2,4} ]

Playing with small pieces of code like this is an essential way to understand and master the Mathematica syntax, you should be tinkering all the time when using Mathematica. In fact, I often keep two or more Mathematica notebooks open at any time (you can create a new notebook window by just selecting the menu File/New/Notebook) and use the second notebook as a scratchpad, to try out various small pieces of code. Do note that, when you open up new notebook windows, the symbols and their values that you have defined in one notebook are automatically carried over to the new notebook and vice versa. So if you change the value of some variable in one notebook, say execute "x=2", that will change the value of x in all open notebooks which can sometimes be a source of trouble. (There are ways to prevent this sharing of changes but that is an advanced topic for another day.)

 I created the new table Bvalues so that I could easily plot the answer using the following ListLinePlot (execute the following cell):

g4 = ListLinePlot[
		Bvalues ,
		PlotRange -> All ,         (* plot full vertical range of data *)
		Frame -> True ,
		FrameLabel -> {"Temperature T (K)", "B/(V/n) = BP/(RT)"}
]

We immediately see that the correction B /(V /n) due to the second virial coefficient is always small compared to one in this temperature range, ranging from the negative value -0.02 for T=100 to just above zero for the highest temperatures. (If I want, I  can see the smallest and largest  numerical values directly by evaluating Bvalues[[1,2]] and Bvalues[[6,2]], which give -0.0194 and 0.00043 respectively.)

Note that I used a new option "PlotRange -> All" to tell Mathematica to show all data in the vertical range. If you run this plot without this option by executing the following cell:

ListLinePlot[
	Bvalues ,
	Frame -> True ,
    FrameLabel -> { "Temperature T (K)", "B/(V/n) = BP/(RT)" }
]

you will see that Mathematica's plotting heuristics chooses poorly and throws out what seems to be an outlier for T=100 K. 

I recommend that you always use PlotRange -> All the first time around when using ListLinePlot or Plot (which plots functions rather than data) so you see all the data. You can then modify the plot accordingly, for example by using the Mathematica default or by using PlotRange to select a range of interest (the notation would be PlotRange -> {ymin,ymax}). You can not always trust the Mathematica graphical defaults.

This completes my discussion of 1.17(a). Collapse this cell by double clicking on the second most left of the cell brackets on the right, then expand the next cell to continue.

Solution to Problem 1.17(b)

1.17(b): I will let you solve part (b) on your own, this involves thinking but no calculation. As a hint, in your introductory physics or chemistry book, read about the long-range van der Waals force between any two neutral atoms or molecules, and review the consequences of the Pauli exclusion principle, which explains why atoms and molecules repel each other when they are so close that their inner electron clouds overlap.

Solution to Problem 1.17(c)

1.17(c):  This is the kind of problem that you should work out analytically on paper rather than get distracted by the machinery of Mathematica or by typesetting in Mathematica. I do outline the solution here and make a few comments, since we will need the answer for 1.17(d), for which we will return to using Mathematica.

The virial expansion given near the top of page 9 is an expansion in successive powers of the small quantity 1 /(V /n) = n /V  where n is the number of moles and V  is the volume of the gas. For n /V  to be small, the density of the gas must be small which is the same as saying that the average distance between molecules is large compared to the size of the molecules. (If you have already solved Problem 1.12 on page 8, as requested in Homework 1, you will know how to estimate the average molecular distance for an ideal gas.) 

To push the van der Waals equation on page 9 toward the form of the virial equation, we want to rewrite it in terms of the small quantity n /V where all the coefficients are expressed as functions of temperature.  By solving for P on the left side, you should be able to confirm that the van der Waals equation (page 9, near bottom, in the text) can be written in the form

P  = ((n R T)/(V - b n)) - a ((((n/V)))^2)  =  ((n R T)/V)( (1/(1 - b(n/V)))  -  ((a )/(R T))(n/V)) .

The second term inside the parentheses, -(a/RT)(n/V) is already of the form expected in the virial expansion, some function of temperature times some power of the small quantity n/V.  The first term inside the parentheses is of the form 1/(1-x)  where x=b(n /V) is a small quantity if n /V  is sufficiently small. We can then use a widely used approximation in physics and applied mathematics, the first few terms of the infinite geometric series expansion

(1/(1 - x)) = 1  +  x  +  (x^2) +  …  ,

to approximate 1/(1-x) to second order (up to second powers) of the small quantity x. (Although we could carry out the expansion to higher-order in x, Problem 1.17 only concerns the second and third virial coefficients B and C so we stop at second-order for that reason.) If x is small enough, the neglected higher-order terms can be ignored to any desired accuracy. 

Using Eq. (2) with x=b(n /V), please verify that the right side of Eq. (1) now becomes:

P  ≃  ((n R T)/V)( 1 + (b (n/V)) + (((b (n/V)))^2) -  ((a )/(R T))(n/V))  =  ((n R T)/V)( 1 + (b - (a/(R T)))(n/V)  +  (b^2) ((((n/V)))^2) ) ,

from which we can directly read off the second and third virial coefficients by comparing with the virial expansion:

B(T) = b - (a/(R T)) ,        C(T) = (b^2) .

Note how the expression for B has a form roughly consistent with the graph g2 that you produced in part (b): since the van der Waals parameters a and b are both positive, for sufficiently small temperatures T,  B(T) is indeed negative, and for sufficiently large temperatures, B becomes positive. There is a special temperature T=a /(bR) for which B(T) vanishes (known as the Boyle temperature) in which case the ideal gas equation becomes particularly accurate since the lowest-order correction in the virial expansion vanishes. We can't tell if the expression for C(T) makes any sense at this point without comparing with some data.

To practice approximating expressions to some order via the geometric series, approximate to third order in the small quantity ϵ the expression

((1 )/(1 - ϵ + 2(ϵ^2))) = 1  +  (?)ϵ  + (?) (ϵ^2) + (?(() ϵ)^3) + … ,

i.e., determine the values of the three coefficients indicated by question marks. Use Eq. (2) for this, don't use a Taylor series. For ϵ=0.1, you can use Mathematica to compare the value of your third-order approximation to the "exact" value 1.08696 of the left side.

I point out briefly, so you can explore on your own, that Mathematica knows how to calculate Taylor series symbolically of arbitrarily complicated mathematical expressions using the function Series. For example, execute the following two code cells:

Series[ (1/(1 - x)), {x,0,6} ]    (* calculate power series of 1/(1-x) around x=0 to 6th order *)

Series[ Sin[x], {x,0,8} ]    (* calculate power series of sin(x) around x=0 to 8th order *)

Series[ Exp[Cos[Sin[x]]], {x,0,8} ] (* series around x=0 to 5th order *)

You can easily compare the accuracy of a Taylor series to some order about some point of some function with the function itself as shown in the next cell. (Execute the next cell and experiment with different values of the order and expansion point). Here you see that the series of tan(x) about x=0 to order 7 accurately approximates tan(x) except where tan(x)  starts to diverge at x = ±π/2.

expansionPoint = 0. ;  (* value around which Taylor series is expanded *)
orderOfExpansion = 7 ;

tanSeries = Normal[   (* use Normal to convert series to polynomial that can be plotted *)
     Series[ Tan[x] , {x, expansionPoint, orderOfExpansion} ] 
]  

Plot[ { Tan[x], tanSeries }, {x, -2, 2} ]    (* plot Tan[x] and series approximation over the interval [-2,2] *)

Plot[ Tan[x] - tanSeries , {x, -1, 1}, PlotRange -> All ]  (* plot difference of function with series to examine accuracy *)

Solution to Problem 1.17(d) using the Manipulate command

We now use Eq. (4) of our answer in part (c) (virial coefficients B in terms of the van der Waals parameters a and b) to answer part (d), which is to "Plot a graph of the van der Waals prediction for B(T) , choosing a and b so as to approximately match the data given above for nitrogen." That is, we want to choose the van der Waals parameters a and b so that the expression for B(T) in Eq. (4) of our answer in 1.17(c) above "best" matches the B vs T table of data. 

To determine two parameters a and b from data, we need two independent conditions that lead to two equations that we can then solve for a and b. With six data points, it is generally not possible  to make a two-parameter curve go exactly through all the data points (unless we are lucky enough that the points all lie exactly on the curve we are trying to fit) so one needs to use scientific or statistical arguments to deduce two conditions. Fitting data is a complicated subtle subject that is best addressed in the context of specific problems. Rather than digress (we'll discuss some data fitting later in the course), let me show you how to use the Mathematica  command Manipulate to explore visually how parameters affect some mathematical expression. This allows you to intuitively (and often quickly) fit a function to some data, or to get some insight why it may not be possible to succeed with fitting a particular function.

Before doing any kind of mathematical analysis or computer analysis, it is useful to get a sense of the order of magnitude of the parameters of interest. We can get rough estimates of the parameters a and b by thinking a bit. For example, the expression B(T) = b - a /(RT) approaches the constant b in the limit T -> ∞ and we have seen from the above graph g2 that B(T) in the table has saturated (become roughly constant) at large temperatures. So we can guess that an estimate of b would be the large-T  limit of the available data, say 30 × (10^(-6)) (m^3)/mol (a bit bigger than the table value 21.3 × (10^(-6)) since the graph g2 shows that the curve is still increasing at T=600 K).

To get an estimate for a, we just need one additional condition, e.g, that the curve b  - a /(RT) pass through a given point of interest or perhaps that this curve pass through zero (for which the ideal gas law is especially accurate, since the lowest-order correction then vanishes) at about the same place that a linear interpolation of the data passes through zero.  Let's try the latter condition, not to get an accurate estimate but just to get a sense of how big a might be. The plot g2 shows that the data passes through zero around T  ≃ 330 K. (If you right-click on the plot, you get a menu that allows you to choose an option "Get Coordinates", which will then you give the x-y coordinates of any point that you click on with your mouse.) Now if function b - a/(RT) vanishes then a = b R T . We know the values of T ≃ 330 K,  b ≃ 30 × (10^(-6)), and R ≃  8.3  and so find that a  ≃ 0.06 J·(m^3)/ (mol^2) to one digit. For these values, let's compare B(T) with the data (please execute the following cell):

b = 30. (10^(-6)) ;             (* give symbol b a value *)
a = b R 330. ;             (* define a by condition for zero of B(T) at T=330 *)
bfn[ T_ ] :=  b - (a/(R T))    (* define function B(T) = b - a / (RT) ; note how the dummy variable T is green *)

g5 = Plot[ bfn[T], {T, 50, 650} ]  (* Plot function B(T) over the range 50 ≤ T ≤ 650 *)

g6 = ListPlot[       (* plot data points without connecting adjacent points with line segments *)
		tableTB,
		PlotStyle -> PointSize[ 0.02 ] , (* make data points a bit bigger than default to see more easily *)
		Frame -> True ,
		FrameLabel -> {"Temperature T (in K)", "Second Virial Coeff B (m^3)/mol"}
    ]

g7 = Show[ g5, g6 , PlotRange -> Automatic]  (* combine the two plots into a single plot to compare function with data *)

We see from the g7 graph that our choice for B(T) does OK for T ≥ 300 and indeed passes through zero at T = 330 K, but does not do well for lower temperatures, especially for the data point at T = 100 K. The graph does confirm that our values for a and b are not too crazy.

Some comments about the Mathematica code in the above code box: I stored a numerical value in the symbol b, then I let Mathematica calculate the value of a from the known values of b, T, and R, so a will update if we choose a different value of b. I defined my own mathematical function bfn to represent the mathematical formula b - a/(RT). For now, just mimic this notation without worrying how it works but quickly: any variable on the  left side of the "delayed assignment" operator := that has an underscore like T_ means a dummy variable, whose value will not be used until the function is evaluated elsewhere with a specific argument. I used a function ListPlot (not ListLinePlot) to plot the data points without connecting line segments so that it is easier to compare a guess for B(T) with the data. I used a new option PointSize to make the data points a bit bigger and so easier to see. Finally, when defining the variable g7, I used the Show command to combine two previous plots into a single plot. Mathematica knows how to choose one set of axes and tick marks so both plots can be displayed together in a reasonable way. Generating several plots and combining them with Show is a common procedure in Mathematica.

The following code uses a new command Manipulate to allow us to vary the values of a and b at will (using a mouse) over specified ranges of interest. Please look over the code until it makes sense, then execute the code and read the next text cell below.

Manipulate[
	g8 = Plot[        (* create but suppress plot B(T) by using semicolon at end of statement *)
			b - a / (R T) ,    (* this is the mathematical expression to plot *)
			{ T, 50, 650 } ,   (* this indicates what variable (T) to vary and its range, [50,650] *)
			PlotRange -> { {50,650}, {-0.00020, 0.00005} }  (* specify explicit x, y ranges so graph doesn't jump about *)
	] ;

	g9 = ListPlot[    (* create but suppress plot of data points by using semicolon *)
		tableTB,
		PlotStyle -> PointSize[ 0.02 ] , (* make points a bit bigger than default to see clearly *)
		Frame -> True ,
		FrameLabel -> {"Temperature T (in K)", "Second Virial Coeff B (m^3)/mol"}
	] ;

	Show[ g8, g9 ] , (* combine plots, do not suppress output with semicolon so can see *)

	{ a, 0.01, 0.20, 0.01, Appearance -> "Labeled" },   (* allow a to vary over range 0.01 to 0.15 in steps of 0.01 *)
	{ b, 20. (10^(-6 )), 100. (10^(-6)), 5. (10^(-6)) , Appearance -> "Labeled"} (* allow b to vary over range 20 (10^(-6)) to 1000 (10^(-6)) in steps of 5 (10^(-6)) *)
]

The idea of Manipulate is fairly straightforward: you have some mathematical expression, here b - a/(RT), that depends on some knobs, here the two parameters a and b. You then tell Manipulate over what range of values you would like to vary the knobs and in what step sizes, these are the last two lines of the Manipulate function. You can use any sequence of Mathematica code between the Manipulate command and the specification of the parameter values with one important qualification: the code can only display a single result, which is what ends up being displayed by Manipulate.

If you look at the above code carefully, you will see that the Plot command and the ListPlot commands end with semicolons ";". The semicolon tells Mathematica to execute the statement but not display the output. Since the Show command does not end with a semicolon, that is the only command displayed, and so you get a single plot for each choice of the variables a and b. As you become more experienced with Mathematica, you will likely start to use semicolons on most statements since most outputs are not interesting to look at and take up a lot of space.

The output of Manipulate is a box with two "sliders" that you can click and drag with your mouse to vary the values of the corresponding parameters a and b. Try this right now, varying a and b to see how the function b - a/(RT) changes its shape.

Homework exercise: give me your best choice of a and b with a print-out of your best plot

Here is the only part of this tutorial you need to hand in: by varying the parameters a and b in the above Manipulate command, do your best to match the curve visually to the data T versus B(T) given in the table.

Write down your "best" values of a and b in your homework assignment.

Print out the Manipulate window that shows the function a-b/(RT) and the data for your best choice of parameters and include that plot in your homework assignment. To do this, select the output cell that contains just the Manipulate window, i.e., click on the left-most cell bracket that just spans the output cell. Then select the menu File/Print Selection to print to your favorite printer. Note: please do not give me a print out of this entire Mathematica notebook, the TA and I have no interest in seeing so much information.

Summary of Mathematica functions and syntax mentioned

The above discussion used the following Mathematica functions: TableForm, Dimensions, Length, ListLinePlot, Table, Plot, ListPlot, Show, and Manipulate. As you have time, you should look up these functions in the Documentation Center under the Help menu to become more familiar with them. The functions ListLinePlot, Plot, Table, and Manipulate will be especially useful during the semester.

The above discussion also introduced you to some syntax: the list notation, the double-bracket notation a[[i]] to access an element of a list, the "=" or assignment operator, the storing of numbers, expressions or even graphs in some symbolic name, the idea of using options to modify the behavior of a function, how to define a function via the syntax f[x_] := x^2, and the concept of an iterator like {i,1,10,2} that generates a sequence of numbers that can be used by some function like Table.

To test your understanding, see if you can deduce the result of each command before executing the following statements:

Clear[ a, b, c, d, e, f, list, gcos, gsin]  (* make sure these variables have no values *)

list = { a, {b, c}, {d, e, f} }

Length[ list ]
Dimensions[ list ]
Dimensions[ list[[3]] ]

list[[1]]
list[[2]]
list[[3]]
list[[3,2]]

list[[1]]^list[[2,2]]
{ list[[1]], list[[2,1]], list[[2,2]] }

Table[ list[[i]], {i,3,1,-1} ]
Table[ list[[3,i]], {i,3} ]

gcos = Plot[ Cos[x], {x,0,2π} ]
gsin = Plot[ Sin[x], {x,0,2π} ]
Show[ gcos, gsin ]