#topicpath ** van der Waals の状態方程式について [#y98a85e9] -- Maxima で計算 --- 状態方程式; P(V,RT,N,a,b) := N*RT/( V- b*N) - a * N^2/V^2; := で関数を定義している。関数の引数は複数でもよい。~ ~ ---「臨界状態の V, RT, P」 ''を'' 「a, b」 ''で'' 表現; result1 : solve( [diff(P(V,RT,N,a,b),V)=0, diff(diff(P(V,RT,N,a,b),V),V)=0], [V,RT])[1]; solve() の部分で、∂P/∂V=0 と、∂^2P/∂V^2=0 $ \left( \frac{\partial p}{\partial V}=0 , \frac{\partial^2 p}{\partial V^2}=0 \right) $ の両方を満たす V, RT を求めている。~ 連立方程式は [ ] 内で、"," でつないで書き連ねる。解く変数も同様に[ ] 内に書く。~ 解いた結果は result1 に代入される。: は代入を表す。~ ~ P(rhs(result1[1]), rhs(result1[2]), N, a, b); result1 には複数(V, RT)の変数=値という式が並んでいる。その1番めは result1[1] で取り出し、右辺の値は、rhs() で取り出している。~ それを P() に代入している。~ ~ ---「臨界状態の V, RT」 ''で'' 「P, a, b」 ''を'' 表現; result2 : solve( [diff(P(V,RT,N,a,b),V)=0, diff(diff(P(V,RT,N,a,b),V),V)=0], [a,b])[1]; P(V, RT, N, rhs(result2[1]), rhs(result2[2])); -- gnuplot で描画 # # van der Waals の状態方程式を # 臨界点での温度と体積で規格化する。 # # パラメタ # N = 1.0 RTc = 1.0 Vc = 1.0 a = 9.0*RTc*Vc/(8.0*N) b = Vc/(3.0*N) set yrange [0:1] set xrange [0.4:5] # # 状態方程式 PR:実在気体 PI:理想気体 # PR(V,RT) = N*RT/(V-b*N) - a*N**2/V**2 PI(V,RT) = N*RT/V # # 確認のため表示 # ( gnuplot で整数の割り算には注意。truncation が起こる) # print 1/3 print PR(1,1) # # X11 # set term x11 0 plot PR(x,2), PR(x,1.5), PR(x,1.2), PR(x,1.1), PR(x,1.0), PR(x,0.90) set term x11 1 plot PI(x,2), PI(x,1.5), PI(x,1.2), PI(x,1.1), PI(x,1.0), PI(x,0.90) # # gif # # set term gif small x808080 x909090 xA0A0A0 xB0B0B0 xC0C0C0 xD0D0D0 ; set output "VanDerWaals_PR.gif" set term gif small; set output "VanDerWaals_PR.gif" plot PR(x,2), PR(x,1.5), PR(x,1.2), PR(x,1.1), PR(x,1.0), PR(x,0.90) set term gif small; set output "VanDerWaals_PI.gif" plot PI(x,2), PI(x,1.5), PI(x,1.2), PI(x,1.1), PI(x,1.0), PI(x,0.90) quit