数値計算プログラミング

Visual Studio Code

　Visual Studio Codeを用いてC言語の学習をするための準備をします。Visual Studio Code(以下、VSCode)は、Microsoftが提供しているフリーのエディタです。初心者にも扱いやすいと思います。また、Mac、Linux、Windowsで使用可能なので大変便利です。

ダウンロード

　ダウンロードは こちら から。各自のOSに合わせてダウンロード、インストールしてください。

新規ファイルを作成する

　プログラムが書き込んであるファイルをソースファイルと呼びます。今回は、初めてプログラムを作成するときによく使われる「Helloworld」のプログラムを例に実行まで説明します。まず、空のソースファイルを作成しましょう。

注釈

拡張子の表示・非表示 拡張子とは、ファイルの種類を識別するための文字列です。ファイル名には必ず拡張子があり、ドット以下に記述されます。普段あまり意識することはないかもしれませんが、プログラムを扱う上では拡張子は重要になります。拡張子が非表示になっている場合には、リンクの記事をご覧ください。 [Windows]拡張子の表示・非表示

  1. VSCodeのファイルメニューより、新規ファイルを選択します。そうすると、新規ファイルへの記述が可能になります。

  2. [名前をつけて保存]を選択します。このとき、ファイルの拡張子(.c)をつけるのを忘れないでください。

  3. 保存が完了したら、そのファイルへの記述が可能になります。その後は定期的に保存をするようにしましょう。

  4. 「Helloworld」のプログラムは以下の通りです。コピーすることもできますが、せっかくなので手で打ってみることをお勧めします。

#include<stdio.h>

int main(void){
  printf("Helloworld\n");

  return 0;
}

拡張機能のインストール

VSCodeには、豊富な拡張機能があります。今回はC言語用の拡張機能をインストールしましょう。

image

まず、拡張機能のパネルを表示します。赤い丸で示されたところをクリックします。

検索フォームに"C/C++"と入力します。一番上に出てきた拡張機能をインストールします。インストールが終了したら、有効化するために[再度読み込む]をクリックしてください。

image

以上で、拡張機能のインストールは終了です。

Microsoft Visual Studio

準備

　Visual Studio 2008 で簡単な C言語プログラムを作成する手順を説明します。まず、メニューからファイル\(\rightarrow\)新規作成\(\rightarrow\)プロジェクトを選びます。そこで、プロジェクトの種類、テンプレートを

  • プロジェクトの種類：全般

  • テンプレート：空のプロジェクト(NET)

のように指定して、プロジェクト名を適当につけます (プロジェクト名はフォルダ名になります)。「場所」は変更する必要がないので、そのまま「OK」を選びます。すると、左の「ソリューション エクスプローラ」にプロジェクトが表示されます。例えばプロジェクト名を sample とした場合、sample以下にソースファイル、ヘッダーファイル、リソースファイルが作られます。更に、この「ソリューション エクスプローラ」の「ソースファイル」の上で右クリックします。メニューから「追加」\(\rightarrow\)「新しい項目の追加」を選びます。カテゴリーは「Visual C++」となっており、右のテンプレートから「テキスト ファイル」を選択し、ファイル名として例えば「sample.c」とします。

ソースファイルの作成

　ソースファイルに次のようなテキストを入力します。

#include <stdio.h>

void main(){
      printf("Hello world!\n");
}

プログラム説明

　上で作成したソースファイルの簡単な説明をします。

  • #include <stdio.h>

    C 言語から標準入出力 (キーボードからの入力, ディスプレーへの出力など)の機能を使うには、標準入出力用のヘッダーファイル stdio.h (standard IO header = 標準入出力ヘッダー, IO = input and output)の組み込みをします。

  • main(){  }

    C言語のメイン部分です。これより前の部分をヘッダーと呼び、通常は各種のヘッダーファイルを組み込みます。{} で囲まれた部分に動作させたいプログラムを書き込みます。

  • printf( );

    printf(= print format) はprintf("\(\sim\)")の構文で使用します。"\(\sim\)"の部分に表示させたい文章を書き込みます。\nは改行(new line) を意味し、文章を出力したあとで改行することを指定しています。またC言語のプログラムでは、各行の最後には必ずセミコロン;が必要となります。

ビルド

　プログラムは「Hello world!」を出力するためのものです。ソースファイルを作成したらメニューの「ビルド」\(\rightarrow\)「sample のビルド」を選びます(プロジェクト名が sample の場合)。ビルドするとソースファイルから実行ファイルが作成されます。ビルドに成功すればウィンドウの下の部分 (出力用のウィンドウ) に正常終了の表示が出ます。もしもソースファイルにエラーがあればエラーメッセージが表示されます。この場合には、ウィンドウのエラーを表示している箇所をダブルクリックすると、エラーの行を指示してくれます。(実際はこの行、あるいは直前の行にエラーの原因があります。)原因を考えて修正を施し、再度ビルドします。ビルドに成功したら、メニューの「デバッグ」\(\rightarrow\)「デバッグなしで開始」を選んで実行してみましょう。（計算機プログラムのエラーのことを「バグ」(bug = 虫) という。昔、計算機が大きな機械であったとき、動作が不良となって調べたところ回路に「虫」がまぎれこんでいたことが言葉の由来である。プログラムのバグ (虫) をとり、正常に動作するようにすることを「デバッグ」(debug)とか「虫取り」と言う。）

コンパイルとリンク

C言語のソースファイルから実行ファイルを作成するには、コンパイルとリンクをします。コンパイル (compile) はソースファイルの文法チェックと機械語への翻訳をします。この際には、例えばprintfのような関数の手順は読み込みません。コンパイルに成功すれば各種関数の手順を書いたファイル(= ライブラリ) との結合をします。これをリンク (link)といいます。Visual C++ では、ビルド (build) することでコンパイルの後リンクが行われます。

基本的な文法の説明

四則演算
算術演算子と意味

演算子

種別

例

意味

たし算

x+y

xにyを加える。

ひき算

x+y

xからyを引く。

\(\star\)

かけ算

x*y

xにyをかける。

/

わり算

x/y

xをyで割る。

%

余

x % y

xをyで割った余りを求める。
forループ

C言語で、繰り返し処理を行いたい時にfor文を使用します。書式は以下のとおりです。

for (i=0;i<n;i++) {
          /* 処理 */;
          /* 処理 */;
}
定数

C言語において、定数（ていすう、constant）とは、プログラム中で変化することのない一定の値を持つデータのことです。例えば、プログラム中で「100」「3.14」「’A’」「"Hello"」等と、データを直接書いた場合、定数となります。また、定数はリテラルとも呼ばれます。定数式とは定数のみからなる式のことです。

#define START 2
ファイル

fopen()はfilenameで指定されたファイルをオープンし、ファイル・ポインタを返します。エラーが発生した場合はNULLが返されます。ファイル・ポインタはファイルに関する様々な情報を含む構造体へのポインタですが、使用に際して、その内容を知る必要はありません。

FILE *fopen(const char *filename, const char *mode);
型の定義

C言語において、変数（へんすう、variable）とは、処理を行うために必要なデータを記憶しておく領域に名前をつけたものを意味します。名前がついているので、たくさんの種類のデータを扱う場合でも、簡単に区別がつきます。 変数を使用するためには変数の宣言をする必要があります。変数の宣言で最も簡単な方法は、変数の型と変数名を記述するような以下の書き方です

変数の型 変数名;
配列

C言語における配列とは、同じ型のデータの集まりで、それらのデータへ変数名に添え字をつけることによりアクセス可能となるものです。関連のあるデータが複数ある時に、それらに対して一つ一つ変数を割り当てていては非効率的であるという場合に使用されます。添え字には任意の整数の式を使用することができるので、ループ処理を使って一括処理を行うのに便利です。C言語の配列の添え字は0から始まります。他の言語では1から始まるものもありますので、混同しないように気をつけましょう。

int data[10];

熱伝導方程式

一次元の熱伝導方程式

_images/netudendou_gaiyou.png

長い棒の一端を温めたときに、熱がどのように伝わって、温度分布が時間変化していくかという問題を、プログラミングを使ってシミュレーションしましょう。

image

一次元の熱伝導方程式とは上のような式になります．簡単のため、ここではa=1の場合について考えます。ここで問題になるのは“微分”です。微分を差分というかたちにして、熱伝導方程式を書き換えましょう。

　差分法（さぶんほう、difference method）とは、微分方程式を解く数値解析の方法のひとつである。関数が2つの変数値に対してとる値の間の有限な差を差分（さぶん、difference）といい、この差分を変数値の差で割って得られる商を差分商（さぶんしょう、difference quotient）という。微分を差分商で近似することにより微分方程式を解くものである。http://ja.wikipedia.org/wiki/差分法）

問題

　以下に示すテイラー級数展開を利用して、無次元化した方程式を差分化してみましょう。

\[f(x + \Delta x)=f(x)+f'(x)\Delta x + \dfrac{f''(x)}{2}(\Delta x)^{2}+\dfrac{f'''(x)}{6}(\Delta x)^{3} +\dfrac{f''''(x)}{24}(\Delta x)^{4} + ...\]

プログラミング

　方程式が差分化できたら、次は具体的なプログラミングにうつりましょう。

　データ温度分布の時間変化を下の表のようにして考えます。

_images/netudendou_hairetu.png

　初期状態であるt=0のときは、x=0のみ温度の値は1、他は全て0にします。xは0 ~ 100まで考えます.

　これが熱伝導方程式に従う場合、t=nやt=n+1でどういった値をとるかを考えます。

陽解法の説明

　普通に差分化して式を変形すると、下の式になります．この式を使って解を求める方法を陽解法といいます．

_images/explicit.png

　この式をよくみていくと、時刻n+1、場所iのポイントの温度（下の表の赤いセル）を求めるには、時刻nのi-1、i、i+1の三つの場所の温度が必要になることが分かります．

_images/netudendou_explicit.png

　すなわち、前の時刻の温度分布が分かっていれば、それらを使って次の時刻の温度が分かることになります．もっというと、初めの時刻t=0の温度分布が与えられれば、t=1の温度、t=2の温度という具合に、次々と計算できることが分かります。

　この方法は直感的には大変分かりやすく、プログラミングしやすい方法です。

陽解法ソースファイル

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>

#define IMAX    101                                 /* x方向の計算領域 */
#define TMAX    10000                               /* 計算終了時間設定*/
#define dt  0.000001                                /* 時間刻み */
#define dx  0.01                                    /* x方向の格子間隔刻み */

void main(void){
    FILE *fp;                                       /*ファイルのため*/
    int i,t;                                        /* 時間と温度の更新のための整数 */
    double time,T[IMAX],Tnew[IMAX];                 /* 無次元時間と新旧無次元温度 */

    if((fp=fopen("result.dat","w"))==NULL){         /*保存用ファイルを開く*/
        printf("Can not creat file.\n");
        exit(-1);
    }


    /* 境界条件 */
    T[0]=1;Tnew[0]=1;T[IMAX-1]=0;Tnew[IMAX-1]=0;    /* 両境界とも等温条件 */

    /* 初期条件 */
    for(i=1;i<IMAX;i++){                            /* 最初は無次元温度０*/
        T[i]=0;
    }

    time=0;                                         /* 無次元時間の初期化*/

    for(t=0;t<TMAX;t++){                            /* 設定した時間まで計算 */
    time=time+dt;                                   /* 時間を更新*/
        for(i=1;i<IMAX-1;i++){                      /* 境界以外を陰解法の式で計算 */
            Tnew[i]=T[i]+dt/dx/dx*(T[i+1]+T[i-1]-2*T[i]);
        }

        for(i=0;i<IMAX;i++){                         /* 得られた温度を古い温度として置き換え*/
            // fprintf(fp,"%lf,%d,%lf\n",time,i,T[i]);
            T[i]=Tnew[i];
        }
    }

    for(i=0;i<IMAX;i++){
        fprintf(fp,"%lf,%d,%lf\n",time,i,T[i]);      /* 計算結果を保存 */
    }

    fclose(fp);
    printf("FINISH!\n");
}
_images/netudendou_res_ex.png

　\(\Delta\)tを10\(^{-6}\)にとったときの計算結果を示します。このような結果になればプログラミングが成功したと思います．

　この方法は大変分かりやすくてよいのですが、欠点があることが知られています．プログラミングではtとxのそれぞれのステップ幅を\(\Delta\)t、\(\Delta\)xとして、自由に値を決められますが、これが\(\Delta\)t/\(\Delta\)x\(^2\) \(>\) 0.5を満たすようになると値がおかしくなるというものです．プログラミングが終わった人はこれについても試し、原因について考えてみましょう。

陰解法の説明

　陽解法では直感的に分かりやすい反面、計算結果が正しくない場合があるということでしたが、ここではもう一つ別の解法について考えてみましょう。この方法は陰解法といい、収束性が良いのですが、プロブラミングは少し複雑になります．

　普通に差分化した場合、xについての２回微分の右辺の時刻は、左辺のtと同じにしたと思います．ここでは、t+\(\Delta\)tにしてみましょう。（プログラミング上はt=n+1にすることに対応します．）このとき式は下のようになります．

_images/implicit.png

　この式を詳しく見ると、時刻n+1、場所iのポイントの温度（下の表の赤いセル）を求めるには、時刻nのiの他に、時刻n+1の場所i-1、i、i+1の三つの場所の温度が必要になることが分かります．つまり現在の温度を求めるのに、現在の隣の温度を知らなければならないことが分かります。これでは先ほどのような考え方では解が求まりません。ここでは各時刻の温度分布にを求めるために、各時刻のステップについて連立方程式を解くことが必要になります。

_images/netudendou_implicit.png

　それでは解くべき連立方程式を書き出してみます．時刻t=nの温度は既に分かっているものとして変数から外し、ここでは\(D_j\)とおきます。また\(T_{j-1}\)、\(T_j\)、\(T_{j+1}\)の係数をそれぞれ\(A_j\)、\(B_j\)、\(C_j\)とおけば下の式になります．

_images/implicit2.png

　さらに行列表示にすると下のようになります。

_images/implicit2_matrix.png

　この行列の逆行列を求めることができれば、温度分布が求まるということが分かります。

解法例

_images/implicit3.png

と仮定して，差分化した先の式に代入し、\(T^{*}_{j-1}\)を消去する。以下の式が得られる。

_images/implicit4.png

式を見比べると、

_images/implicit5.png

であることがわかる。ここで境界条件として、

_images/boundary.png

を設定し、以下の式について考える。

_images/implicit6.png

j=0のとき、

_images/implicit6_0.png

この式がいつも成り立つには

_images/boundary2.png

となる。

陰解法ソースファイル

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>

#define IMAX    101                                    /* x方向の計算領域 */
#define TMAX    50                                    /* 計算終了時間 */
#define dt  0.01                                       /* 時間刻み */
#define dx  0.01                                       /* x方向の格子間隔 */

void main(void){
    FILE *fp;                                          /*ファイルのため*/
    int i,t;
    double time,T[IMAX],Tnew[IMAX];                    /* T[100]  */
    double A[IMAX],B[IMAX],C[IMAX],D[IMAX];            /* TDMAの変数 */
    double P[IMAX],Q[IMAX];                            /* 逆行列用  */

    if((fp=fopen("result.dat","w"))==NULL){
        printf("Can not creat file.\n");
        exit(-1);
    }

    /* 境界条件 */
    T[0]=1;Tnew[0]=1;T[IMAX-1]=0;Tnew[IMAX-1]=0;
    P[0]=1;
    Q[0]=2;

    /* 初期条件*/
    for(i=1;i<IMAX-1;i++){
        T[i]=0;
    }

    /* TDMA */
    for(i=1;i<IMAX-1;i++){
        A[i]=-dt/dx/dx;
        B[i]=(1+2*dt/dx/dx);
        C[i]=-dt/dx/dx;
        D[i]=0;
    }

    time=0;

    for(t=0;t<TMAX;t++){
    time=time+dt;

        for(i=1;i<IMAX-1;i++){                                  /*TDMAを解く*/
            P[i]=C[i]/(B[i]-A[i]*P[i-1]);
            Q[i]=(D[i]-A[i]*Q[i-1])/(B[i]-A[i]*P[i-1]);
        }

        for(i=IMAX-1;i>1;i--){                                  /*温度を得る*/
            Tnew[i-1]=Q[i-1]-P[i-1]*Tnew[i];
        }

        for(i=0;i<IMAX;i++){                                    /*値の更新*/
            //fprintf(fp,"%lf,%d,%lf\n",time,i,Tnew[i]);
            T[i]=Tnew[i];
            D[i]=Tnew[i];
        }
    }

    for(i=0;i<IMAX;i++){
        // printf("%d,%lf\n",i,T[i]);
        fprintf(fp,"%lf,%d,%lf\n",time,i,Tnew[i]);
    }
    fclose(fp);
    printf("FINISH!\n");
}
_images/netudendou_res_imp.png

\(\Delta\)tを10\(^{-2}\)にとったときの計算結果を示します。このような結果になればプログラミングが成功したと思います．

\(\Delta\)t/\(\Delta\)x\(^2\) \(>\) 0.5を満たす場合でも計算がちゃんと終了することを確かめてみましょう。

分子動力学計算

分子動力学とは

_images/gainen.png

　分子（原子）を質点とし，運動をニュートンの運動方程式で規定する．分子 i の質量をm，位置ベクトル\({\mathbf r}_i\), 分子 i に作用する力を\({\mathbf f}_i\)としたときの運動方程式は、

\[m_{i} \frac{d^2 {\mathbf r}_i}{d t^2} = {\mathbf f}_{i}\]

　ある領域内にN個の分子が存在すれば，N個の運動方程式が存在し(i=1,...,N)，\({\mathbf f}_i\)を介して互いの運動に影響を及ぼす。

　微小時間 \(\Delta\)t ごとの粒子（分子）の位置\({\mathbf r}_i\) と速度 \({\mathbf v}_i\) を求める。

差分

　時刻 t から微小時間\(\Delta\)t 後の粒子の位置\({\mathbf r}_i (t + \Delta t)\) はテイラー展開より、

\[{\mathbf r}_{i} (t + \Delta t) = {\mathbf r}_{i} (t) + \frac{d {\mathbf r}_i}{d t} \Delta t + \frac{1}{2!}\frac{d^2 {\mathbf r}_i}{d t^2} (\Delta t)^2 + \frac{1}{3!}\frac{d^3 {\mathbf r}_i}{d t^3} (\Delta t)^3 + \cdot\cdot\cdot\]

　時刻 t から微小時間\(\Delta\)t 前の粒子の位置\({\mathbf r}_i (t - \Delta t)\)はテイラー展開より、

\[{\mathbf r}_{i} (t - \Delta t) = {\mathbf r}_{i} (t) - \frac{d {\mathbf r}_i}{d t} \Delta t + \frac{1}{2!}\frac{d^2 {\mathbf r}_i}{d t^2} (\Delta t)^2 - \frac{1}{3!}\frac{d^3 {\mathbf r}_i}{d t^3} (\Delta t)^3 + \cdot\cdot\cdot\]

　２つの式の両辺を加えて1階の微分項を消去する。

\[{\mathbf r}_{i} (t + \Delta t) = 2{\mathbf r}_{i} (t) - {\mathbf r}_{i} (t - \Delta t) + \frac{1}{2!}\frac{d^2 {\mathbf r}_i}{d t^2} (\Delta t)^2 + O((\Delta t)^4)\]

　上の式にニュートンの運動方程式を代入すると、

\[{\mathbf r}_{i} (t + \Delta t) = 2{\mathbf r}_{i} (t) - {\mathbf r}_{i} (t - \Delta t) + \frac{f_i}{m_i}(\Delta t)^2\]

　または、

\[{\mathbf r}_{i} (t + \Delta t) = {\mathbf r}_{i} (t) - v_i (t)\Delta t + \frac{f_i}{m_i}(\Delta t)^2\]

が得られる。

位置が条件として与えられる場合

_images/kaihou1.png

　時刻 t - Δt での位置および時刻tにおける位置と力から時刻 t + Δt における位置がわかる。このとき速度は図中の式で与えられる。

位置と速度が条件として与えられる場合

_images/kaihou2.png

時刻 t における位置と力および速度が与えられるとき、速度を図中の式として式(5)に代入し、式(6)をつかう．

力とポテンシャル

　粒子 i に働く力\({\mathbf f}_{i}\)は、

\[{\mathbf f}_{i} (t) = \sum_{j=1, j\neq i}^N \frac{\partial \Phi (r_{ij})}{\partial r_{ij}} \frac{{\mathbf r}_{ij}}{r_{ij}}\]

　ここで、\(\Phi\)は2体間ポテンシャル（ポテンシャルが2個の原子間の距離のみを独立変数とする），r\(_{ij}\)は粒子間の距離を表す。この2体間ポテンシャルは経験的、つまり量子力学の厳密な理論ではなく，ポテンシャルを微分可能な関数形で仮定し，実験的事実に合致するように係数を決めたものが利用される。例えば以下のようなポテンシャルが提案されている。

Lennard-Jonesポテンシャル

　Van-der-Waals力により結合されている不活性ガスに対して提案された経験的ポテンシャル。

_images/LJ.png _images/LJ_eq.png

Morseポテンシャル

　2原子分子の相互作用を記述するものとして提案された経験的ポテンシャル。

_images/morse.png _images/morse_eq.png

課題：3粒子モデルの変形問題

　以下に示す３粒子モデルについて、粒子間にMorseポテンシャルを仮定し，粒子２の運動を数値計算により求めよ。

_images/report.png
Morseポテンシャルのパラメータ
Cu

  • D = 0.3429 [eV]

  • A = 2.866 [\({\rm A} ^{-1}\)]

  • r\(_0\) = 2.886 [\({\rm A}\)]

  • v\(_0\) = 1.00*10\(^{-4}\) [\({\rm A}\) /fs]

  • r(0) =1.05 r\(_0\)

Ca

  • D = 0.1623 [eV]

  • A = 1.5079 [\({\rm A} ^{-1}\)]

  • r\(_0\) = 4.569 [\({\rm A}\)]

  • v\(_0\) = 1.00*10\(^{-4}\) [\({\rm A}\) /fs]

  • r(0) =0.9r\(_0\)

分子動力学ソースファイル

/*  Satoshi Iikubo    Apr. 21, 2011 */

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>

void main(void){
    /*使用する変数を定義*/
    double f;
    double x12;
    double r6,r12;
    double x1_now;
    double x2_now;
    double dt;
    double vx1;
    double vx2;
    double x1_next;
    double x2_next;
    FILE *fp1;
    int t;

    /*計算に使う初期値を指定*/
    t=0;//時間
    dt=;//時間ステップ
    x1_now=0;//原子1の位置
    x2_now=;//原子２の位置
    vx1=;//原子１の速度
    vx2=;//原子２の速度
    //x1_nexti=0;
    //x2_nexti=2;

    /*データ出力に使うファイルの設定*/
    if((fp1=fopen("result.dat","w"))==NULL){
        printf("can not open file\n");
        exit(-1);
    }

    /*ファイルに出力*/
    fprintf(fp1,"%d,%lf,%lf,%lf\n",t,x1_now,x2_now,r2);

    /*シミュレーションの開始*/
    for(t=1;t<5000;t++){
        /*ポテンシャルを使って力の計算*/
        x12=;
        r6=;
        r12=;
        f=;

        /*Newton 方程式に代入*/
        x1_next=;
        x2_next=;

        /*ファイルに出力*/
        fprintf(fp1,"%d,%lf,%lf\n",t,x1_next,x2_next);

        /*次のサイクルのために位置と速さを更新*/
        vx1=;
        vx2=;
        x1_now=;
        x2_now=;

    }
    fclose(fp1);
}