【科学技術計算講座1-5】Pythonでオイラー法をプログラミングする

前回、数値計算法としてオイラー法を説明しました。

いよいよ今回から、Pythonでプログラミングして、問題を解いていきます。

問題のおさらい

今回解くべき問題をおさらいしておきます。

問題：　Aという物質が時間がたつとBという物質に変化するとき、AとBの物質の濃度の時間変化を求めたい。

問題は一分子反応の化学反応で、AとBの濃度変化を求めるというものです。そして、その方程式は、

$$\left\{ \begin{aligned} \frac{d [\rm{A}]}{dt} &= -k [\rm{A}]\\ \frac{d [\rm{B}]}{dt} &= k [\rm{A}] \end{aligned} \right. \tag{1}$$

という微分方程式でした。さらに、この微分方程式は、

$$f_{A,n} = -k [{\rm{A}}]_n \tag{2}$$

$$f_{B,n}= k [{\rm{A}}]_n \tag{3}$$

$$[{\rm{A}}]_{n+1} = [{\rm{A}}]_n + \Delta t \ f_{A,n} \tag{4}$$

$$[{\rm{B}}]_{n+1} = [{\rm{B}}]_n + \Delta t f_{B,n} \tag{5}$$

$$t_{n+1} = t_n + \Delta t \tag{6}$$

という形でオイラー法という数値計算法で解くことができるんでしたね。そして、初期値を与えてやれば、時々刻々芋づる式に計算していくことができるということでした。

プログラミングの前に、問題の条件を決めておきましょう。まず、初期（時刻０秒）に、Aの濃度は 1.0 [mol/L]、Bの濃度は 0.0 [mol/L] あったとしましょう。そして、100秒までA、Bがどう変化していくか計算してみます。反応速度定数 $k$ は、0.1 [1/s]（１秒あたり10%変化する）とします。

Pythonでプログラミングしてみる

ではPythonで、この問題をプログラミングしてみましょう。

初期値を入力する

まず、テキストエディタ（メモ帳でも何でもよいです）を開いてください。はじめに、初期値を与えてやるので、次のように書いてください。

a = 1.0
b = 0.0
t = 0.0

まず、Aの濃度を a という変数で表します。Aの初期の濃度は 1.0 なので、a に 1.0 を代入します。これが、a = 1.0 のところです。

次に、Bの濃度を b という変数で表します。Bの初期濃度は 0.0 なので、b = 0.0 です。

そして、時刻は t という変数にします。初期は時刻０なので、 t = 0.0 とします。

他のパラメータを入力する

続いて、反応速度定数 $k$ を与えてやります。この変数は、k としましょう。これは 0.1 なので次のように書きます。

k = 0.1

そしてあと一つ、時間刻み $\Delta t$ を与える必要があります。時間刻みは前回説明したように、大きいと粗い計算になってしまいます。今はとりあえず0.1秒としておきます。100秒間の計算をするので、1000回ステップ計算をすることになります。これだと十分滑らかそうです。時間刻みは、dt という変数にします。なので、次のように書きます。

dt = 0.1

ステップ計算部分を書く

ここまでで、オイラー法の計算に必要な条件が設定できました。次にオイラー法の時々刻々ステップ式に計算する部分を書いていきます。

for文による繰り返し処理

オイラー法では、(2)～(6)式を初期値から始めて、次々に繰り返し計算して行けばいいんでした。今回は1000回繰り返すのですが、まさか1000回これらを書いていくということはしません。大変なので。。。

プログラミング言語には、このような繰り返しをさせるための命令が用意されています。Pythonでは for文というのが用意されています。C言語とか多くの言語でもfor文というのがあります。

for文は指定回数、処理を繰り返すという命令です。例えば1000回処理をくりかえしたい場合は、

for i in range(1000):
    print(i)

のように書きます。range()の回数だけfor文の下の処理を繰り返しなさいというものです。ここで i は変数です。i は繰り返しのたびに 1 ずつ増えていきます。この例では、print()で i を表示するという処理を1000回繰り返すことになります。

ちなみに、print(i)は字下げされていますね。Pythonでは、for文の下の字下げされている部分を繰り返すという決まりになっています。字下げは４つのスペースを入れることが多いようですが、２つとかでも構いません。ただ、複数行を繰り返す場合は、字下げの数は揃えておいてください。

また、for文の最後の : （コロン）を忘れないようにつけてください。

それでは、途中ですが、ここまでを保存してプログラムを実行してみましょう。ファイル名は何でもよいですが、拡張子を.pyとしておいてください。とりあえず、test.py としておきます。

Anacondaのインストールで説明したように、Anaconda Promptのショートカットを test.py と同じフォルダーにコピーしてきます。（Anaconda Promptのプロパティの設定を忘れないようにしておいてください。）

Anaconda Promptを起動して、pythonコマンドでファイルを実行します。

python test.py

実行結果は、

のように、0 から 999 まで数字が出力されたと思います。これが、for文で出力された i の値です。i が1000回繰り返し出力されているのがわかります。i は０から始まり１ずつ増えていきます。

ちなみに、このようなfor文による繰り返しを forループ と言ったりします。

オイラー法の繰り返し処理

では、引き続きオイラー法のプログラムに戻りましょう。先程のfor文のところを、次のように書きます。

for i in range(1000):
    fa = -k * a
    fb = k * a
    
    a = a + dt * fa
    b = b + dt * fb
    t = t + dt

    print(i, t, a, b)

(2)式の関数 $f_A$ が、初めの fa に相当します。k と a は前にに入力した数値が入ります。k と a は掛け算ですが、これは * という記号です。マイナスは、- です。これは式そのままですね。fa = -k * a であれば、-k * a の計算結果が fa に代入されるという意味になります。同様に(3)式の関数 $f_B$ を fb として計算しています。

次に(4)式ですが、これは a = a + dt * fa の所です。まず、右辺ですがこれは(4)式の右辺そのままです。+ は足し算ですね。a と dt は最初に入力した値が入ります。fa は直前で計算した結果が入ります。

ここまではいいと思いますが、問題は左辺です。ここでは、a + dt * fa をまた a に代入しています。(4)式では、今のステップを $n$、次のステップを $n+1$ としています。これだと、ステップごとに別の変数を用意しなくてはいけないので大変です。プログラミングではこういう場合、同じ変数に代入していくという書き方をよくします。

例えば次のようなプログラムを例にとります。

n = 1
n = n + 2

このプログラムでは、１行目で変数 n に１を入力しています。なので n は１です。２行目では、右辺で n + 2 としているので、計算すると３です。これを左辺の変数 n に代入することになります。これで、n は３になりました。１行目ではn は１でしたが、２行目で n は３になるのです。普通の数学ではこんな書き方はしませんが、プログラミングではごく普通の書き方です。慣れない方は、最初とまどうかもしれませんが、理解しておいてください。

さてオイラー法プログラムに戻ります。(5)式、(6)式も全く同じように、同じ変数に代入されています。これで、a、b、t は今の時刻の値を使って、次の時刻の a、b、t を求めたことになります。

最後に、print()で、i（ステップ数）、t（時刻）、a（Aの濃度）、b（Bの濃度）を出力します。

これをfor文で繰り返します。すると次のステップでは、前に計算された a、b、tを使って、また次のステップを計算することになります。これをfor文の指定回数くりかえします。これで、時々刻々、値が求まっていくのです。

ここまでで、一分子反応をオイラー法で解くプログラムが出来ました。ここまでを保存しましょう。名前を euler.py としておきます。

念の為、ここまでのプログラムを載せておきます。

プログラムを表示

# parameter
a = 1.0
b = 0.0
t = 0.0
k = 0.1
dt = 0.1

# time loop
for i in range(1000):
    fa = -k * a
    fb = k * a
    
    a = a + dt * fa
    b = b + dt * fb
    t = t + dt

    print(i, t, a, b)

では、Anaconda Promptで、これを実行してましょう。

python euler.py

実行結果は、こんな感じです。

0 0.1 0.99 0.010000000000000002
1 0.2 0.9801 0.0199
2 0.30000000000000004 0.970299 0.029701
3 0.4 0.9605960100000001 0.03940399
4 0.5 0.9509900499 0.049009950100000005
5 0.6 0.941480149401 0.05851985059900001
6 0.7 0.93206534790699 0.06793465209301
7 0.7999999999999999 0.9227446944279202 0.07725530557207991
8 0.8999999999999999 0.9135172474836409 0.0864827525163591
...
997 99.7999999999986 4.4047798602855e-05 0.9999559522013993
998 99.8999999999986 4.3607320616826455e-05 0.9999563926793853
999 99.9999999999986 4.317124741065819e-05 0.9999568287525915

1000回繰り返し出力されてますね。これで、各時刻でのAとBの濃度が計算できました。（コンピュータ環境によって微妙に数値が違っているかもしれませんが、だいたい合ってればOKです。大きく違う場合は、どこかプログラムが間違っているかもしれませんので、もう一度見直してください。）

数値誤差

しかし、、、とっても見にくいですね。それに、２列目は時刻 t を出力したのですが、0.30000000000000004 とか 0.7999999999999999 とかおかしな数値になっています。時間刻みは0.1秒だったので、0.1 ずつきれいに出力されていてほしいものです。

これは数値計算には付き物の数値誤差というやつです。コンピュータの中では２進数で計算が行われます。このとき小数点はちょっとやっかいです。10進数で0.1 という数値を２進数に直すと、0.000110011001100... という 1100 が無限に続く循環小数というものになります。無限に続くものは扱えないので、どこか適当なところで打ち切って計算されています。つまり、0.3 が0.30000000000000004 とか、0.8 が 0.7999999999999999 とかの結果になってしまうのです。これは丸め誤差というものですが、数値誤差には他にも様々な種類があります。ここでは、詳しくは説明しませんが、数値計算を行うとこのような誤差を必ず含んでいると知っておいてください。

出力のフォーマット

このままでは見にくいですし、小さな誤差は切り良い数値にしたいと思います。そこで、最後のprint(i, t, a, b)のところを次のように書き換えてみてください。

print('i: {0:4d}, t: {1:6.2f}, A: {2:9.6f}, B: {3:9.6f}'.format(i, t, a, b))

以前、Hello world!を出力したときのように、' ' で囲まれている部分は文字として出力されます。そのとき {0:4d} のように {} で囲まれた部分に、後ろの .format() で指定されている変数が表示されます。

{0:4d} の 0 のところは番号になっていて、formatの中で指定されている順番（0番から始まる）になっています。つまり、図のように i、t、a、b はそれぞれ０、１、２、３の番号の {} に対応しています。{} の番号の次の：以下が数値をどのような形式で表示するかの指定です。

4d は、整数を４桁で表示するという意味です。dは整数を表します。

6.2f のfは浮動小数点（いわゆる小数値）を表し、全部で６桁で小数点以下２桁まで表示します。値は指定した小数点の桁になるよう四捨五入して表示されます。

9.6f も同様に、浮動小数点を全部で９桁、小数点以下を６桁まで表示します。

これでもう一度 Python で実行してみてください。

プログラムを表示

# parameter
a = 1.0
b = 0.0
t = 0.0
k = 0.1
dt = 0.1

# time loop
for i in range(1000):
    fa = -k * a
    fb = k * a
    
    a = a + dt * fa
    b = b + dt * fb
    t = t + dt

    print('i: {0:4d}, t: {1:6.2f}, A: {2:9.6f}, B: {3:9.6f}'.format(i, t, a, b))

計算結果は次のようになります。

i:    0, t:   0.10, A:  0.990000, B:  0.010000
i:    1, t:   0.20, A:  0.980100, B:  0.019900
i:    2, t:   0.30, A:  0.970299, B:  0.029701
i:    3, t:   0.40, A:  0.960596, B:  0.039404
...
i:  997, t:  99.80, A:  0.000044, B:  0.999956
i:  998, t:  99.90, A:  0.000044, B:  0.999956
i:  999, t: 100.00, A:  0.000043, B:  0.999957

今回は見やすくスッキリしました。どの数値がどの変数に相当するかも表示されていますし、微小な数値誤差も切りの良い数値にされています。

これで、AとBの濃度がどのように変化していくのかわかりました。Aの濃度は1.0[mol/L] から徐々に小さくなっていき、7秒後には半分の 0.5[mol/L] より小さくなって、46秒後には1%以下の 0.01[mol/L] より小さくなっていることがわかります。Bの濃度は逆に 0.0[mol/L] から徐々に増えて、同じく７秒後には 0.5[mol/L] になり、46秒後には 0.99[mol/L] を超えています。

このように、オイラー法を使った数値計算でAとBの濃度がどのように変化するのか計算することができました。数値データだけが欲しい場合は、これで十分です。しかし、このデータをグラフにして、視覚的に見ることができればもっと有益です。データが表示されているので、これをコピペしてExcelにでも貼り付けてグラフにしてもいいのですが、いちいち面倒です。せっかくなので、Pythonでやってしまいたいものです。次回は、Pythonでグラフを作成してみたいと思います。

まとめ

今回は、Pythonでプログラミングして、実際に計算してみました。次回は、グラフ表示して、視覚的に確認するところを実施していきます。

補足

方程式の数

数値データをみて賢明な方はお気づきかと思います。今回考えている問題は、Aが減った分だけBが増えるというものです。これは、ある時刻のBの濃度は、Aの初期値からその時刻のAの濃度を引いてやれば簡単に求めることができるのです。AとB両方とも微分方程式をたてて計算していますが、Aのみ計算すれば単純な引き算でBの濃度は求まります。今回は、連立方程式の練習も兼ねてあえて両方とも計算しています。

累算代入演算子

同じ変数に代入していく、

n = n + 2

のような書き方をしましたが、これは次のように書くこともできます。

n += 2

+= というのは累算代入演算子とよばれ、左辺の変数に右辺を足して、再び左辺の変数に代入するというものです。これを使うと

a += dt * fa

のように書くことができます。この書き方はよく使いますが、どちらを使ってもかまいません。

コメント行

サンプルのプログラムに、# で始まる行が書かれています。これはコメント行といい、この行は実行されません。何かコメントを残しておきたいときに使います。また、一度プログラムで書いた行を # をつけて実行しないようにすることをコメントアウトといったりします。

また、何も書いていない空行をプログラム中に入れることもあります。これはプログラムを見やすくするために入れます。

←前回　数式を計算する方法（オイラー法）

→次回　Pythonでグラフを表示する

全体の目次