Pythonで数値計算①：Population growth

学生のころからときおり数値計算プログラムを作っていて、もっとわかりやすい教材ないのかなーといつも思っていました。

数値計算のプログラム作成って、プログラミング言語の基本構文や関数を知るだけでは作れないんですよね。

細かいノウハウや式変形、論理思考が必要です。
こういったものは実際に人の書いたコードを見ながら、習得していくのが一番手っ取り早いと思います。

このシリーズでは、簡単な数値計算の問題を解いてみて、その解放やノウハウを皆さんに紹介したいと思います。

このノートの対象者

・Pythonで科学計算をしたい（することになってしまった）初心者の人
　　例えば、理工系の大学生や理系科目が好きな人
・Pythonの基本文法がわかる人
・Pythonで外部モジュールがインストールできる人
　　pipかcondaどちらかできればOK
・微分や偏微分方程式を見たことがある人
　　少しわかればOK

かなりかみ砕いて説明するので、本職の人からしたら物足りないかも…。

【問題】Population growth

今回は、Population growth（人口増加）の簡単なモデル、
マルサスモデルを数値的に解いてみましょう。

人口の増加速度（dP/dt）が、人口に比例するという式ですね。
簡単に言えば、人口が少ないうちはあまり増加しませんが、
時間が経って人口が増加すると、どんどん増加速度も大きくなるというモデルです。

解析解

この微分方程式は、簡単なので解析解があります。
解析解というのは式を変形したり積分、微分したりして、導出した解です。

こちらが解です。疑っている方は、元の式に代入してみてください。
右辺＝左辺になります！

時間 vs 人口をグラフにするとこんな感じです。

ちなみに、この図はpythonで出力しています。

import matplotlib.pyplot as plt
import numpy as np

# １．解析解を定義しています。
def analytical_solution(t, k, p0):
   return p0 * np.exp(k * t)

# ２．t=0から50までの配列作成
t = np.linspace(0, 5)

# ３. matplotlibでプロット
plt.plot(t, analytical_solution(t, 1, 10), label="k=1.0, p0=10")

plt.title("population growth")
plt.xlabel("time")
plt.ylabel("population")
plt.legend()
plt.show()

１．解析解の関数を定義

２．numpyのlinspaceでt=0から50までの配列を作成
　　⇒NumPy linspace | 等間隔の数列を作成する！

３．matplotlibで簡単にプロット
　　⇒matplotlib | pyplotで基本的なグラフを簡単に作れるようになるには？

今回のような簡単な微分方程式では解析解を導出することができることがあります。
解析解を用いて、方程式の答えを出すことを「解析的に解く」とか言ったりします。

一方で、複雑な微分方程式は解析的に解けないことが多いです。
というか、解析的に解ける微分方程式では、複雑な現象を記述できないことが多いです。

そういう時は、微分方程式を「数値的に解く」ことになります。

数値解

今回は差分法という手法で、微分方程式を数値的に解いてみましょう。
差分法は直感的で、入門にぴったりです。

私は特に差分法との付き合いが長く、学生時代は差分法で１万行近いコードの数値シミュレータを作成していました（長けりゃいいってわけじゃないですよね・・・）。

差分法では、微分を次のように近似します。

微分dP/dtは、時間当たりのPの変化量と言えます。

今、変数の値がP(n)だったものが、Δt秒後にP(n+1)になっていたとしたら、
・変化量＝P(n+1)ーP(n)
・経過時間=Δt
なので上の式で近似できそうですよね。

もう少しだけ正確に言うと、微分項をテーラー展開して二次以降の項を無視しているんですが、そういう細かいことは、後から勉強しても良いと思います（持論を展開）

時間の経過をわかりやすいように時間を表すインデックスnを変数Pの上に書くことが多いです。これはn乗ではないので注意しましょう（慣例で上につけます）。
・t=0のときはn=0
・t=Δtのときは、n=1
・t=2Δtのときは、n=2…
としましょう。Δt秒経つごとに、nが1増えます。

では、この近似式を使って初めの微分方程式を解いてみましょう。

もう普通の漸化式ですね。n=0のときはP0ですから、nを1ずつ増やせばとけてしまいます。

では、実際にpythonで解いてみましょう。

１．計算を行う時刻tを配列で用意
　　同時にタイムステップ幅Δtを用意
　　⇒NumPy linspace | 等間隔の数列を作成する！
※Δtを十分小さくしましょう！近似式はΔtが小さい時に精度が上がります。

２．各時刻のPを格納できるように配列Pを用意
　　（今回はグラフを描くために全時刻でのPの値を保存します）
　　⇒numpy zeros, ones | 初期化した配列を生成

３．n=1から順に漸化式を解きます。

４．プロットします。
　　⇒matplotlib | pyplotで基本的なグラフを簡単に作れるようになるには？

さてできましたか？こんな感じのグラフが出力できればOKです！

私のコードを紹介します。実行するだけで、上のようなグラフが作れます。

解析解のコードを実行できていれば、こちらのコードも動くはずです。