こちらに PDF が置いてあります.
サイトでは見づらい方はこちらをご覧ください.
近似とは?
今回アニメーションはありませんが引き続きシミュレーションの話をします.
今回のテーマは近似についてです.
シミュレーションというのはいいし,
近似というのもいいけど,
実際どのくらいよく近似できてるの?
こういう問題を考えてみましょう.
ちゃんと使えばすごく精度はいいですよ,
というのが今回の主張です.
- 微分方程式というのはいいけど近似なんて雑なことしたくない!
- 近似なんてやって意味あるの? どうせ大雑把で使いものにならないんじゃないの?
- 百歩譲って微分方程式は役に立つとしてもその近似計算は本当に役に立つの?
こんな疑問に答えるのが今回のテーマです.
放射性物質の崩壊に関する微分方程式
で, ちょっと不吉な例ではあります.
しかしというか残念ながらというか,
役に立つ例になってしまったので放射性物質の崩壊に関する微分方程式を考えてみましょう.
いちおう微分方程式から書いておきますね.
\begin{align}\label{math_refuge00002}
\frac{dx}{dt}
=
– c x.
\end{align}
気分だけ説明しておくと放射性物質が単位時間あたりに崩壊するスピードはそのときの放射性物質の量に比例する,
というのを式で書いた形です.
$c$ には物理的に大事な意味がありますが式や計算を軽くするためにここでは $c = 1$ にして話を進めます.
そして前回のようにこの微分方程式を近似した式を書いてみます.
\begin{align}\label{math_refuge00003}
\frac{x_{n+1} – x_{n}}{h}
=
– x_{n}.
\end{align}
これは時間間隔 $h$ が十分小さければ時間 $h$ だけ離れた放射性物質の量 $x_{n+1}$ と $x_{n}$ の関係が上の式で書けることを意味しています.
近似計算はどのくらい正確なの?
物理的な話はこのくらいにして,
この式にしたがって近似計算 (シミュレーション) したとき,
本当の答えと近似した解がどのくらい近いのかを考えてみます.
何でこれにしたかというともとの方程式の答えが厳密に書けるからです.
あなたが指数関数, 特に自然対数の底 $e$ をご存知なら,
式 \eqref{math_refuge00002} の答えは $x(t) = C e^{- t}$ と書けます.
指数関数の微分をご存知ならこれを \eqref{math_refuge00002} に代入してみて本当に方程式の答えになっていることを確認してみてください.
この答えと式 \eqref{math_refuge00003} にしたがって計算した答えがどのくらい一致するか,
コンピュータに計算させてみましょう.
結果はこの記事の後半,
プログラムパートに載せてあります.
近似なんていい加減ことしたくない!
あなたは「近似なんて雑なことしていいの?」なんて思っているかもしれません.
でも記事の後半にある「比較のために重ねる」節のグラフを見るとわかるように,
厳密な様子ともよく合っているのでちゃんと使えば問題ありません.
「ちゃんと使う」の「ちゃんと」が難しいと言われればそれはもちろんそうなんですが,
それは本当に難しい話です.
せっかくなので分割をどんどん大きくしていって近似の精度が上がっていく様子もグラフにしておきました.
プログラムが書けるとこういうのもささっとやれます.
遊びやすくなるのは本当にありがたいですね.
中高生の頃にこれを知っていたら数学や物理だけじゃなくてプログラミングにもはまっていたと思います.
やってみるとわかるんですが,
書いたプログラムで意図通りの結果が出るとけっこう嬉しいです.
パラメータをいじったり方程式を変えたり,
ちょろっといじって遊べる要素も増えて遊べる範囲が広がるのもいいですね.
Python プログラミングに関する資料,
数値計算・シミュレーションに関する資料は GitHub に上げてありますし,
今後も講座の進展に合わせて少しずつ増やしていくのでぜひあなたも遊んで遊んでみてください.
このサイトにもいろいろな情報がありますよ.
数学的なポイントまとめ
あと今回の中高数学上のポイントを復習しておきます.
指数関数とか自然対数の底とか,
はたまたその微分が出てきました.
放射性物質の崩壊をちゃんと調べるにはこういう数学が必要なんです.
次回は経済や生物ネタです!
次回はまた別の微分方程式を紹介します.
経済学や生物学で出てくる微分方程式ですよ.
これが終わったらもう少し数学的にちゃんとした話をしましょう.
まずは微分方程式の射程距離の長さを知ってほしいからです.
もっとちゃんと数学の説明してほしいというあなた,
もうちょっと我慢してください.
今回もアンケートがあるのでぜひ回答をお願いします.
ではまた次回をお楽しみに!
プログラミングパート
放射性物質の崩壊
一番単純でしかも実際に使われる微分方程式として,
まずは 1 階の線型常微分方程式を考えよう.
ちょっと不吉な例であるが放射性物質の崩壊の方程式を紹介する.
導出をしたければちゃんと物理を勉強してもらう必要がある.
ここでは物理は省略して数学に集中する.
\begin{align}
\frac{dx}{dt}
=
– c u.
\end{align}
厳密解は $x(t) = C_0 e^{-ct}$ だ.
初期値を設定すれば $C_0$ はそこから決まる.
微分を単純に離散化すると次のようになる.
\begin{align}
\frac{x_{n+1} – x_{n}}{\Delta t}
=
-c x_{n}.
\end{align}
時間っぽく見えるよう $\Delta t$ と書いてみた.
タイプが面倒なので $h$ と書くこともある.
整理すると次の通り.
\begin{align}
x_{n+1}
=
x_{n} – c (\Delta t) x_{n}.
\end{align}
これに沿って計算したのがいわゆるオイラー法.
コードに落としていこう.
厳密解のグラフ
まずは厳密解 $x(t) = e^{-t}$ をグラフに描こう.
初期値は 1, 上の定数も $c = 1$ としている.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
c = 1
init = 1
nt = 100
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['exact'])
plt.show()
RESULT
近次解
今度は近次解を計算してグラフにしてみよう.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def radioactive_euler(nt, init = 10):
dt = 2 / (nt - 1)
# 初期条件設定
x = np.zeros(nt)
x[0] = init
for i in range(1, nt):
x[i] = x[i-1] - c * dt * x[i-1]
# ベクトル計算で書き直したい
return x
c = 1
init = 1
nt = 101
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
plt.legend(['approximation'])
plt.show()
RESULT
比較のために重ねる
これだけではよくわからないので重ねて描いてみる.
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
ほぼピッタリ重なって見える.
シミュレーションの精度はけっこう良さそうだ.
もちろん細かいことを言えばいろいろあるけれども,
それはもっと進んだお話なので一旦不問にしておこう.
細かく刻んだ方がいい近似になるはずだ
上のグラフは区間を 100 分割して計算した.
nt = 101 と言うところで分割を指定している.
101 と言うのは分点数の指定なので実際の分割としては 100 等分なのだ.
ここで刻みが荒いといい近似にならないことを大雑把に調べておこう.
2 等分
nt = 3
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
4 等分
何となくそれっぽい形にはなっている.
もちろん精度は大したことない.
nt = 5
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULTS
10 等分
あまり適当なことを言うのも良くないが, けっこう滑らかに見える.
精度はまだまだ.
nt = 11
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
50 等分
まだ差が目で見える.
nt = 51
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
1000 等分
100 は最初にやってあるから一気に大きくしてみた.
nt = 1001
x_approx = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)
nt = 101
t = np.linspace(0, 2, nt)
x_exact = init * np.exp(- c * t)
plt.plot(t, x_exact)
plt.legend(['approximation', 'exact'])
plt.show()
RESULT
参考
コンピュータで計算するとき,
いろいろな都合があって分割数を多くすれば単純に精度が上がっていくと言うわけではない.
桁落ちや丸め誤差のような問題がある.
それはそれで別途検討する必要があるのだ.