{VERSION 4 0 "IBM INTEL NT" "4.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 2 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 276 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 277 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 278 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 279 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 282 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 283 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 286 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 287 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 288 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 289 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 } {CSTYLE "" -1 290 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 291 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 292 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 293 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 294 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 295 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 296 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 297 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 298 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 299 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 300 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 301 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 302 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 303 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 304 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 305 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 306 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 307 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 308 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 309 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 310 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 311 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 312 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 313 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 314 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 16 0 0 0 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 256 1 {CSTYLE " " -1 -1 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT 256 20 "Ideal Gas Simulation" } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 1164 "In thi s worksheet a simple 2-dimensional ideal gas simulation is performed. \+ Particles are moving in a 2D container and collide with the walls. Sin ce we are simulating a dilute, ideal gas, we can ignore the collisions between the gas molecules. This is also called a collisionless plasma . The effect of collisions with the walls is to couple the x and y mot ions, and to define the pressure of the gas in the container. The temp erature of the gas is related to the total kinetic energy. Our gas is \+ simplified in one important aspect: the molecules are treated as point particles, and thus are able to store kinetic energy only through the ir translational motion. In reality, diatomic (or polytaomic) molecule s can convert heat into electronic excitations (eV energy scale), vibr ations (0.1 eV scale), and rotations (0.01 eV scale). There are quantu m manifestations of these phenomena, i.e., abrupt changes in the speci fic heat as a function of temperature when a storage mechanism becomes unavailable as the temperature T is lowered. It is useful to know the energy associated with temperature, i.e., kT , which at room temperat ure corresponds to 1/40 eV (0.025 eV)." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{EXCHG {PARA 0 "" 0 "" {TEXT 267 34 "The Maxwell-Boltzmann distri bution" }}{PARA 0 "" 0 "" {TEXT -1 963 "We begin with the Boltzmann di stribution and learn how to generate a random sample with a non-unifor m distribution. Our simulation cannot be completely free of assumption s, since it does not describe the relaxation process that takes place \+ over a long time scale: much longer than the time it takes for a parti cle to traverse the box, i.e., the time scale of our simulation. The r elaxation to a steady state from any intial state takes place as a res ult of collisions between the molecules. The result is the Maxwell-Bol tzmann distribution which describes how the particles of the ensemble \+ are distributed over kinetic energy (magnitude of velocity) for given \+ temperature. First-year physics texts provide the main principle behin d a derivation. Collisions between the gas molecules imply that the in dividual molecules can have any amount of kinetic energy, since they e xchange it all the time amongst each other. The likelihood to observe \+ a given kinetic energy " }{XPPEDIT 18 0 "mv^2/2" "6#*&%#mvG\"\"#F%!\" \"" }{TEXT -1 239 " depends on the number of possibilities to reach th at state. To observe a very high, or a very low kinetic energy at give n temperature is less likely than some value in between. The average k inetic energy is given by the following relation" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "avg:=m/2*u^2=d/2*k*T;" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 136 "Here d is the number of dimensions (d=3 in nature, d=2 for our simulation), k is Boltzmann's constant, T is the temperature \+ of the gas, " }{TEXT 257 1 "m" }{TEXT -1 44 " the mass of the particle s (molecules), and " }{TEXT 258 2 "u " }{TEXT -1 57 "is the average(!) velocity. The distribution has to have " }{TEXT 259 1 "u" }{TEXT -1 124 " for its average speed, but note that the average velocity itself (as a vector) vanishes. For d=3 the distribution reads as:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "fB:=(vx,vy,vz,T)->(m/(2*Pi*k *T))^(3/2)*exp(-m*(vx^2+vy^2+vz^2)/(2*k*T));" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 121 "Note how the symmetry of the distribution will render \+ a zero result for an expectation value of either of the components " } {XPPEDIT 18 0 "v[x ], v[y] , v[z] " "6%&%\"vG6#%\"xG&F$6#%\"yG&F$6#%\" zG" }{TEXT -1 78 "; it is more convenient to define a function that de pends on the magnitude of " }{XPPEDIT 18 0 "[v[x], v[y], v[z]]" "6#7%& %\"vG6#%\"xG&F%6#%\"yG&F%6#%\"zG" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "fB:=(v,T)->(m/(2*Pi*k*T))^(3/2)*exp(-m*v^2/(2 *k*T));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "int(fB(v,T)*v^2, v=0..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "assume(m >0): assume(k>0): assume(T>0): interface(showassumed=0):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 126 "We perform now the calculation of the av erage in 1D, 2D, 3D: the volume integration in velocity space contribu tes a factor 1, " }{XPPEDIT 18 0 "v" "6#%\"vG" }{TEXT -1 6 ", and " } {XPPEDIT 18 0 "v^2" "6#*$%\"vG\"\"#" }{TEXT -1 14 " respectively." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "int(fB(v,T)*v^2,v=0..infinit y)/int(fB(v,T),v=0..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "int(fB(v,T)*v^3,v=0..infinity)/int(fB(v,T)*v,v=0..infinity);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "int(fB(v,T)*v^4,v=0..infi nity)/int(fB(v,T)*v^2,v=0..infinity);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 768 "Now we need to make a choice of units. It is possible (a lbeit cumbersome) to always carry the units in Maple (multiplied onto \+ the numerical values for the constants). One can define a natural syst em of units that renders the equations dimensionless. Alternatively on e can fix a popular system of units (cgs or MKSA=SI), ignore the units and supply them at the end of the calculation. This is less safe than either of the previous methods, and requires the diligence to supply \+ all variables and constants in the same unit system. We choose the MKS A system of units, for k we omit the J/K (J = Nm = Ws). The Boltzmann \+ constant can be determined from the observation of Brownian motion (Mi llikan experiment to determine the unit of charge), and is given appro ximately as:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "k:=1.38*10^ (-23);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 80 "The mass of the gas par ticles in kg is assumed as a multiple of the proton mass:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "m:=2*1.67*10^(-27);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 118 "We are ready to graph M-B distributions \+ for a few temperatures. We include the volume element for a 2D gas (fa ctor of " }{TEXT 260 1 "v" }{TEXT -1 61 ") to observe the probability \+ to find particles in the range [" }{TEXT 265 1 "v" }{TEXT -1 2 ", " } {TEXT 264 1 "v" }{TEXT -1 2 "+d" }{TEXT 263 1 "v" }{TEXT -1 52 "]. Bef ore we graph we need to figure out ranges for " }{TEXT 262 1 "x" } {TEXT -1 5 " and " }{TEXT 261 1 "y" }{TEXT -1 1 ":" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "evalf(fB(v,100));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "plot([10^8*v*fB(v,100),10^8*v*fB(v,300),10^8*v*f B(v,1000)],v=0..0.5*10^4,color=[red,blue,green]);" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 273 "Note the increase in the most likely speed (in m/ sec) as the temperature changes from 100 to 1000 degrees Kelvin. One c an see that the distributions are not properly normalized yet, the are a changes as the temperature is increased. For this reason the expecta tion value of " }{XPPEDIT 18 0 "v^2" "6#*$%\"vG\"\"#" }{TEXT -1 254 " \+ was normalized by dividing through the area under the distribution cur ve. Note that there is a difference between the position of the maximu m, and the average value (We leave it as a small exercise to determine both values for a few given temperatures)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 266 41 "Simulations in Maple u sing random numbers" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 190 "We proceed with random number generation. Maple has a \+ built-in generator of random integers which can be easily modified to \+ produce equidistributed random numbers on the unit interval [0,1]." }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "# info can be obtained from ?rand" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "randomize(Seed); \+ # this can be changed to generate another sequence" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "randx:=()->evalf(rand()/999999999999);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "randx();" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 142 "A way to check the quality of uniform random numb ers is to generate a sequence of points, to display them and observe w hether patterns appear." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 " for i from 1 to 500 do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "xv[i]:=r andx(): yv[i]:=randx(): od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "i:='i': plot([[xv[i],yv[i]]$i=1..500],0..1,0..1,style=point, axes= boxed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 631 "Ideally we should fin d a uniformly covered square. To assess the quality of the generator o ne has to observe how the square is covered as the number of randoms i s increased. The development and testing of pseudorandom number genera tors is a science (or an art) in its own right. Note that we are using seeded reproducible (deterministic) sequences (thus, pseudorandom num bers) with certain statistical properties. More quantitative ways to t est these generators involve the observation of the mean (does it appr oach 0.5 with a 1/sqrt(N) deviation behaviour), and more generally the quality of the histogram over the unity interval." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "histo:=proc (xv,N,xi,xf,dx) local bin,n,i; global nm;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "nm:=trunc((xf-xi)/dx); for n from 1 to nm do: bin[n]: =0; od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to N do:" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "n:=trunc((xv[i]-xi)/dx)+1;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "bin[n]:=bin[n]+1; od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "bin;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "res:=histo(xv,500,0,1,0 .1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "seq(res[i],i=1..nm) ;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 89 "It is not appropriate to gra ph these numbers directly. A bar graph is carried below. The " }{TEXT 19 7 "op([,])" }{TEXT -1 142 " construct is used to place two neighbor ing points beside each other in the sequence construct, which is then \+ turned into a list for graphing." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "i:='i': plot([seq(op([[i-1,res[i]],[i,res[i]]]),i=1.. nm)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 268 11 "Exercise 1:" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 110 "Observe the deviations from t he expected average value (50 for a sample of 500) for each bin as a f unction of " }{TEXT 275 1 "N" }{TEXT -1 28 ". How does it decrease wit h " }{TEXT 274 1 "N" }{TEXT -1 1 "?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 175 "We can turn the histogram gener ation into a single procedure (to which we supply the maximum count va lue as well), and try it out on the other chosen sample of random numb ers:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "histo:=proc(xv,N,xi ,xf,dx,ymax) local bin,n,i,nm;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "n m:=trunc((xf-xi)/dx); for n from 1 to nm do: bin[n]:=0; od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "for i from 1 to N do:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 26 "n:=trunc((xv[i]-xi)/dx)+1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "bin[n]:=bin[n]+1; od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "i:='i': plot([seq(op([[i-1,bin[i]],[i,bin[i]]]),i=1..nm)],0..n m,0..ymax);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "histo(yv,500,0.,1.,0.1,80);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 465 "To dial up an ensemble that is no t equidistributed but follows a given distribution we use a rejection \+ method. First we truncate the infinite range of velocities to a finite one with a cut-off appropriate for the given temperature (see the gra ph). For finite ensembles this constitutes a small error. Alternativel y we could remap the independent variable v , using e.g. the tangent f unction, such that the infinite range would be mapped to a finite rang e (0..pi/2)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 155 "For T<1000 K we should be able to truncate the velocities at 1 0,000 m/s. We dial up a velocity in this range and call it v_i. Next w e dial a random number " }{XPPEDIT 18 0 "xi " "6#%#xiG" }{TEXT -1 60 " in the range 0..1, and calculate the distribution function " } {XPPEDIT 18 0 "fB(v[i])" "6#-%#fBG6#&%\"vG6#%\"iG" }{TEXT -1 84 " norm alized such that at the maximum of fB a value of 1 is obtained. Now we compare " }{XPPEDIT 18 0 "xi " "6#%#xiG" }{TEXT -1 5 " and " } {XPPEDIT 18 0 "fB(v[i])" "6#-%#fBG6#&%\"vG6#%\"iG" }{TEXT -1 54 " : we accept v_i as a member of the ensemble only if " }{XPPEDIT 18 0 "fB( v[i]) " 0 "" {MPLTEXT 1 0 9 "T0 :=1000;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "f0:=evalf(v*fB(v ,T0));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "v0:=solve(diff(f0 ,v)=0,v);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "f0max:=evalf(s ubs(v=abs(v0[1]),f0));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "f 0n:=f0/f0max;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "vmax:=1000 0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "Np:=100;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 109 "The rejection method requires an infinit e loop. For practical reasons we truncate that loop at 1000 attempts. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "#begin loop:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to Np do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "vi:=vmax*randx();" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "f0i:=subs(v=vi,f0n);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "xi:=r andx();" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "for j from 1 to 1000 whi le xi>f0i do: #1000 attempts" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "vi: =vmax*randx();" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "f0i:=subs(v=vi,f0 n);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "xi:=randx(); od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "vB[i]:=vi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "histo(vB,Np,0,v max,vmax/20,Np/4);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 276 12 "Exercise 2 : " }}{PARA 0 "" 0 "" {TEXT -1 80 "Modify the procedure histo to displ ay the velocity instead of bin number on the " }{TEXT 277 1 "x" } {TEXT -1 151 " axis, and compare the graph with the distribution. The \+ maximum in the fourth bin appears to be in the correct place, i.e., ne ar 2000 m/s for T=1000 K." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 90 "In our next step we should dial up veloci ty vectors. A randomly dialed polar angle in the " }{TEXT 272 1 "x" } {TEXT -1 1 "-" }{TEXT 271 1 "y" }{TEXT -1 33 " plane should allow us t o define " }{TEXT 270 1 "x" }{TEXT -1 5 " and " }{TEXT 269 1 "y" } {TEXT -1 174 " components of the initial velocities. Theta varies betw een -pi and pi, and the uniform distribution is obtained by stretching the random numbers obtained on a unit interval." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to Np do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "theta:=evalf((randx()-0.5)*2*Pi);" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 53 "vx[i]:=vB[i]*cos(theta); vy[i]:=vB[i]*sin(theta); o d:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "We determine some averages \+ for our discrete ensemble." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "vxavg:=0: vyavg:=0: vavg:=0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to Np do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "vxavg:=v xavg+vx[i]: vyavg:=vyavg+vy[i]: vavg:=vavg+sqrt(vx[i]^2+vy[i]^2): od: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "vxavg:=vxavg/Np: vyavg:=vyavg/N p: vavg:=vavg/Np:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "print( `vx,vy,v: `,vxavg,vyavg,vavg);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 216 "Note that the components of the average velocity vector should go to zero as the number of test particles is increased, but that the ma gnitude of the average velocity vector should be given by the ideal ga s estimate:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "evalf(sqrt(v xavg^2+vyavg^2),3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "eval f(sqrt(2*k*T0/m),5); evalf(vavg,5);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 106 "A proper histogram of the velocity components can be obtained \+ by using Maple's statistics package. On the " }{TEXT 278 1 "x" }{TEXT -1 52 "-axis we will find velocities (in m/sec), while the " }{TEXT 279 1 "y" }{TEXT -1 108 "-axis does not show counts for the bins, but \+ the ratio of counts to total particle number (count frequency)." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "with(stats):\nwith(stats[sta tplots]):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "histogram(conv ert(vx,list),color=cyan,area=1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 273 33 "Simulation of the particle motion" }}{PARA 0 "" 0 "" {TEXT -1 50 " The container is defined through its dimensions: -" }{TEXT 280 1 "L" } {TEXT -1 3 " < " }{TEXT 281 1 "x" }{TEXT -1 3 " < " }{TEXT 282 1 "L" } {TEXT -1 8 " , and -" }{TEXT 283 1 "L" }{TEXT -1 3 " < " }{TEXT 284 1 "y" }{TEXT -1 3 " < " }{TEXT 285 1 "L" }{TEXT -1 162 ". To propagate t he ensemble of molecules we can use the solution to Newton's equation \+ for rectilinear motion with possible reflections from one of the four \+ walls." }}{PARA 0 "" 0 "" {TEXT -1 78 "W choose a size for the contain er in meters and pick random initial positions." }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 7 "L:=0.1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to Np do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "x [i]:=evalf((randx()-0.5)*2*L); y[i]:=evalf((randx()-0.5)*2*L); od:" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "i:='i': plot([[x[i],y[i]]$i =1..Np],-L..L,-L..L,style=point,axes=boxed);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 201 "To propagate the particles we need to determine the ti me at which they collide with one of the walls. A scan has to be perfo rmed and the particles are updated at each time that one of them hits \+ a wall." }}{PARA 0 "" 0 "" {TEXT -1 116 "We can limit the run time of \+ the simulation by fixing the total number of collisions with the walls (our choice for " }{TEXT 19 5 "ntmax" }{TEXT -1 49 " corresponds to 2 0 wall encounters per particle)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 72 "An average collision interval or time sca le for our problem is given by:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "tscale:=2*L/vavg;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 38 "Insi de the loop we determine the time " }{TEXT 19 5 "tcoll" }{TEXT -1 102 " at which the next particle will encounter one of the walls. This che ck is done independently for the " }{TEXT 296 1 "x" }{TEXT -1 6 "- and " }{TEXT 295 1 "y" }{TEXT -1 131 "-components of each particle's moti on. A compact expression can be found that is valid irrespective of th e particle direction: for " }{TEXT 297 1 "v" }{TEXT -1 15 ">0 one rela tes " }{TEXT 298 3 "v t" }{TEXT -1 8 "_c with " }{TEXT 302 1 "L" } {TEXT -1 1 "-" }{TEXT 301 1 "x" }{TEXT -1 6 ", for " }{TEXT 300 1 "v" }{TEXT -1 15 "<0 one relates " }{TEXT 299 3 "v t" }{TEXT -1 9 "_c with -" }{TEXT 303 1 "L" }{TEXT -1 1 "-" }{TEXT 304 1 "x" }{TEXT -1 50 ". \+ This can be combined into one expression using " }{TEXT 306 1 "L" } {TEXT -1 6 " sign(" }{TEXT 305 1 "v" }{TEXT -1 2 ")." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 199 "All other particle pos itions are updated, and the velocity component perpendicular to the wa ll is reversed for the colliding particle. In addition we sum the mome ntum transfer to the container (mass " }{TEXT 287 1 "m" }{TEXT -1 118 " times twice the velocity component perpendicular to the wall - this \+ comes from elastic collisions of objects of mass " }{TEXT 286 1 "m" } {TEXT -1 229 " with an infinitely heavy object). The force impulse due to the collision can be converted to an incremental pressure increase on the container (by dividing the force by the container area, and av eraging over time). The variable " }{TEXT 19 4 "Delv" }{TEXT -1 59 " a ccumulates the momentum transfers (apart from the factor " }{TEXT 288 1 "m" }{TEXT -1 46 "). The pressure information is accumulated in " } {TEXT 19 1 "P" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "tcoll:=tscale: P[0]:=0:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "ti m:=0: Delv:=0: Area:=8*L: Vol:=(2*L)^2: ntmax:=Np*20:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "for ntim from 1 to ntmax do: #step through colli sions" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "for i from 1 to Np do: #de termine shortest collision time" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 " tx[i]:=(L*sign(vx[i])-x[i])/vx[i]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "if tx[i] " 0 "" {MPLTEXT 1 0 85 "ty[i]:=(L*sign(vy[i])-y[i])/vy[i]:if ty[i] " 0 "" {MPLTEXT 1 0 38 "tim:=tim+tcoll: #now update positions:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to Np do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 80 "x[i]:=evalf(x[i]+vx[i]*tcoll): y[i]:=evalf(y[i]+vy[i]*tcoll): \+ od: tcoll:=tscale:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "#now reverse \+ the momentum (velocity) of colliding TP:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 142 "if tx[icoll] " 0 "" {MPLTEXT 1 0 27 "P[ntim]:=m* Delv/(Area*tim):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "P[ntmax-1000],P[ntmax-500],P[ntmax] ;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 165 "Evidently the value of the \+ pressure is stable in the later parts of the simulation. To display th e time evolution we perform averages over a certain number of steps." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "Pa:=[seq(add(P[(i-1)*20+n u],nu=1..20)/20,i=1..ntmax/20)]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "i:='i': plot([[i,10^(17)*Pa[i]]$i=2..nops(Pa)],style= point);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 178 "Now that we have conf idence that the measurement of the gas pressure on the container walls is reasonably stable we can compare it with the expected value from t he ideal gas law:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "Pr:=Np /Vol*k*T0;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 261 "We are getting a c lose answer from the simulation. It is possible that the truncation of the velocity interval results in a distribution that corresponds to a somewhat lower temperature. In fact, the average speed in the sample \+ fell short of the expected value." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 636 "One should not be too surprised that the simulation comes up with a value for the pressure that agrees with th e ideal gas law. The derivation of the ideal gas law is given in first -year physics texts which emphasize the connection between descriptive thermodynamics and kinetic theory (e.g., Serway, 4th ed., chapter 21) . The relationship between the gas pressure and the average molecular \+ kinetic energy (which is connected to the temperature) follows from a \+ straightforward argument about the momentum transfer for a particle co llision at the wall, and the time it takes to traverse the container b etween two molecule-wall collisions." }}{PARA 0 "" 0 "" {TEXT -1 127 " One has to sum over all molecules (thus, the dependence on the particl e number, and to convert the force impulse to a pressure." }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 572 "The simpl e programming exercise allows one to verify the ideal gas law by simul ation. Note that we have not included collisions between the gas molec ules (hence ideal gas), and that we have not discussed the mechanism b y which the gas acquires and maintains its temperature. In a real sett ing the heating is provided through (some of) the walls. Maintaining t he walls at some temperature implies that the walls are not stationary and impart kinetic energy on the gas molecules during the collisions. Nevertheless, we can verify a few predictions made by the ideal gas l aw." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 311 11 "E xercise 3:" }}{PARA 0 "" 0 "" {TEXT -1 274 "Repeat the present simulat ion with a new ensemble of particles. How consistent is the result for the pressure? Increase/decrease the ensemble size by a factor of two, and observe the effect on the scatter in the results for the gas pres sure for the chosen temperature value." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT 289 11 "Exercise 4:" }}{PARA 0 "" 0 "" {TEXT -1 290 "Repeat the simulation for different values of the temper ature T0. Graph the result for the pressure as a function of T0, and c ompare with the expected result from the ideal gas law. Explain why di screpancies will appear at high T0, unless the algorithm for the sampl e selection is modified." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT 290 11 "Exercise 5:" }}{PARA 0 "" 0 "" {TEXT -1 182 "Chan ge the volume of the container. Repeat simulations at temperature valu es T0 for which results for the pressure were obtained before. Check t he consistency with the ideal gas law." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT 291 11 "Exercise 6:" }}{PARA 0 "" 0 "" {TEXT -1 97 "At constant volume and temperature observe the dependence of the pressure on the particle number." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 292 16 "Pressure balance" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 155 "An inter esting extension of the simulation is to replace the top wall of the c ontainer by a piston which experiences the force of gravity due to som e mass " }{TEXT 293 1 "M" }{TEXT -1 102 ". The educational program pac kage Albert (Springer-Verlag) contains a simulation of this set-up. Th is " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 295 "W hat modifications do we require from our program to accomodate this ch ange? We can imagine that the height of the piston is a direct functio n of the gas pressure: if the external pressure exceeds the gas pressu re, the piston moves down; if the gas pressure exceeds the pressure du e to the mass " }{TEXT 294 1 "M" }{TEXT -1 245 " pushing the piston do wn as a result of gravity, the piston moves up. Eventually the piston \+ reaches a height at which the two pressures are equal (equilibrium). I t is probably useful to add some friction to the equation of motion fo r the piston." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 188 "The main addition to the simulation is the equation of m otion for the piston. This can be handled by a very naive step-by-step integration of Newton's equation for piston height (variable " } {TEXT 19 2 "H0" }{TEXT -1 41 " in the program below) and piston speed \+ (" }{TEXT 19 2 "V0" }{TEXT -1 144 "). In addition one has to modify th e collision time criterion by checking separately for collisions with \+ the variable upper wall (piston height " }{TEXT 308 1 "H" }{TEXT -1 25 " instead of fixed height " }{TEXT 307 1 "L" }{TEXT -1 2 ")." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 122 "The Newt on equation for the piston height depends on the pressure of the gas t imes the surface area of the piston (length " }{TEXT 309 1 "L" }{TEXT -1 355 " in the 2D simulation). When the external mass on the piston i s increased, and the piston moves down, the gas pressure increases on \+ account of the fact that more collisions occur in the vertical directi on as a result of the reduced path length. At some point the gas press ure is sufficient to overcome the downward pressure, and some equilibr ium is found." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 70 "We start the simulation at the end of a previous cycle (p iston height " }{TEXT 310 1 "L" }{TEXT -1 53 "). We have a well-define d value for the gas pressure." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 76 "When choosing the values of parame ters we have to keep a few things in mind:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 51 "1) this is a 2D simulation, and not a real-life one" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 316 "2) our sample of gas molecules is extremely small (100 instead of Avogadro's number); it might be better to think of the par ticles as being bunches of particles, i.e., to increase their mass app ropriately. Alternatively, we can adjust the external mass to be very \+ small [of the order of 100/N_A, where N_A=6*10^(23)]" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 260 "3) the time scale over which simulations need to be run to reach equilibrium may be very lon g. To adjust for this we chose a strong value for the gravitational ac celeration (increase by a factor of 100), and adjusted also the fricti on constant correspondingly." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 276 "The simulation below runs for a very long time in Maple (hours). Simulations are better performed in a compiled floa ting-point environment (Fortran, C, or Matlab). Nevertheless, one has \+ to make compromises, and limit the realism of a simulation due to lack of computing power." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "M:=20*10^(-21);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "b0:=100.; # friction constant for the piston" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "g:=9.81*100; #artificial in crease of gravity" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "Prt:=Pr ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "H0:=L; V0:=0;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "tcoll:=tscale: P[0]:=Prt:" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "ntmax:=Np*50:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 112 "for ntim from 1 to ntmax do: tim:=0: Delv:=0: Are a:=4*L+4*H0: for iit from 1 to 40 do: #step through collisions" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "for i from 1 to Np do: #determine s hortest collision time" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "tx[i]:=(L *sign(vx[i])-x[i])/vx[i]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "if tx[ i] " 0 "" {MPLTEXT 1 0 124 "if vy[i]>0 then ty[i]:=(H0-y[i])/vy[i]; else ty[i]:= (-L-y[i])/vy[i]: fi: if ty[i] " 0 "" {MPLTEXT 1 0 38 "tim:=tim+tcoll: #now update p ositions:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for i from 1 to Np do: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 80 "x[i]:=evalf(x[i]+vx[i]*tcoll): \+ y[i]:=evalf(y[i]+vy[i]*tcoll): od: tcoll:=tscale:" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 53 "#now reverse the momentum (velocity) of colliding T P:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 142 "if tx[icoll] " 0 "" {MPLTEXT 1 0 31 "od: P[ntim]:=m*Delv/(Area*tim):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "H0:=evalf(H0+tim*V0); V0:=evalf(V0+tim*(-g+P[ntim]*L/ M-b0*V0)); Hv[ntim]:=H0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "P[2],P[3],P[4],P[5],P[6],P[4 0],P[100];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "seq(Hv[i],i=1 ..20);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "i:='i': plot([[i,10^(17)*P[i]]$i=2..ntmax ],style=point,view=[0..ntmax,0..40]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "i:='i': plot([[i,Hv[i]]$i=1..ntmax],style=point,view= [0..ntmax,-0.1..0.1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 227 "Please not that the horizontal axis is not labeled by real time, but by the \+ number of collisions. To look at the evolution as a function of time o ne would have to store the time in a table or list (by accumulating th e variable " }{TEXT 19 3 "tim" }{TEXT -1 90 " in the loop), and plot t he pressure and displacement of the piston against this variable." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 340 "The fric tion constant was chosen to be sufficiently large to suppress much lar ge-scale oscillations. It is interesting to observe that despite the e rratic input from the pressure measurements the piston height turns ou t to be a rather smooth function. This is the result of the smoothing \+ properties of the integration of Newton's equations." }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 314 11 "Exercise 7:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 165 "If your comput er runs fast enough (or you can leave it to run unattended for hours) \+ look at the results for different values of the input parameters: (ext ernal mass " }{TEXT 313 1 "M" }{TEXT -1 20 ", friction constant " } {TEXT 312 1 "b" }{TEXT -1 153 "). You can also look at runs for differ ent initial temperatures (in this case the first part of the worksheet has to be re-run with a different value of " }{TEXT 19 2 "T0" }{TEXT -1 29 "). For the present choice of " }{TEXT 19 7 "T0=1000" }{TEXT -1 71 " (Kelvin) we found the approximate equilibrium positions and press ures:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 45 " Mass Height H Pressure P" }}{PARA 0 "" 0 "" {TEXT -1 18 "5*10^(-21) 0.06 " }}{PARA 0 "" 0 "" {TEXT -1 43 "1*10^( -20) 0.0 1*10^(-16)" }}{PARA 0 "" 0 "" {TEXT -1 42 " 2*10^(-20) -0.04 2*10^(-16)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 103 "Verify that pressure time volu me (area of the rectangle in the 2D case) remains approximately consta nt." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 20 }{VIEWOPTS 0 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }