このサイトは学部では早稲田で物理を, 修士では東大で数学を専攻し, 今も非アカデミックの立場で数学や物理と向き合っている一市民の奮闘の記録です. 運営者情報および運営理念についてはこちらをご覧ください.

理系のための総合語学・リベラルアーツの視点から数学・物理・プログラミング・語学 (特に英語) の情報を発信しています. コンテンツアーカイブに見やすくまとめているのでぜひご覧ください.

こちらに 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 に上げてありますし,
今後も講座の進展に合わせて少しずつ増やしていくのでぜひあなたも遊んで遊んでみてください.
このサイトにもいろいろな情報がありますよ.

• https://github.com/phasetr/mathcodes

数学的なポイントまとめ

あと今回の中高数学上のポイントを復習しておきます.
指数関数とか自然対数の底とか,
はたまたその微分が出てきました.
放射性物質の崩壊をちゃんと調べるにはこういう数学が必要なんです.

次回は経済や生物ネタです!

次回はまた別の微分方程式を紹介します.
経済学や生物学で出てくる微分方程式ですよ.
これが終わったらもう少し数学的にちゃんとした話をしましょう.
まずは微分方程式の射程距離の長さを知ってほしいからです.
もっとちゃんと数学の説明してほしいというあなた,
もうちょっと我慢してください.

今回もアンケートがあるのでぜひ回答をお願いします.

• https://goo.gl/forms/1a7OmsfetjJCgVfB3

ではまた次回をお楽しみに!

プログラミングパート

放射性物質の崩壊

一番単純でしかも実際に使われる微分方程式として,
まずは 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

参考

コンピュータで計算するとき,
いろいろな都合があって分割数を多くすれば単純に精度が上がっていくと言うわけではない.
桁落ちや丸め誤差のような問題がある.
それはそれで別途検討する必要があるのだ.