Archive for the ‘educational technology’ Category

h1

Notebook 2.1: Work in Reversible Isothermal Expansions

January 16, 2010

This short exercise gives students more practice using Maxima’s integration and comma operators.

  • Download and install the Euler Math Toolbox.
  • Cut and paste the code below into a plain text file.
  • Save the file with an .en file extension.
  • Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

--------------------------------------------------------snip here----------------------------------------------------
 Notebook 2.1: Work in Reversible Isothermal Expansions
% 
% In reversible expansions or compressions, the external pressure constantly changes so that 
% it remains exactly equal to the pressure of the gas. In this special case, we can substitute
% the pressure calculated from the gas's equation of state for the external pressure when calculating work.
% For example, the work done to compress n moles of an ideal gas from V1 to V2 under reversible
% conditions at temperature T is
% 
>: assume(V1>0, V2>0, V1>V2, V>0)
 
                        [V1 > 0, V2 > 0, V1 > V2, V > 0]
 
>: work_ideal: -integrate(n*R*T/V,V,V1,V2)
 
                            n R T (log(V1) - log(V2))
 
% If the "assume" had been omitted, Maxima would have asked a series of questions about 
% the signs of the volumes during the integration.
% 
% If we want to calculate the work for specific temperatures and volumes, we can use Maxima's
% comma operator. For example,
>: work_ideal, n=1, R=8.314, T=298, V1=1, V2=0.5, numer
 
                                1717.322046434265
 
% 
% Integrating -n*R*T/V with respect to V from V1 to V2 at constant T gives -n*R*T*ln(V2/V1), so we could
% also have calculated the work as -n*R*T*log(V2/V1). 
%  --> Notice that the Maxima function for natural logarithms is log(),  NOT ln()!)
%  --> The final ",numer" tells Maxima that we want a number as an answer. Maxima will sometimes 
%        just display logs of whole numbers without calculating them in a final result if we don't include this.
% ---------------------------------------------------------------------------------------------------------
% In the space below, create an expression called work_vdw that calculates the work for a reversible isothermal 
% expansion of a van der Waals gas as it expands from volume V1 to volume V2 at temperature T.
% 
>: assume(V>b, V1>b, V2>b)
 
                             [V > b, V1 > b, V2 > b]
 
>: work_vdw: -integrate(n*R*T/(V-b)-n^2*a/V^2,V,V1,V2)
 
                                       2                          2
                                    a n                        a n
              - n R T log(V2 - b) - ---- + n R T log(V1 - b) + ----
                                     V2                         V1
 
% ---------------------------------------------------------------------------------------------------------
% In the space below, use the comma operator to calculate the work for one mole of  xenon expanding from 1 m^3
% to a volume of 10 m^3 at a temperature of 298 K.
% 
>: work_vdw, a=5.125, b=1.06e-2, T=298, n=1, V1=1, V2=10, R=8.314,numer
 
                               - 5723.982679505956
 
% ---------------------------------------------------------------------------------------------------------
% In the space below, use the comma operator to calculate the work for one mole of  xenon expanding from 1 m^3
% to a volume of 10 m^3 at a temperature of 298 K IGNORING INTERMOLECULAR ATTRACTIONS.
% 
>: work_vdw, a=0, b=1.06e-2, T=298, n=1, V1=1, V2=10, R=8.314, numer
 
                               - 5728.595179505956
 
% ---------------------------------------------------------------------------------------------------------
% In the space below, use the comma operator to calculate the work for one mole of  xenon expanding from 1 m^3
% to a volume of 10 m^3 at a temperature of 298 K IGNORING INTERMOLECULAR REPULSIONS.
% 
>: work_vdw, a=5.125, b=0, T=298, n=1, V1=1, V2=10, R=8.314, numer
 
                               - 5700.207854019444
 
% ---------------------------------------------------------------------------------------------------------
% In the space below, use the comma operator to calculate the work for one mole of  xenon expanding from 1 m^3
% to a volume of 10 m^3 at a temperature of 298 K IGNORING ALL INTERMOLECULAR FORCES.
% 
>: work_vdw, a=0, b=0, T=298, n=1, V1=1, V2=10, R=8.314, numer
 
                               - 5704.820354019445
Advertisements
h1

Notebook 2.0: Introducing integration with Maxima

January 16, 2010

This very simple Euler notebook introduces integration with Maxima by asking students to set up simple integrals for calculating expansion work.

To use the notebook:

  • Download and install the Euler Math Toolbox.
  • Cut and paste the code below into a plain text file.
  • Save the file with an .en file extension.
  • Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

--------------------------------------------------------snip here----------------------------------------------------

 Notebook 2.0 - Integration with Maxima
% 
% When we're finding derivatives, we are finding the instantaneous rate
% of change for a function. For example, we can see how P changes with 
% V for an ideal gas by calculating the derivative dP/dV at constant T:
 
>: Pideal : R*T/V
>: dPdV_ideal : diff(Pideal,V)
 
                                        R T
                                      - ---
                                         2
                                        V
 
>: Pvdw : R*T/(V-b)-a/V^2
>: dPdV_vdw : diff(Pvdw, V)
 
                                 2 a     R T
                                 --- - --------
                                  3           2
                                 V     (V - b)
 
% 
% Often we want to calculate the change in a variable from its rate
% of change. For example, if we have dP/dV, we can reverse the process
% above to find P, by integrating the derivative with respect to
% V:
% 
>: Pideal : integrate(dPdV_ideal,V)
 
                                       R T
                                       ---
                                        V
 
>: Pvdw : integrate(dPdV_vdw,V)
 
                                    R T    a
                                   ----- - --
                                   V - b    2
                                           V
 
% 
% Integration undoes differentiation, and vice versa:
% 
>: integrate(diff(Pideal,V),V)
 
                                       R T
                                       ---
                                        V
 
>: diff(integrate(Pideal, V),V)
 
                                       R T
                                       ---
                                        V
 
% -------------------------------------------------------------------------------------
% In the space below, calculate the derivative of P for a Dieterici gas with respect
% to V (with T held constant. Call it dPdV_Dieterici.
% 
 
>: P_Dieterici : R*T*exp(-a/(R*T*V))/(V-b)
>: dPdV_Dieterici : diff(P_Dieterici,V)
 
% -------------------------------------------------------------------------------------
% In the space below, show that integrating dPdV_Dieterici gives back
% P for a Dieterici gas..
% 
>: integrate(dPdV_Dieterici,V)
% -------------------------------------------------------------------------------------
% We can use integration to calculate the work done during gas expansions, or the work
% required to compress a gas. This is done by integrating the external pressure (NOT the
% gas pressure) with respect to volume, over some volume range. 
% 
% In Maxima, you can integrate over a range by adding the beginning and end of the range as 
% third and fourth arguments for the integrate() function. For example, the work done to compress a gas from 4 m^3 to 2 m^3 under an external
% pressure of 101000 Pa is
 
>: work: -integrate(101000,V,4,2)
 
% When Pext is constant, the work is simply -Pext times the volume change, so we could have
% calculated the work more simply as
% 
>: work: -101000*(2-4)
% 
% Both calculations give work in joules, because we've used SI units (pascals and cubic meters)
% throughout the calculation. 
% 
% In the space below, calculate the work done (in joules) when a sample of 1.70 mol of argon gas
% expands from 20 L to 50 L against a constant external pressure equal to the final
% pressure of the gas. The temperature of the gas is 25 degrees C.
% 
>: n: 1.70
>: R: 8.314
>: T: 25+273.15
>: V2: 50*1e-3
>: V1: 20*1e-3
>: Pext: n*R*T/V2
>: work: -integrate(Pext,V,V1,V2)
% or, more simply,
>: work: -Pext*(V2-V1)
% -------------------------------------------------------------------------------------
% In the space below, calculate the expansion work done when 1 mole of water is electrolysed
% under a constant 1 atm external pressure at 298 K. Hint: Study example 2.1 in the text!
% 
% Solution: The work is equal to -Pext times the volume change. The initial volume of the water
% is 1 mol * 18 g/mol * 1 mL / 1 g * 1e-3 L/ 1 mL * 1 dm^3 / 1 L * 1e-3 m^3 / 1 dm^3
>:  V1: 1*18*1*1e-3*1*1e-3
% 1/2 mol of O2 and 1 mole of H2 is produced by the reaction. The final pressure of the gases
% will be equal to the external pressure. If the gases are assumed to be ideal, the final volume is
% n*R*T/Pext.
>: n: 1.5
>: R: 8.314
>: T: 298
>: Pext: 1 * 101325
>: V2: n*R*T/Pext
% and the work is
>: work: -integrate(Pext,V,V1,V2)
% or, equivalently, 
>: work: -Pext*(V2-V1)
% We could have simplified the solution by neglecting V1, since gases have about 1000 times
% the molar volume of liquids around room temperature and pressure. Then the work is 
% approximately
>: work: -Pext*V2
------------------------------------snip here------------------------------------
h1

Notebook 1.5: Critical Point Coordinates from Equations of State

January 16, 2010

This notebook uses Euler’s symbolic math capabilities to derive the critical point pressure, volume, temperature, and compression factor for a van der Waals gas. It then asks students to repeat the calculations for a Dieterici gas.

To use the notebook:

  • Download and install the Euler Math Toolbox.
  • Cut and paste the code below into a plain text file.
  • Save the file with an .en file extension.
  • Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

------------------------------------snip here-----------------------
 Notebook 1.5: Calculating Critical Point Coordinates from Equations of State
% 
% In this notebook, we'll use Euler's symbolic math capabilities to
% derive the critical point pressure, volume,
% temperature, and compression factor for a van der Waals gas. You'll
% then repeat the calculations yourself for a Dieterici gas.
% 
% We saw in the last notebook that Euler calls a separate package 
% called Maxima to do symbolic math.
% When you want numerical answers, use Euler; but if you want to 
% solve or manipulate equations, send your commands to
% Maxima by typing a colon (':') after the prompt.
% 
% Here is the van der Waals equation, in Maxima:
 
>: P : R*T/(V-b)-a/V^2
 
                                    R T    a
                                   ----- - --
                                   V - b    2
                                           V
 
% 
% Notice that Maxima uses a colon rather than an equals sign to set
% variables. Notice also that the result is symbolic, rather than 
% numerical, so you don't have to assign numbers to R, T, a, etc. as 
% you must from a regular Euler prompt.
% 
% Maxima doesn't simplify expressions automatically. You have to be
% explicit about how you want an expression 'simplified'. For example,
% you can remove parent
 
% Maxima can differentiate expressions with the diff(expr,x) function, 
% where expr is the expression and x is the variable you want to 
% differentiate. Let's calculate the first derivative of P with respect
% to V, and call it dPdV:
% 
>: dPdV : diff(P,V)
 
                                 2 a     R T
                                 --- - --------
                                  3           2
                                 V     (V - b)
 
% 
% Call the second derivative d2PdV2:
% 
>: d2PdV2 : diff(dPdV,V)
 
                                  2 R T     6 a
                                 -------- - ---
                                        3    4
                                 (V - b)    V
 
% 
% We could have calculated the second derivative more simply as
% 
>: diff(P,V,2)
 
                                  2 R T     6 a
                                 -------- - ---
                                        3    4
                                 (V - b)    V
 
% 
% where the third argument to diff is the order of the derivative.
% 
% Now let's calculate the temperature and volume of a van der Waals
% gas at the critical point, using Maxima's solve function.
% The first argument of solve is an array of equations to solve.
% The second argument is an array of unknowns to solve for.
% 
% 
>: solve([dPdV=0,d2PdV2=0],[T,V])
 
                                            8 a
                     [[T = 0, V = b], [T = ------, V = 3 b]]
                                           27 b R
 
% 
% We could also have typed 
% 
>: solve([diff(P,V)=0, diff(P,V,2)=0], [T, V])
% 
% for the same result. We can't put P into the array of unknowns,
% because we've assigned an expression to P already.
% 
% Notice that there are two solutions- one at absolute zero, where all
% of the volume is taken up by the molecules- and the "real" critical 
% point. Let's call the second solution "criticalpoint":
>: criticalpoint : %[2]
 
                                    8 a
                              [T = ------, V = 3 b]
                                   27 b R
 
% 
% The % alone means "the last thing calculated in Maxima". The last
% thing was an array, and we are selecting the second element with the 
% bracketed 2.
% 
% We can now use Maxima's "at" function to calculate the coordinates
% of the critical point. The at function's first argument is the 
% variable we want to solve for; the second argument is an array of
% equations that hold at the point where the solution exists.
% 
>: Pc : at(P,criticalpoint)
 
                                        a
                                      -----
                                          2
                                      27 b
 
>: Tc : at(T,criticalpoint)
 
                                      8 a
                                     ------
                                     27 b R
 
>: Vc : at(V,criticalpoint)
 
                                       3 b
 
% 
% Finally, we calculate Zc, the compressibility factor at the critical
% point:
% 
>: Zc : Pc*Vc/(R*Tc)
 
                                        3
                                        -
                                        8
 
% 
% -----------------------------------------------------------------
% In the space below, derive expressions for Pc, Vc, Tc, and Zc for
% a Dieterici gas. You can check your results against Table 1.7 in
% Atkins. 
% 
% We'll clear all the variables above using 'reset':
>reset;
>: P : R*T*exp(-a/(R*T*V))/(V-b)
 
                                            a
                                        - -----
                                          R T V
                                  R T %e
                                  -------------
                                      V - b
 
>: dPdV : diff(P,V)
 
                                   a               a
                               - -----         - -----
                                 R T V           R T V
                           a %e          R T %e
                           ----------- - -------------
                            2                     2
                           V  (V - b)      (V - b)
 
>: d2PdV2 : diff(dPdV, V)
 
                    a               a                a                 a
                - -----         - -----          - -----           - -----
                  R T V     2     R T V            R T V             R T V
          2 a %e           a  %e           2 a %e          2 R T %e
        - ------------- + -------------- - ------------- + ---------------
            3                  4             2        2              3
           V  (V - b)     R T V  (V - b)    V  (V - b)        (V - b)
 
>: solve([dPdV=0,d2PdV2=0],[T,V])
 
                                             a
                     [[T = 0, V = b], [T = -----, V = 2 b]]
                                           4 b R
 
>: criticalpoint : %[2]
 
                                     a
                              [T = -----, V = 2 b]
                                   4 b R
 
>: Pc : at(P,criticalpoint)
 
                                       - 2
                                     %e    a
                                     -------
                                         2
                                      4 b
 
>: Tc : at(T,criticalpoint)
 
                                        a
                                      -----
                                      4 b R
 
>: Vc : at(V,criticalpoint)
 
                                       2 b
 
>: Zc : Pc*Vc/(R*Tc)
 
                                         - 2
                                     2 %e
 
% 
% Is Zc for the Dieterici equation different than Zc for the van der
% Waals eqution? Which is closer to the experimental values of Zc given
% in your text?
% 
>3/8
                   0.375 
>2*%e^-2
         0.2706705664732 
% ANSWER: The experimental values in Table 1.5 range from 0.27 to about 
% 0.31; The Dieterici equation clearly gives more realistic results than
% the van der Waals equation.
------------------------------------snip here----------------
h1

Notebook 1.4: Trends in Boyle temperatures

January 15, 2010

All of those ugly, clunky-looking equations of state for real gases have to simplify to PV = nRT at very low pressure. It’s tempting to say that real gases behave like ideal gases under those conditions, but that isn’t quite right.

Even at pressures near zero, all of the properties of real gases won’t become ideal. In particular, any property that depends on the second (and higher) virial coefficient won’t necessarily converge to its ideal gas value as pressure drops. That convergence occurs only at the Boyle temperature, where the second virial coefficient is zero …as long as we assume all the higher order virial coefficients are zero, too. Usually that isn’t a bad assumption.

The Euler notebook posted below asks students to find the Boyle temperatures of noble gases from the van der Waals equation of state. They then critically examine trends in their calculated Boyle temperatures. The notebook introduces algebraic calculations with Maxima, including limits, differentiation, and equation solving.

To use the notebook:

  • Download and install the Euler Math Toolbox.
  • Cut and paste the code below into a plain text file.
  • Save the file with an .en file extension.
  • Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

—————————-snip here———————————————
 Notebook 1.4: Boyle Temperatures of Noble Gases
% 
% In this notebook, we'll use Euler's symbolic math capabilities to
% derive an expression for the Boyle temperature in terms of the 
% van der Waals constants a and b. You'll then use this expression
% to try to discover periodic trends in the Boyle temperature.
% 
% We can derive algebraic expressions in a notebook using Euler.
% Symbolic algebra and calculus can be done with a software package 
% called Maxima which is built into Euler. We can talk to Maxima
% directly by typing a colon (':') after the '>' prompt.
% For example, here is the van der Waals equation in Maxima:
% 
>: (P+a/V^2)*(V-b)=R*T
 
                              a
                             (-- + P) (V - b) = R T
                               2
                              V
 
% 
% If you forget to type : after the > prompt, the command gets
% sent to Euler, not to Maxima. Euler will give you an error 
% message (because it can't compute a result unless you assign numbers
% to each variable). 
% 
% It's useful to name equations, variables, and expressions so
% we don't have to type them over and over again. Maxima lets us
% type a name (say, 'vdw') followed by a colon to name an 
% expression: 
% 
>: vdw : (P+a/V^2)*(V-b)=R*T
 
                              a
                             (-- + P) (V - b) = R T
                               2
                              V
 
% 
% Use a colon to set a name, not an equal sign- you can see 
% from the previous line that '=' means equality, not assignment!
% 
% Now we can just type vdw when we want the van der Waals equation:
% 
>: vdw
 
                              a
                             (-- + P) (V - b) = R T
                               2
                              V
 
% 
% We can make substitutions in the equation using the at() function,
% which looks like at(expression,[conditions]). For example, to see
% that the van der Waals equation is the ideal gas 
% equation when a and b are zero, we could type
% 
>: at(vdw,[a=0,b=0])
 
                                    P V = R T
 
% 
% To solve an equation in Maxima, type solve(equation,x), where
% equation is the equation you want to solve and x is the variable
% to solve for. For example, to solve the equation 
% 2*x^2 + 3*x + 1 = 0 for x, you could type 
% 
>: solve(2*x^2 + 3*x + 1 = 0, x)
 
                                      1
                               [x = - -, x = - 1]
                                      2
 
% 
% Notice that the result is an array; there are two 
% solutions to this quadratic equation.  
% 
% To solve the van der Waals equation for P, you could type
% 
>: solve(vdw,P)
 
                                      2
                                 R T V  - a V + a b
                            [P = ------------------]
                                      3      2
                                     V  - b V
 
% 
% Notice that solve returns an array of solutions - even when 
% there is only one solution, as is the case here. We can 
% extract the % solution from the array by putting a subscript
% 1 (which is written as [1]) after solve:
% 
>: solve(vdw,P)[1]
 
                                      2
                                 R T V  - a V + a b
                             P = ------------------
                                      3      2
                                     V  - b V
 
% 
% The square brackets are gone. Suppose we're just interested
% in the right hand side of this equation. We can use 
% Maxima's rhs function to extract the right hand side:
% 
>: Pvdw : rhs(solve(vdw,P)[1])
 
                                    2
                               R T V  - a V + a b
                               ------------------
                                    3      2
                                   V  - b V
 
% 
% This doesn't look like equation 1.21a from the textbook-
% but it is! You can check to see that Pvdw really is 
% equivalent to R*T/(V-b) - a/V^2 by typing
% 
>: is(equal(Pvdw, R*T/(V-b) - a/V^2))
 
                                      true
 
% 
% Maxima reports and uses expressions in the form of rational polynomials
% (ratios of polynomials). As physical chemists, we would rather
% write expressions broken into separate terms that have some physical
% meaning.
% 
% Maxima can rearrange or simplify the expression, but we have to
% tell it what we want. Here, the partfrac() function can be used.
% The call looks like partfrac(expression,x) where the expression
% is expanded into fractions with powers of x in the denominator.
% Here we can use
% 
>: partfrac(Pvdw,V)
 
                                    R T    a
                                   ----- - --
                                   V - b    2
                                           V
 
% 
% which is exactly how the textbook writes the van der Waals
% equation.
% 
% We can also check to see if Pvdw is correct by seeing that Pvdw
% becomes equal to Pideal as V approaches infinity (which is typed
% as 'inf' in Maxima). To calculate the limit of y as x approaches
% a, type limit(y,x,a). If you enter
>: limit(Pvdw,V,inf)
 
 Is  R T  positive, negative, or zero?
 
>: positive
 
                                        0
 
% Maxima asks whether RT is positive, negative, or zero, and
% you must answer to get the limit. (When you answer, type : after
% the prompt, or the answer goes to Euler, not to Maxima) 
% Maxima then tells us that 
% Pvdw approaches 0 as V approaches infinity, which isn't the same
% as saying it approaches Pideal. Instead, we can define Pideal
% and see whether Pideal/Pvdw approaches 1 as V approaches infinity:
% 
>: Pideal : R*T/V
 
                                       R T
                                       ---
                                        V
 
>: limit(Pideal/Pvdw,V,inf)
 
                                        1
 
% We can also check derivatives of Pvdw this way. As Pvdw approaches
% Pideal, all of the derivatives of Pvdw should approach corresponding
% derivatives of Pideal. 
% 
% It's easy to calculate derivatives in Maxima. The derivative of
%%y with respect to x would be typed diff(y,x). To calculate the
% derivative of Pvdw with respect to V, we can type
>: dPdV_vdw : diff(Pvdw,V)
 
                                  2                2
                2 R T V - a   (3 V  - 2 b V) (R T V  - a V + a b)
                ----------- - -----------------------------------
                  3      2                 3      2 2
                 V  - b V                (V  - b V )
 
% In the space below, calculate the derivative of Pideal with respect
% to V and name it dPdV_ideal.
% 
>: dPdV_ideal : diff(Pideal,V)
 
                                        R T
                                      - ---
                                         2
                                        V
 
% In the space below, show that dP/dV for a van der Waals gas approaches
% dP/dV for an ideal gas as V approaches infinity.
% 
>: limit(dPdV_ideal/dPdV_vdw,V,inf)
 
                                        1
 
% 
% ----------------------------------------------------------------
% 
% Now let's derive an exptression for the Boyle temperature of a
% van der Waals gas. Our approach will be to:
%%1. Calculate dZ/d(1/V) for a van der Waals gas.
%%2. Calculate B as the limit of dZ/d(1/V) as V approaches
%%infinity.
%%3. Find the Boyle temperature as the the temperature where
%%B = 0.
% 
% In the space below, calculate the compressibility factor for a
% van der Waals gas and call it Z_vdw.
% 
>: Z_vdw : Pvdw*V/(R*T)
 
                                     2
                             V (R T V  - a V + a b)
                             ----------------------
                                      3      2
                                R T (V  - b V )
 
% In the space below, calculate the derivative of Z with respect
% to V, and call it dZdV.
% 
>: dZdV : diff(Z_vdw,V)
 
             2
        R T V  - a V + a b   V (2 R T V - a)
        ------------------ + ---------------
               3      2            3      2
         R T (V  - b V )     R T (V  - b V )
                                                 2                2
                                           V (3 V  - 2 b V) (R T V  - a V + a b)
                                         - -------------------------------------
                                                           3      2 2
                                                     R T (V  - b V )
 
% In the space below, calculate the derivative dZ/d(1/V) and call
% it dZd1V. HINT: dy/dx = -(dy/dz)/(dx/dz). 
% 
>: dZd1V : -dZdV/diff(1/V,V)
 
                     2
           2    R T V  - a V + a b   V (2 R T V - a)
        - V  (- ------------------ - ---------------
                       3      2            3      2
                 R T (V  - b V )     R T (V  - b V )
                                                2                2
                                          V (3 V  - 2 b V) (R T V  - a V + a b)
                                        + -------------------------------------)
                                                          3      2 2
                                                    R T (V  - b V )
 
% In the space below, calculate B as the limit of dZ/d(1/V) as V
% approaches infinity. Name it B.
% 
>: B : limit(dZd1V,V,inf)
 
                                     b R T - a
                                   - ---------
                                        R T
 
% In the space below, calculate the Boyle temperature as the temperature
% where B = 0. Call the Boyle temperature TB.
% 
>: TB: rhs(solve(B=0,T)[1])
 
                                        a
                                       ---
                                       b R
 
% Now let's see how well this equation works by looking at trends
% in Boyle temperature for the noble gases. Because we'll be calculating
% numbers and making graphs, we'll switch from Maxima back to Euler.
% 
% In the space below, define arrays a and b that contain the van
% der Waals constants for the noble gases He, Ne, Ar, Kr, and Xe,
% in that order. Take the numbers from Atkins Table 1.6 in the Data
% Section at the end of the book.
% 
>a = [0.0341,0.205,1.337,5.125,4.137]; 
>b = [2.38e-2,1.67e-2,3.20e-2,1.06e-2,5.16e-2];
% In the space below, assign R in units consistent with the units
% for a and b.
% 
>R = 0.082059;
% Next we will calculate an array of predicted Boyle temperatures
% for the noble gases, using the formula we derived in Maxima. Although
% the formula is so simple that it's easy to just retype it in Euler,
% let's see how you would move equations and expressions from Maxima
% into Euler.
% 
% The Euler function mxm() returns a named object in Maxima as a
% text string:
% 
>mxm("TB")
 a/(b*R)
% 
% To evaluate a text string in Euler, use the evaluate() function:
% 
>TB = evaluate(mxm("TB"))
      17.46027991133      149.5929648399      509.1610914098      5891.987058552      977.0338245001 
% 
% In the space below, plot the Boyle temperatures on the y axis
% with atomic number on the x axis. Insert the plot into this notebook.
% As always, title the plot and label the axes.
>atomicnumber = [2,10,18,36,54];
>plot2d(atomicnumber,TB);
>title("Boyle Temperature Trend for the Noble Gases");
>xlabel("Atomic Number");
>ylabel("Boyle Temperature (K)");
>insimg;

% 
% We'd expect the Boyle temperature to increase with increasing
% intermolecular attraction (a) and with decreasing molecular volume
% (b), so it's conceivable that the Boyle temperature might gradually
% rise and then fall going from He to Xe- but the huge spike in
% Boyle temperature at Kr is hard to miss. Either something is wrong
% with the data or calculations, or we have discovered something
% interesting about Kr.
% 
% In the space below, draw the plot again, but this time use van
% der Waals parameters from some other source (such as the CRC handbook).
>a = [0.03412,0.2107,1.345,2.318,4.194];
>b = [2.37e-2,1.709e-2,3.219e-2,3.78e-2,5.105e-2];
>TB = a/(b*R);
>plot2d(atomicnumber,TB);
>title("Boyle Temperature Trend for the Noble Gases");
>xlabel("Atomic Number");
>ylabel("Boyle Temperature (K)");
>insimg;

% 
% The moral of this last plot is: always double-check data points
% that don't fit the pattern implied by the rest of the data. Discordant
% points indicate either an error or an interesting discovery. 
% 
% The little cusp at Ar is interesting. To see if it's real, add
% a curve that shows the true Boyle temperatures as a thicker line
% on the plot. Insert the plot into this notebook again in the space
% below.
% 
>TB = [22.64,122.1,411.5,575,768];
>plot2d(atomicnumber,TB,add=1,thickness=2); insimg;

% 
% Explain the deviation between the true and van der Waals Boyle
% temperatures in terms of intermolecular forces. Does the van der
% Waals equation overestimate or underestimate the strength of intermolecular
% attractions?
% 
% ANSWER: The equation for the Boyle temperature we derived above
% shows that the larger a is, the higher the Boyle temperature will
% be. This suggests that the that the stronger intermolecular attractions
% are, the higher the Boyle temperature must be.
% 
% However the van der Waals predictions of the Boyle temperature
% are too high. That suggests that the van der Waals equation overestimates
% the strength of intermolecular attractions.
% 
—————————-snip here———————————————
h1

Notebook 1.3: Properties of Real Gases from the Ideal Gas Law

January 15, 2010

Real molecules attract and repel each other. Ideal gas molecules don’t. That’s what causes nearly all deviations from ideality.

The ideal gas law is derived by ignoring attractions and repulsions between molecules. It doesn’t work at all for solids and liquids, which have very strong intermolecular forces. It generally isn’t accurate for real gases, either.

We can still apply the ideal gas law in the limit of low pressure for real gases. At low pressure, the gap between molecules widens, weakening the forces acting between molecules until they are essentially negligible, so the equation applies.

The Euler notebook posted below asks students to connect data for real gases collected at several different pressures with the ideal gas law by extrapolating to zero pressure. It introduces linear regression with Euler.

To use the notebook:

* Download and install the Euler Math Toolbox.
* Cut and paste the code below into a plain text file.
* Save the file with an .en file extension.
* Double-click the file. It should open up in the Toolbox.

I used this in my Fall ‘09 PChem class (it’s the fourth notebook in the series). Italics mark the items students had to fill in themselves. They seemed to have little trouble with this one.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition.

—————————-snip here———————————————
% Notebook 1.3: Extrapolation to zero pressure
%
% In this notebook, we’ll calculate properties of real gases from the
% ideal gas law by extrapolating to zero pressure.
%
% Consider the following experimental data collected for gaseous
% dimethyl ether at 25 degrees Celsius. The pressures P are in pascals,
% and the densities d are in kg m^-3. R is in Pa m^3/mol K and T is
% K.
%
>P = [12.223,25.20,36.97,60.37,85.23,101.3]*1000;
>d = [0.225,0.456,0.664,1.062,1.468,1.734];
>R = 8.314;
>T = 25 + 273.15;
%
% Let’s find the molar mass of dimethyl ether from this data.
%
% The molar mass M of an ideal gas is mass per mole. The mass of the gas
% is its density times its volume, so
%
% M = dV/n
%
% For an ideal gas, PV=nRT or V/n = RT/P, so
%
% M = dRT/P
%
% The ideal gas law holds when P approaches zero, so for a real gas,
%
% M = lim dRT/P
% P->0
%
% This is called a limiting law. How can we apply the limiting law with the
% pressure and density measurements in P and d?
%
% First, calculate an array of dRT/P values.
%
>dRToverP = d*R*T/P;
%
% Let’s plot these values against P to see if we can linearly extrapolate
% them to zero.
% * When plotting experimental data, it is good practice to show
% the points on the graph (use the points=1 option).
%
>plot2d(P,dRToverP,points=1,style=”[]”);
>title(“Determining the Molar Mass of Dimethyl Ether”);
>xlabel(“Pressure, in Pa”);
>ylabel(“dRT/P, in kg/mol”);
>insimg;

%
% The data looks reasonably linear. It’s worth trying to fit a
% line to it. We can use Euler’s polyfit() function to do this. The
% call looks like
% a = polyfit(x,y,n)
% where:
%
% a is an array of coefficients in the fitted polynomial
% a[1] + a[2]*x + a[3]*x^2 + … + a[n+1]*x^n
% x and y are the arrays of x and y data
% n is the degree of the polynomial to fit (1 for a line)
%
>fit = polyfit(P,dRToverP,1)
0.04587503092209 -3.583072735655e-008
%
% The first number is the intercept; the second is the slope.
%
% To calculate an array of points along the regression line, use the
% polyval() function. The call looks like
% y = polyval(a,x)
% where x, y, and a have been defined above.
%
>x = [0,max(P)];
>y = polyval(fit,x);
%
% Now insert the plot, with the raw data points marked as boxes,
% and the regression line displayed as a solid line.
%
>plot2d(x,y);
>plot2d(P,dRToverP,points=1,style=”[]”,add=1);
>title(“Determining the Molar Mass of Dimethyl Ether”);
>xlabel(“Pressure, in Pa”);
>ylabel(“dRT/P, in kg/mol”);
%
% The y-intercept is M, the value of dRT/P at P=0.
%
>M = fit[1]
0.04587503092209
%
% We can label the y-intercept using the label function. The call
% looks like label(text,x,y,offset). You can play around with the
% offset parameter until the label is where you want it.
%
>label(“M”,0,M,0);
>insimg;

% —————————————————————–
%
% Apply these techniques to solve problem 1.3 in the text, which
% finds absolute zero on the Celsius temperature
% scale from experimental data.
%
% Enter the P and alpha data in the space below.
% (HINT: look carefully at the units.)
%
>P = [749.7, 599.6, 333.1, 98.6];
>alpha = [3.6717, 3.6697, 3.6665,3.6643]*1e-3;

%
% It will be easiest to find alpha for an ideal gas first.
% Do this with the polyfit() function in the space below.
% Call the alpha value for an ideal gas ‘alphaideal’.
>fit = polyfit(P,alpha,1)
0.003662977448051 1.139259281072e-008
>alphaideal = fit[1];

%
% In the space below, insert a plot that extrapolates alpha to zero
% pressure. Mark the data on the line as circles (use style=”o”) and
% add a regression line to the plot. The y-intercept should be shown
% and labeled as “alpha ideal”.
%
>x = [0,max(P)];
>y = polyval(fit,x);
>plot2d(x,y);
>plot2d(P,alpha,points=1,style=”o”,add=1);
>title(“Determining alpha for an ideal gas”);
>xlabel(“Pressure, in torr”);
>ylabel(“alpha, in deg. C^-1″);
>label(” alpha ideal”,20,fit[1],-25);
>insimg;


%
% Now solve Charles’ Law in the form V=Vo(1+alpha*theta) for
% absolute zero. Call the solution ‘absolutezero’.
% HINT: You won’t need an actual value for Vo, if you do your
% algebra correctly.
% HINT: At V=0, theta is absolute zero on the Celsius scale.
% HINT: Charles’ Law applies exactly only to ideal gases.
%
>absolutezero = -1/alphaideal
-273.0019537882

%
% ———————————————————————–
%
% Apply these techniques again to solve problem 1.7(b) in the text,
% which estimates the value of the ideal gas law constant R from
% experimental data.
%
% Enter the P, V, and density(d) data in the space below.
%
>P = [0.750000, 0.500000, 0.250000];
>V = [29.9649, 44.8090, 89.6384];
>d = [1.07144, 0.714110, 0.356975];

%
% Calculate an array of R values calculated at each of the three points.
% Call it RP (for “R at pressure P”).
%
>M = 2*15.9994;
>T = 273.15;
>RP = P*M/(d*T);

%
% In the space below, insert a plot of the RP values as a function
% of pressure. Show the data points as boxes (use style=”[]”). Draw
% the regression line on the plot. AS ALWAYS title the plot and
% label the axes. The axis labels should include the units.
%
>fit = polyfit(P,RP,1)
0.08206187029697 -7.886858741291e-005
>x = [0,max(P)];
>y = polyval(fit,x);
>plot2d(x,y);
>plot2d(P,RP,points=1,style=”[]”,add=1);
>title(“Determination of R”);
>xlabel(“Pressure, in atm”);
>ylabel(“P*M/(density*T), in dm^3 atm/mol K”);
>label(“R”,0,fit[1],0);
>insimg;


%
% In the space below, print your calculated R value.
>R = fit[1]
0.08206187029697

h1

Notebook 1.2 – Solving Equations of State (Molar Volumes of Model Gases)

January 14, 2010

The Euler notebook posted below asks students to solve various gaseous equations of state for molar volume. This is a common task that’s quite easily done with a computer algebra system, but difficult and tedious to do with pencil and paper.

To use the notebook:

  • Download and install the Euler Math Toolbox.
  • Cut and paste the code below into a plain text file.
  • Save the file with an .en file extension.
  • Double-click the file. It should open up in the Toolbox.

I used this in my Fall ‘09 PChem class (the third in the series). Italics mark the items students had to fill in themselves. They had less trouble with this one, although command line syntax was still a source of confusion and frustration. They wanted to type something like root(R*T/V-P,V) or root(RT/V-P,V) instead of root(“R*T/V-P”,V).

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition.


----------------------------snip here---------------------------------------------
Notebook 1.2 - Molar Volume of Model Gases (Solving Equations of State)

%
% In this task, we'll use Euler to numerically solve equations of state.
% Study Example 1.4 on page 18 of Atkins, which estimates the molar
% volume of CO2 at 500 K and 100 atm by treating CO2 as a van der Waals
% gas.
%
% To solve this problem in Euler, set up constants in any convenient
% units (in this case, dm^3 for volume, atm for pressure, and K for
% temperature:
>R = 0.082059;
>a = 3.592;
>b = 4.267e-2;
>T = 500;
>P = 100;
%
% Now cast the equation into a form that looks like f(x) = 0.
% Atkins rearranges the van der Waals equation into polynomial
% in V:
% V^3-(b+R*T/P)*V^2 + (a/P)*V - a*b/P = 0.
% but that really isn't necessary. *Any* rearrangement of the equation
% we want to solve that puts zero on one side of the equation will work.
%
% To solve the equation in Euler, we'll need an initial guess for V.
% The ideal gas volume will do.
>V = R*T/P
0.410295
%
% Now use Euler's 'root' function, which looks like
% root("expression",x)
% where expression (enclosed in double quotes) is the nonzero side of the
% equation we want to solve, and x is the variable in that expression we
% want to solve for.
%
>root("V^3-(b+R*T/P)*V^2 + (a/P)*V - a*b/P",V)
0.3663332616871
%
% Note that root has reset the value of V:
%
>V
0.3663332616871
%
% Atkins started with the equation
% P = R*T/(V-b) - a/V^2
% We could have subtracted P from both sides of this equation to more
% simply cast it into the f(x) = 0 form:
%
>root("R*T/(V-b)-a/V^2-P",V)
0.3663332616871
%
% Atkins claims that a perfect gas under these conditions has a molar
% volume of 0.410 dm^3/mol. We can verify this by typing
%
>root("R*T/V - P", V)
0.4102949999999
% ---------------------------------------------------------------------
% TASK 3. Use this second approach to calculate the molar volume of
% argon in dm^3 per mole at 100 degrees Celsius and 100 atm, on the
% assumption that it's
% A) an ideal gas;
% B) a hard sphere gas (assuming that each molecule is a
% sphere with radius of 0.15 nm);
% C) a van der Waals gas;
% D) a Dieterici gas;
% E) a virial gas (ignoring all virial coefficients beyond B).
%
>T = 100+273.15;
>P = 100;
%
% 3A) The molar volume as an ideal gas:
%
>V = R*T/P
0.3062031585

%
% 3B) The molar volume as a hard sphere gas, with radius 0.15 nm:
% HINT: In Euler, the constant pi is written as %pi
% HINT: What units should you put b into
% if V is in dm^3/mol?
%
>r = 0.15*1e-9/1e-1;
>NA = 6.022e23;
>b = 4*(4/3)*%pi*r^3*NA;
>root("R*T/(V-b)-P",V)
0.3402567662277

%
% 3C) The molar volume as a van der Waals gas:
%
>a = 1.337;
>b = 3.20e-2;
>root("(P+a/V^2)*(V-b)-R*T",V)
0.2981759820445

%
% 3D) The molar volume as a Dieterici gas:
% HINT: See Table 1.7 in Atkins.
% HINT: a and b for the Dieterici gas are NOT the same as
% a and b for a van der Waals gas!
% HINT: All of the expressions in the Critical Constants columns
% in Table 1.7 are equal.
% HINT: In Euler, the constant e is written as %e
% HINT: In Euler, e^x can be written as exp(x)
%
>B = 3*b/2;
>A = 4*%e^2*B^2*a/(27*b^2);
>root("R*T*exp(-A/(R*T*V))/(V-B)-P",V)
0.245635997406

%
% 3E) The molar volume as a virial gas (with all virial coefficients
% beyond B ignored)
% HINT: Watch your units; what should the units of B be, if V
% should be in dm^3/mol?
% HINT: B varies with temperature.
%
>B = -4.2/1000;
>root("(R*T/V)*(1+B*P)-P",V)
0.17759783193

----------------------------snip here---------------------------------------------

h1

Notebook 1.1: Atmospheric Modeling

January 14, 2010

Famous right-wing radio personality Rush Limbaugh claims that chlorofluorocarbons can’t possibly deplete ozone in the stratosphere, because they’re heavier than air. They’ll settle close to the ground, says Rush, where they can’t do any harm.

The Euler notebook posted below asks students to apply the barometric equation to see whether or not gases really stratify in the atmosphere based on their molar masses. To use it, download the Euler Math Toolbox. After installing the toolbox, just cut and paste the following code into a plain text file, save it with an .en file extension, and double-click the file. It should open up in the Toolbox.

I used this in my Fall ‘09 PChem class (the second in the series). Italics mark the items students had to fill in themselves. Some of my students struggled with this assignment because they couldn’t fathom doing calculations with arrays. I really should have had a notebook before this one that focused exclusively on arithmetic with arrays.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition.


---------------------------------snip here-----------------------------
Notebook 1.1 - Atmosphere Modeling
%
% As before, hit return on each line to execute the calculation on that
% line. Use ALT-Insert to insert a new line, and ALT-Delete to remove a
% line.
%
% -----------------------------------------------------------------------
%
% In this notebook, we'll compare two models for describing the
% change in air pressure with altitude. Both models are based on the
% barometric formula
% P = Po*exp(-h*M*g/(R*T))
% where P is the pressure, Po is the pressure at altitude 0, h is the
% altitude, M is the molar mass, g is the acceleration due to gravity,
% R is the ideal gas constant, and T is the temperature in kelvins.
% * MODEL A assumes that air behaves like an ideal gas with a
% single molar mass, Mo.
% * MODEL B assumes that each gas in air behaves as an ideal
% gas, and each gas develops its own separate altitude
% profile. In this model, we should see that air is richer in
% heavy gases at lower altitudes.
%
% TASK 2A. Let's begin with model A by using the barometric formula to plot
% air pressure (in atm, on the x-axis) against altitude (on the y-axis).
%
% First, we need to know how temperature varies with altitude.
% We can put the temperature vs. altitude data into lists called
% 'arrays'. Arrays can be entered into Euler notebooks by using
% square brackets. For example, here are the data arrays for altitude
% (h, in meters) and temperature (t, in degrees Celsius).

>h=[-1000,0,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000,15000,20000,25000,30000,40000,50000,60000,70000,80000];
>t=[21.50,15.00,8.50,2.00,-4.49,-10.98,-17.47,-23.96,-30.45,-36.94,-43.42,-49.90,-56.50,-56.50,-51.60,-46.64,-22.80,-2.5,-26.13,-53.57,-74.51];

% Each altitude in h corresponds to a temperature in t.
%
% We can plot arrays and show the data as separate
% points using plot2d's 'points=1' option. We can choose the style of the
% points using the 'style' option. Style can be set to "o", ".", "x",
% or "[]" with points=1.
%
>plot2d(t,h);
>plot2d(t,h,points=1,add=1,style="o");
%
% Always title your plots, and label the axes. Include units in the
% axis labels.
%
>title("Variation of Temperature with Altitude");
>xlabel("temperature (degrees Celsius)");
>ylabel("altitude (m)");
>insimg;

%
% Next, we need the average molar mass of air at sea level, Mo (in
% kg/mol). The molecule fractions in air are 78.09% N2, 20.95%
% O2, 0.93% Ar, and 0.03% CO2 at sea level, with negligible
% amounts of other gases in dry air. We can use arrays again, with
% each separate gas as an 'element' in the array (N2 is element 1,
% O2 is element 2, Ar is element 3, and CO2 is element 4):
%
>percent = [78.09, 20.95, 0.93, 0.03];
>abundance = percent/100;
>mass = [2*14.0067, 2*15.9994, 39.948, 12.0107+2*15.9994];
%
% If we're going to do our calculations in SI units, we need the
% masses in kg/mol, not g/mol. Note that we can do the conversion
% on all elements of the array at once.
%
>mass = mass/1000;
%
% We can access the i-th element of array x by writing x[i].
% For example, the mass of N2 is mass[1]. The abundance of Ar is
% abundance[3].
%
>Mo = abundance[1]*mass[1] + abundance[2]*mass[2] + abundance[3]*mass[3] + abundance[4]*mass[4]
0.02896413191
%
% It's much easier to use the sum function, which adds up all of the
% elements of an array.
%
>Mo = sum(abundance * mass)
0.02896413191
%
% We also need the acceleration due to gravity (g, in m/s^2). We can't
% get accurate results by just plugging g=9.8 into the barometric
% equation because g decreases with altitude. Each number in the
% following array gives g at the corresponding altitude in the h array.
%
>g=[9.810,9.807,9.804,9.801,9.797,9.794,9.791,9.788,9.785,9.782,9.779,9.776,9.761,9.745,9.730,9.715,9.684,9.654,9.624,9.594,9.564];
%
% Also assume that the pressure Po at sea level is 1 atm.
%
>Po=1;
%
% Now calculate an array P that contains the pressures predicted
% by the barometric formula at each altitude h.
%
>R=8.314;
%
% Remember that names in Euler are case sensitive. We can use T
% for the kelvin temperature, and t for the Celsius temperature.
%
>T = t+273.15;
%
% Now use the barometric equation to calculate the total pressure,
% assuming that the average molar mass of air is the same at every
% altitude.
%
>P=Po*exp(-h*Mo*g/(R*T));
>plot2d(P,h);
>plot2d(P,h,points=1,add=1,style="o");
>title("MODEL A: Barometric Pressure Profile with Fixed Air Molar Mass");
>xlabel("Barometric Pressure (atm)");
>ylabel("altitude (m)");
>insimg;


%
% ----------------------------------------------------------------------
%
% TASK 2B. In MODEL A, we assumed that the average molar mass of air
% is the same at every altitude. However, each separate gas might
% have its own altitude profile, with heavier gases being more
% abundant at lower altitudes than light gases.
%
% On a second plot, use the barometric equation to show
% the partial pressures for each of the four gases as a function of
% altitude. Put the four curves on the same plot, in different colors.
%
% Call the arrays holding the partial pressures of each gas
% PN2, PO2, PAr, and PCO2.
%
>PN2 = (Po*abundance[1])*exp(-h*mass[1]*g/(R*T));
>plot2d(PN2,h,color=1);
>PO2 = (Po*abundance[2])*exp(-h*mass[2]*g/(R*T));
>plot2d(PO2,h,add=1,color=2);
>PAr = (Po*abundance[3])*exp(-h*mass[3]*g/(R*T));
>plot2d(PAr,h,add=1,color=3);
>PCO2 = (Po*abundance[4])*exp(-h*mass[4]*g/(R*T));
>plot2d(PCO2,h,add=1,color=4);

%
% Then use Dalton's Law to get a fifth curve,
% the total pressure profile. Call it Ptot.
%
>Ptot = PN2 + PO2 + PAr + PCO2;
%
% Plot the curve for the total pressure with a heavier line using
% 'thickness'. For example, to plot y as a function of x with a
% heavier curve, you could type plot2d(x,y,thickness=2).
%
>plot2d(Ptot,h,add=1,thickness=2,color=5);
>title("MODEL B: Barometric Pressure Profile with Variable Molar Mass");
>xlabel("Pressure (atm)");
>ylabel("Altitude (m)");
>insimg;


%
% ----------------------------------------------------------------------
%
% TASK 2C. Calculate the variation in average molar mass with altitude
% for MODEL B. You can do this by solving the barometric equation
% for average molar mass M, using pressure Ptot (from Task 2B).
% HINT: You may run into some practical difficulties in your calculation
% of M at h = 0, because
% * division by zero is not allowed
% * ln (Ptot/Po) = 0 when h=0, so you may get molar mass = 0 at
% altitude zero.
% The average molar mass at h=0 should be Mo. You can set the second
% element of M to Mo by typing M[2] = Mo.
%
>h[2]=1; M = -R*T/(h*g)*ln(Ptot/Po);
>M[2]=Mo;
>M
Column 1 to 3:
0.02897184063066 0.02896413191 0.02895617819929
Column 4 to 6:
0.02894796718375 0.02893948906187 0.02893072545567
Column 7 to 9:
0.02892166297037 0.02891228549505 0.02890257572915
Column 10 to 12:
0.02889251507121 0.0288820868212 0.02887126708741
Column 13 to 15:
0.02882661789134 0.02878778498486 0.02875566381133
Column 16 to 18:
0.02872691182118 0.02868986042015 0.02865867484957
Column 19 to 21:
0.02859332001119 0.02852196225203 0.02845761067037

% Insert a plot the average molar mass M for model B as a function of
% altitude.
%
>plot2d(M,h);
>title("MODEL B: Variation of Molar Mass with Altitude");
>xlabel("Average Molar Mass (kg/mol)");
>ylabel("Altitude (m)");
>insimg;


%
% ----------------------------------------------------------------------
%
% TASK 2D. Assess the accuracy of the two methods for calculating
% the absolute error in total pressure for the two methods against
% altitude, using measured pressures in the array Ptrue for the
% altitudes in h. Ptrue is in atmospheres.
%
>Ptrue = [1.1229908, 1.000000, 0.885790, 0.780201, 0.683079, 0.594156, 0.513204, 0.439950, 0.374090, 0.315296, 0.263226, 0.217485, 0.0949375, 0.0435330, 0.0218130, 0.0113010, 0.00455875, 0.00200189, 0.000290389, 2.35653e-5, 1.48565e-6];
>plot2d(Ptrue - Ptot, h, color=2, thickness=1);
>plot2d(Ptrue - P, h, add=1, color=3, thickness=5);
>title("Absolute Error in Model B (thinner curve) and Model A (thicker curve)");
>xlabel("Pressure (atm)");
>ylabel("Altitude (m)");
>insimg;


%
% Which of the two models is more accurate? Explain.
%
>"Model A is more accurate. This suggests that the gases in air don't develop separate altitude profiles.";
>"A reasonable hypothesis is that winds remix the gases that gravity would separate.";
>