\caption{van der Waals の状態方程式(p vs. V):
\ 理想気体の場合、He,Arの場合をあわせて表示.T=200K.}
\label{fig:vdw1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
コンピュータ演習 : プログラムの作成
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{はじめに}
現代のコンピュータの利用法には、大きく分けて次の二つがあると言える。
* 数値計算
* 情報処理
化学の分野では、日々新たに作られていく化合物を対象としなければなら
ないので、後者の意味のコンピュータ利用はたいへん重要である。しかし
ここでは前者についてとりあげ、しかも計算のためのプログラムを自分で
作ることを実習する。将来出来合いのプログラムパッケージを利用する場
合にも、プログラムの作成を経験していれば、ブラックボックスの中身に
ついての感覚を持ち合わせていることになり、応用の幅が広がるであろう。
化学の研究で行われている数値計算には次のようなものがある。
* 実験データの処理
+ 統計的な処理
+ 仮説・モデルへの当てはめ
* 計算機実験・シミュレーション(simulation)
+ 分子軌道計算
+ 分子力場計算
+ 物性値の計算
+ 仮説・モデルによる予測
ここで実験データの処理については、直線回帰の様な定型的なものであれ
ば、現在では表計算ソフトを用いるのが便利である。分子軌道計算などは、
初歩的なものをのぞいて通常はパッケージのプログラムを用いる。その他
のモデルへの当てはめ、予測や、シミュレーションは、個別の問題であり、
研究者がそれぞれ独自に開発したプログラムを用いることが多い。
ここではモデルによる予測に近いものとして、理論式による物性値の計算
をおこない、それを視覚的にわかりやすく図示することを試みる。そのた
めにまず、物理化学の教科書に出てくるような理論式をとりあげる。そし
てこれを用いた計算のプログラムの実例を示して、その流れを解説する。
そのあと別の式による計算のプログラムの作成を、演習としておこなう。
サンプルプログラムは、共通教育の情報処理の簡単なプログラム作成で用
いたc言語により記述する。\footnote{%
B. W. カーニハン, D. M. リッチー 著, 石田晴久 訳, プログラミング言語C, 第2版,
共立出版. }
なお簡単化のために、計算のプログラムは結果をファイルに書き出すものとする。
実際にグラフにプロットするのは、表計算ソフトを用いるものとする。
課題の部分はレポートとして提出すること。なおこの場合に、プログラム
を記述する言語は、各人が使用可能なものを用いれば良い。(c言語でのプ
ログラム作成を強制するものではない。)
\subsection*{注意点}
* 3.5-inch のフロッピーディスク(1.44MB または DOS/V 用と書かれている
もの)を各人が用意し、データファイルはそこに保存すること。(共有
ディスクを用いることもできる。)パソコンのハードディスクには
save しない。
* 昨年度の共通教育:情報処理(理系)の実習用資料をもってくること。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{実在気体の状態方程式 (van der Waals の方程式)}
実在の気体についての状態方程式はvan der Waalsにより次の\ref{eq:vdW}式
のように与えられている。\footnote%
{P. K. Atkins, {\it Physical Chemistry}, 5th ed., Oxford,
p. 44 \& Data section C4.}
p = (nRT/(V-nb))-a(n/V)^2
ここでパラメータ a,b は、表\ref{tbl:vdW}に示すとおりである。
\begin{table}[h]
\caption{van der Waals coefficients.}
\label{tbl:vdW}
-----------------------------------------------------
& a/( atom L^2 mol^{-2}) & b/(10^{-2}L mol^{-1})
-----------------------------------------------------
Ar & 1.345 & 3.220
CO_2 & 3.592 & 4.267
He & 0.034 & 2.37
N_2 & 1.390 & 3.913
-----------------------------------------------------
\end{table}
van der Waals の実在気体の状態方程式による、実在気体の体積と圧力
との関係を計算するプログラムの例を図\ref{prg:vdw}に示す。またこの
プログラムの実行結果を gnuplot で表示させた結果を図\ref{fig:vdw1},
\ref{fig:vdw2}に示す。これらのグラフは、 a=0, b=0 の理想気体の
場合と、He,Ar の場合について計算し、その3つの結果を重ねて表示し
てある。図\ref{fig:vdw1}では少しわかりにくいが、図\ref{fig:vdw2}で
は圧力が高くなるほど理想気体からずれていく様子がよくわかる。
\begin{figure}
\caption{van der Waals の実在気体の状態方程式による体積と圧力との関係を
計算するプログラム.}
\label{prg:vdw}
#include
#include
void main( int argc, char* argv[] ){
float T, a, b;
float n, R;
float p, V;
if ( argc!=4 ) {
fprintf(stderr, " parameter error \n");
fprintf(stderr, " usage: vdw T a b \n");
exit(1);
}
T=atof(argv[1]);
a=atof(argv[2]);
b=atof(argv[3]);
fprintf(stderr, " T/K=%f, a/(atom L^2 mol^-2)=%f, b/(L mol^-1)=%f \n",
T, a, b);
R=8.31451; /* J mol^-1 K^-1 */
n=1.0; /* mol */
for(V=0; V<10; V+=0.01) {
if (V!=0) {
p=(n*R*T)/(V-n*b)-a*(n/V)*(n/V);
if (p>0)
fprintf(stdout, "%f %f %f\n", p, V, 1/V);
}
}
}
\end{figure}
\begin{figure}
\caption{van der Waals の状態方程式(p vs. V):
\ 理想気体の場合、He,Arの場合をあわせて表示.T=200K.}
\label{fig:vdw1}
\caption{van der Waals の状態方程式(p vs. 1/V).}
\label{fig:vdw2}
\end{figure}
%%%%%%%%%%%%%%%%%%%%
\subsection*{プログラムの解説}
以下に図\ref{prg:vdw}のプログラムについて解説する。解説の便宜上各行
に番号がついているが、実際にはこの番号は必要ない。
[line 1-2] プログラム中で用いる標準関数のプロトタイプや定数の宣言の
ファイルをここに読み込む。始めのうちは c 言語の約束の呪文だと思って
もよい。詳しくは解説書、文法書を参照のこと。
[line 4] このプログラムは、メインルーチンだけからなる。起動時にコマ
ンドラインから指定した引数を使用する。
[line 5-7] プログラム中で用いる変数の型を宣言する。ここではすべて浮
動小数点による実数型である。
[line 9-13] 起動時にコマンドラインから指定した引数の数がおかしければ、
簡単な使い方の説明とあわせて表示し、プログラムを終了する。
fprintf 関数をもちいてメッセージを標準エラー出力に表示する。
exit 関数はプログラムを終了し、引数を OS に返す。
[line 15-17] 起動時にコマンドラインから指定した引数を変数に代入する。
[line 18] きちんと値が代入されたかどうかの確認のため、変数の値を標
準エラー出力に表示する。
[line 21-22] プログラム中で使用する物理定数の値を変数に代入する。
n と R は定数変数としても良いが、このプログラムでは通常の変数にし
てある。
[line 24-30] 体積 V を変えながら、その時の圧力 p を計算する。計算は
V \neq 0 の時にのみおこない、また p > 0 の時にのみ結果を標準出力に
表示する。
プログラムの実行時に、コマンドラインで標準出力をリダイレクトすれば、
計算結果だけをファイルに得ることができる。これを表計算ソフトなどで
図示することができる。
[line 26] van der Waals の実在気体の状態方程式により、気体の圧力を
計算する。
[line 28] 計算結果を標準出力に表示する。p,V だけでなく、1/V も出力する。
このプログラムは、25行と27行の if 文でわかるように、特定の条件
を満たした場合にのみ、計算および結果の出力を行うようになっている。
式の上ではどんな計算も可能であるが、意味のある結果であるかどうかの判
断は、利用する人間が考えなければならない。計算機で求めた結果であるこ
とは、結果の正しさとは無関係である。
\begin{exam}
図\ref{prg:vdw}のプログラムを入力し、動作を確かめよ。\footnote{%
backslash と yen の記号は表示の上では異なるが、機械内部のコードとしては
同じものである。適宜読み換えること。}
\end{exam}
\begin{exam}
25行と27行の if 文による制限を取り払うとどんな結果になるか。
\end{exam}
\begin{exam}
実行時にコマンドラインから指定する引数として、表\ref{tbl:vdW}に示す値
とは異なるものを試してみよ。どのような結果が得られるか。
\end{exam}
\begin{exam}
教科書から適当な式を拾いだし、そのグラフを描くプログラムを作ってみよう。
\end{exam}
\begin{prob}
球面調和関数(Spherical Harmonics)の形を、極座標で描け。\footnote{%
P. K. Atkins, {\it Physical Chemistry}, 5th ed., Oxford, p. 413.}
2p, 3d 軌道からそれぞれひとつずつ描け。興味があれば 4f 軌道についても
試みよ。これは例えば y=0 の xz 平面において(\phi=0)、\theta を
0--2\pi まで変化させた時の関数の値の絶対値を、その方向での原点
からの距離に対応させればよい。これは xz 平面による断面図にである。
形がわかれば良いので、規格化定数は省略してよい。また数学関数を用い
るので、math.h を #include すること。
\end{prob}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{群の表現行列の積}
分子の性質を調べる場合に、分子の形を考えにいれておくと、わかりやすい
場合が多い。この様な分子の形をとりあつかう学問分野に、
``群論(group theory)''というものがある。\footnote%
{P. K. Atkins, {\it Physical Chemistry}, 5th ed., Oxford,
Chapter 15, pp.511--538.}
群論では分子をある軸の回りでまわしたり、ある面で鏡映させたときに、
もとの分子とぴったり重なるような操作(対称操作)について考える。
回れ右してから鏡に映すと、その結果がどうなるかを考えるわけである。
この様な対称操作は、数学的には行列をもちいて記述でき、上の例では
それぞれの操作を表す行列の積を求めるのである。
例えばアンモニア分子(NH_3)のようなものは C_{3v} の点群に属すると言い、
ここで考えられる対称操作は、
動かさない(E)、N を通る軸で 120度回す(C_3)、同じ軸で 240度回す(C_3^2)、
N-H のひとつを含む面で鏡映する( 3通り、(\sigma_v,\sigma_{v'},\sigma_{v''}))、
の合計 6通りのものがある。これらを 4 x 4 の行列で表すと
表\ref{tbl:c3v}の様に書くことができる。
またこれらの対称操作の積は表\ref{tbl:mult}の様になる。
ここでは、Firstの操作をあらわす行列が右に、Secondの方が左にくる。
(New\_operation) = (Second)(First)
\begin{table}
\caption{The matrix representation of C_{3v}.}
\label{tbl:c3v}
-------------------------------------------------
E & C_3 & C_3^2
( 1 0 0 0 ) ( 1 0 0 0 ) ( 1 0 0 0 )
( 0 1 0 0 ) ( 0 0 0 1 ) ( 0 0 1 0 )
( 0 0 1 0 ) ( 0 1 0 0 ) ( 0 0 0 1 )
( 0 0 0 1 ) ( 0 0 1 0 ) ( 0 1 0 0 )
-------------------------------------------------
\sigma_v & \sigma_{v'} & \sigma_{v''}
( 1 0 0 0 ) ( 1 0 0 0 ) ( 1 0 0 0 )
( 0 1 0 0 ) ( 0 0 1 0 ) ( 0 0 0 1 )
( 0 0 0 1 ) ( 0 1 0 0 ) ( 0 0 1 0 )
( 0 0 1 0 ) ( 0 0 0 1 ) ( 0 1 0 0 )
-------------------------------------------------
\end{table}
%%%%%
\begin{table}
\caption{The multiplication table of C_{3v}.}
\label{tbl:mult}
----------------------------------------------------------------
First & E & C_3 & C_3^2 & \sigma_v & \sigma_{v'} & \sigma_{v''}
----------------------------------------------------------------
Second
E & E & C_3 & C_3^2 & \sigma_v & \sigma_{v'} & \sigma_{v''}
C_3 & C_3 & C_3^2 & E & \sigma_{v'} & \sigma_{v''} & \sigma_v
C_3^2 & C_3^2 & E & C_3 & \sigma_{v''} & \sigma_v & \sigma_{v'}
\sigma_v & \sigma_v & \sigma_{v''} & \sigma_{v'} & E & C_3^2 & C_3
\sigma_{v'} & \sigma_{v'} & \sigma_v & \sigma_{v''} & C_3 & E$ & C_3^2
\sigma_{v''} & \sigma_{v''} & \sigma_{v'} & \sigma_v & C_3^2 & C_3 & E
----------------------------------------------------------------
\end{table}
%
n次元の行列 A と B があって、それぞれの(i,j)要素が a_{ij},b_{ij}
であるとする。このときの積 C=AB の(i,j)要素 c_{ij} は、
c_{ij} = \sum_{k=1}^{n} a_{ik}b_{kj}
と書くことができる。行列の積を計算するプログラムは、これをそのまま
記述すれば良い。プログラムを図\ref{prg:mat1},\ref{prg:mat2}に示す。
\begin{figure}
\caption{ 4 x 4 行列ふたつの積を計算するプログラム.
メインルーチンだけのファイルにする.}
\label{prg:mat1}
#include
#include
void input_mat44( float m[4][4] );
void print_mat44( float m[4][4] );
void main()
{
int i, j, k;
float w;
float A[4][4], B[4][4], C[4][4];
fprintf(stderr,"Input matrix A \n");
input_mat44( A );
fprintf(stderr,"Input matrix B \n");
input_mat44( B );
fprintf(stderr,"Print matrix A \n");
print_mat44( A );
fprintf(stderr,"Print matrix B \n");
print_mat44( B );
for(i=0; i<4; i++) {
for(j=0; j<4; j++) {
w=0;
for(k=0; k<4; k++) w+=A[i][k]*B[k][j];
C[i][j]=w;
}
}
fprintf(stderr,"Print resulting matrix C=AB \n");
print_mat44( C );
}
#include "mat2.c"
\end{figure}
\begin{figure}
\verbatimlisting{mat2.c}
\caption{ 4 x 4 行列ふたつの積を計算するプログラムの続き.
``mat2.c''というファイル名にすること.}
\label{prg:mat2}
void input_mat44( float m[4][4] )
{
int i, j;
for(i=0; i<4; i++)
for(j=0; j<4; j++) {
fprintf(stderr, "m[%i][%i] = ", i,j);
scanf("%f", &(m[i][j]));
}
}
void print_mat44( float m[4][4] )
{
int i, j;
for(i=0; i<4; i++) {
for(j=0; j<4; j++) fprintf(stderr, "%f ", m[i][j]);
fprintf(stderr, "\n");
}
}
\end{figure}
%%%%%%%%%%%%%%%%%%
\subsection*{プログラムの解説}
以下に図\ref{prg:mat1}のプログラムについて解説する。解説の便宜上各行
に番号がついているが、実際にはこの番号は必要ない。また、テキスト作成
の都合でプログラムは図\ref{prg:mat1}と図\ref{prg:mat2}のふたつのファ
イルに分割されており、後者は mat2.c というファイル名である必要がある。
もちろんひとつのファイルに結合しても良いが、その場合には図\ref{prg:mat1}の
34行目の #include 文は削除する。
[line 1-2] プログラム中で用いる標準関数のプロトタイプや定数の宣言の
ファイルをここに読み込む。
[line 4-5] 用いているサブルーチンのプロトタイプ宣言。
[line 7] このプログラムは、処理の一部をサブルーチンにしているが、中
心となる処理はこのメインルーチンに書いてある。サブルーチンは
図\ref{prg:mat2}に示す。
[line 9-11] メインルーチン中で用いる変数の型を宣言。
[line 13-16] 4 x 4 行列を入力する。実際の入力はサブルーチン
input_mat44 で行う。
[line 17-20] 確認のため、入力された行列を表示する。実際の表示は、サ
ブルーチン print_mat44 で行う。
[line 22-28] 行列の積の計算。c_{ij} がいくつになるか、順次計算する。
[line 29-30] 計算結果の行列を表示する。
[line 34] テキスト作成の都合でサブルーチン部分が別ファイルになってい
るが、それをここに読み込む。分割コンパイルはしない。
\begin{exam}
表\ref{tbl:c3v}の行列どうしの積を計算し、表\ref{tbl:mult}の様になって
いることを確認せよ。
\end{exam}
\begin{exam}
その他のデータを入力し、行列の積が計算できることを確認せよ。
\end{exam}
\begin{prob}
逆行列をもとめるサブルーチンを用い、
4元連立一次方程式を解くプログラムを作れ。
\end{prob}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{課題}
ここまでに出てきた課題とあわせて、以下の課題からはひとつ以上を選んで
プログラムを作成し、レポートとして提出せよ。レポートはプログラムのソ
ースリスト、プログラムの流れ(計算手順、アルゴリズム)の簡単な解説、プ
ログラムの実行結果、目的の図表とそれを得るための手順など、および今回
のプログラム作成に関する感想などを含み、ワードプロセッサで作成すること。
さらにこれらのデータは、各人のフロッピーディスクに記録し、氏名を書い
たラベルをはって提出すること。協力した友人、参考図書などがあれば、も
ちろん記すこと。
\begin{prob}
n次の行列の逆行列を求めるプログラムを作れ。行列の次数はデータで指定
できるようにせよ。逆行列を求めるアルゴリズムとしては、(1)ガウスの消
去法を用いる(掃き出し法)、(2)行列の余因子を用いる、(3)固有値分解を用
いる、等の方法が考えられるが、どんなアルゴリズムを用いたかの説明をつ
けること。また、入力データを確認し、間違いがあれば訂正できるような機
能を持たせよ。逆行列がないような行列を、入力データとして与えた場合に
は、正しく処理を中止して、その旨表示するようにすること。
\end{prob}
\begin{prob}
次に挙げる関数を x=0 の近傍で Taylor 展開し、高次の項までとるにつれて
近似が良くなる様子を図示せよ。
(1) $\sin(x)$; (2) $\exp(-x)$; (3) $1/x$
\end{prob}
\begin{prob}
次に挙げる関数の積分(定積分)を計算し、区間の分割の程度と計算手法に
よる違いや、解析的な値との差を調べ、必要なら図示せよ。
(1) $\sin(x),~[0\sim\pi]$; (2) $\exp(-x),~[0\sim1]$; (3) $1/x,~[1\sim2]$
\end{prob}
\begin{prob}
表\ref{tbl:pvt},\ref{tbl:pvt2}\ に与えるデータを直線にあてはめ、
データ点と直線を図示せよ。また、表\ref{tbl:wtph}\ に与えるデータを
2次曲線にあてはめるとどうなるか。
\begin{table}[h]
\caption{Data1} \label{tbl:pvt}
---------------
t & p
---------------
3.62 & 0.96541
3.59 & 0.91373
3.53 & 0.80672
3.52 & 0.78949
3.49 & 0.76586
3.48 & 0.73067
3.46 & 0.70854
3.45 & 0.67210
3.43 & 0.64124
3.41 & 0.62337
3.40 & 0.59077
3.38 & 0.56463
3.37 & 0.54156
3.35 & 0.51327
3.34 & 0.49933
3.33 & 0.48062
3.32 & 0.46218
3.28 & 0.39637
3.27 & 0.37702
3.26 & 0.35421
3.23 & 0.32186
3.22 & 0.29569
----------------
\caption{Data2} \label{tbl:wtph}
------------
w & t
------------
10 & 45.3
15 & 58.9
20 & 64.0
25 & 66.5
30 & 67.1
36 & 67.2
40 & 67.0
45 & 66.2
55 & 61.0
64 & 48.0
70 & 31.0
------------
\end{table}
\end{prob}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{参考図書}
* プログラム作成に当たっての考え方
+ B. W. カーニハン、P. J. プローガー \ 著、木村泉訳、
「プログラム書法 \ 第2版」、共立出版
+ B. W. カーニハン、P. J. プローガー \ 著、木村泉訳、
「ソフトウェア作法」、共立出版
* アルゴリズムに関するもの
+ N. ヴィルト \ 著、浦昭二・國府方久史 \ 共訳、
「アルゴリズムとデータ構造」、近代科学社
+ 奥村晴彦、「C言語による最新アルゴリズム辞典」、技術評論社
* C言語に関するもの
+ B.W.カーニハン、D.M.リッチー \ 著、石田晴久 \ 訳、
「プログラミング言語 C \ 第2版」、共立出版
+ 「改訂第3版 \ C言語入門」、アスキー
* 数値計算に関するもの
+ 高橋大輔、理工系の基礎数学 \ 8 \ 「数値計算」、岩波書店
+ 森正武、「FORTRAN77 \ 数値計算プログラミング \ 増補版」、岩波書店
+ 小池慎一、「改訂 \ Cによる科学技術計算」、CQ出版
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{付録:Exelでグラフを描く}
数値計算の結果をファイルに出力したら、その次にはこの結果をグラフに
することを考える。このテキストでは図\ref{fig:vdw1},\ref{fig:vdw2}に
示したようなグラフを、gnuplot というソフトを用いて描いた。その他に
もいろいろな方法があるだろうが、ここでは表計算ソフト(Exel)のグラフ
表示機能を利用する。
Exelによるグラフの書き方を、手順を追って示す。なおここで、グラフ化
しようとするデータは、PC98 の B: ドライブにあるフロッピーディスクの
datafile という名前のファイルであるとする。データはx軸の値とy軸の値
がこの順で一行にひと組ずつ記述されており、数値は空白で区切られてい
るものとする。
1) Windows3.1のExelをダブルクリックして起動する
2) Exelで、datafileを読み込む
a) メニューバーの[ファイル]--[開く]を順に選ぶ
b) ドライブ(b:)、ディレクトリ(\yen)、ファイル名(datafile)を順に指定する
3) テキストファイルウィザードが現れ、ファイルの形式を指定する
a) 「区切り文字によるデータ並び」の形式を選ぶ
b) 「列の区切り文字はスペース」をチェック
c) データ形式は「標準」を選ぶ
4) ツールバーのグラフウィザードポタンを押す
a) データ範囲を列ごとに2回指定する
i) まずx軸のデータを指定する。$A$1 から縦に $A$9 までマウスを
ドラッグすると、ウインドウに =$A$1:$A$9 の文字が入っている
ii) カンマ(,)を入力する
iii) y軸のデータとして、$B$1 から縦に $B$9 までマウスをドラッグする
iv) ウインドウに =$A$1:$A$9,$B$1:$B$9 の文字が入っている
b) グラフの種類として、「散布図(S)」を選ぶ
c) グラフの形式を選ぶ。マークを打って、折れ線で結ぶのなら、(2)を選ぶ。
d) グラフのサンプルが表示されるので、「先頭の1列をxの値として
使用」すると設定する
e) 必要ならグラフのタイトルや軸のラベルを設定する
5) 画面上に枠が表示され、グラフが描かれる。枠をマウスでドラッグす
ると縦横好きな大きさに引き延ばされる
6) 枠の中程でマウスをドラッグすると、画面上でグラフの位置を動かす
ことができる
7) 出来上がったら印刷してみよう
a) この方法ではワークシート上にグラフが乗った形になっているので、
グラフの位置などを確認し、必要なら直す
b) [ファイル]--[プレビュウ]で出来上がりを確かめる
c) ([ファイル]--)[印刷]すると、ディスプレイの上部に張ってある色
のシールと同じプリンタから出力される
d) 紙をむだにしないように
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%