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

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

こちらに PDF があります.
サイトでは見づらい方, コンテンツを手元に置いておきたい方は
ダウンロードしてご覧ください.

数列も関数

ベクトル, 関数と来て今回は数列です.
簡単に復習すると,
関数は数と数の割り当てルールのことで,
数列も関数なのでした.

等比数列

いろんな数列があります.
高校でもいろいろやります.
今回の話で 1 番大事な等比数列です.
長々とやっても仕方ないので軽く復習しときましょう.

$(x_n)$ を数列とします.
高校だと {$x_n$} と中括弧で書くと思いますが大学の数学的な都合で小括弧で書くことにします.
等比数列がどんな数列だったか,
つまりどんなルールで与えられているのかというとこんなルールです:
\begin{align}
x_{n+1}
=
\alpha x_{n}.
\end{align}
$x_{n}$ が 0 のときでも成り立つ形で書きました.
$x_{n}$ が 0 ではないならもちろん $x_{n+1} / x_{n} = \alpha$ と書けます.
どんな $n$ に対しても比が一定な数列だから等比数列というのでした.

等比数列の一般項

このルールを使って一般項を出します.
$x_{n+1} = \alpha x_{n}$ で,
$x_{n} = \alpha x_{n-1}$ で,
$x_{n-1} = \alpha x_{n-2}$ で,
と延々続きます.
これが終わるのは $x_1 = \alpha x_{0}$ です.
順々に代入していくと $x_{n} = x_{0} \alpha^{n}$ になりますね.

この計算, はまる人ははまります.
実際私も高校の頃, 時々わけがわからなくなることがありました.
試験で慌てていてパニクったことを思い出します.
それはそれとして,
ここで細々とした計算を詳しく解説してると本筋を見失うので省略します.
数列をどんなところでどう使うのか,
それをちゃんと見てからきちんとやった方がよさそうですしね.

微分方程式と等比数列

これがどこで出てきたかというと微分方程式を近似した方程式で出てきたのでした.
放射性物質の崩壊で出てきたのはこんなやつです.
\begin{align}
\frac{x_{n+1} – x_{n}}{h}
=
– x_{n}.
\end{align}
これを整理すると $x_{n+1} = (1 – h) x_{n}$ で,
$1-h = \alpha$ とすればまさにさっき書いた等比数列ですね.
さっき書いたように等比数列は指数関数で書けます.
で, 元の微分方程式の解も実は指数関数で書けます.

もう一度: 数列も関数

ちょっと先走ったことも書きました.
何はともあれ, 数列は関数で, 関数は数と数の割り当てルールのことで,
その割り当てルールとして等比数列っていう割り当て方があるんだ,
ということです.
微分方程式もその割り当てルールと見なせるっていうのが今書いたことです.
微分方程式と揃えた言い方として上の式は差分方程式と言うこともあります.

漸化式

ちなみに上で紹介した割り当てルールはもっと一般的にできます.
これ, 要は次の項をその前の項で決まると言っているだけなので,
$a_{n+1} = a_{n}^2$ みたいにしてもいいし,
$a_{n+1} = a_{n} + a_{n-1}$ みたいに 2 つ前と関係あるようにしてもいいです.
この一般的なやつがいわゆる漸化式です.

2 つの数列に対する漸化式

もちろんこれ以外にもいろいろな漸化式があります.
あとで使うのでもう 1 つ漸化式を紹介します.
そのために 2 つの数列 $(x_{n})$, $(v_{n})$ を用意します.
で, 次のような漸化式を考えます.
\begin{align}
x_{n+1} – x_{n}
=
v_{n}, \quad
v_{n+1} – v_{n}
=
– \alpha^2 x_{n}.
\end{align}
2 つの数列が絡まりあう漸化式です.
高校でやる言葉を使うなら,
互いに階差数列を取ってると思ってもいいかもしれません.

単振動の漸化式

これがどこから出てくるかというと次の微分方程式からです.
\begin{align}
\frac{d^2 x(t)}{dt^2}
=
– \alpha^{2} x(t).
\end{align}

この微分方程式も物理でよく出てくる方程式で,
いわゆる単振動の運動方程式です.
さっきの漸化式を出すには $v(t) = x'(t)$ としてもとの方程式も $v’ = – \alpha^2 x$ とすればいいです.
これを差分近似すると次の式が出てきます.
\begin{align}
x_{n+1}
=
x_{n} + h v_{n}, \quad
v_{n+1}
=
v_{n} – h \alpha^{2} x_{n}.
\end{align}
$h$ が入っていて係数が微妙に違いますが,
あまり気にしないでください.

この辺の計算プログラムも準備してあります.
この記事後半のプログラムパートを確認してください.

アンケートの回答をお願いします

今回もアンケートがあります.
改善につなげるためぜひ回答をお願いします.

• https://goo.gl/forms/gJSiOBd9vUa9D0Xr2

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

プログラミングパート

1 階の線型常微分方程式

復習でもう 1 回.

一番単純でしかも実際に使われる微分方程式としてまずは 1 階の線型常微分方程式を考えよう.
ちょっと不吉な例であるが放射性物質の崩壊の方程式を紹介する.
導出をしたければちゃんと物理を勉強してもらう必要がある.
ここでは物理は省略して数学に集中する.

\begin{align}
\frac{dx}{dt} = – c x.
\end{align}
厳密解は $x = 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}

これに沿って計算したのがいわゆるオイラー法.
次のセルではこれをコードに落としている.

1 階の方程式のプログラム

%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
nt = 101
init = 5
x1 = radioactive_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x1)

 # 厳密解
t = np.linspace(0, 2, nt)
x2 = init * np.exp(- c * t)
plt.plot(t, x2)

 # 凡例
plt.legend(['approximation', 'exact'])
plt.show()
 # 描画
#plt.savefig(img_file)
#b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
#img_str = "" % (b64)
#print(img_str)

RESULT

2 階の常微分方程式: 単振動

こちらは単振動.

次に 2 階の常微分方程式を紹介しよう.
高校の物理で出てくるばねの振動（単振動）がまさにこの例だ.
項を増やすと減衰振動になったり、外力をつけたりといろいろなケースがある.
まずは一番単純な式を考えよう.

\begin{align}
\frac{d^2 x}{dt^2}
=
– \omega^2 x.
\end{align}

オイラー法やルンゲ-クッタ法は
1 階の方程式に対する計算法なので直接は使えない.

今は中間処理として $v = dx/dt$ を置いて計算すればいい.
これは単なる数値計算の便法ではない.

速度の意味もあるから,
という表面的な理由ではなくもっと深く解析力学の文脈で物理としても大事な視点だ.
もっといえばシンプレクティック計算法などもっといい計算法にも発展する.

オイラー法

とりあえずオイラー法で計算したい.
まずは微分方程式自体を書き直す.

\begin{align}
\frac{dx}{dt}
=
v, \quad
\frac{dv}{dt}
=
– \omega^2 x.
\end{align}

これをオイラー法で近似しよう.

\begin{align}
x_{n+1}
=
x_{n} + h v_{n}, \quad
v_{n+1}
=
v_{n} – h \omega^{2} x_{n}.
\end{align}

オイラー法をコードに落とす.

オイラー法のプログラム

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def harmonic_euler(nt, init = (5, 0)):
    dt = t_range / (nt - 1)
    # 初期条件設定
    x = np.zeros(nt)
    v = np.zeros(nt)
    x[0] = init[0]
    v[0] = init[1]

    for i in range(1, nt):
        x[i] = x[i-1] + dt * v[i-1]
        v[i] = v[i-1] - dt * (omega ** 2) * x[i-1]

    # ベクトル計算で書き直したい

    return (x, v)


omega = 2 * np.pi
nt = 101
t_range = 2
init = (5, 0)
harm = harmonic_euler(nt, init)
t = np.linspace(0, 2, nt)

 # 厳密解
x_exact = init[0] * np.cos(- omega * t)
v_exact = - omega * init[0] * np.sin(omega * t)

 # グラフ描画
plt.subplot(3, 1, 1)
plt.title('x-t graph')
plt.plot(np.linspace(0, 2, nt), harm[0])
plt.plot(t, x_exact)
plt.legend(['x approximation', 'x exact'])

plt.subplot(3, 1, 2)
plt.title('v-t graph')
plt.plot(np.linspace(0, 2, nt), harm[1])
plt.plot(t, v_exact)
plt.legend(['v approximation', 'v exact'])

plt.subplot(3, 1, 3)
plt.title('phase space')
plt.plot(harm[0], harm[1])

 # 描画
plt.tight_layout()
plt.show()
#plt.savefig(img_file)
#b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
#img_str = "" % (b64)
#print(img_str)

RESULT

オイラー法では近似の精度が悪い

見ての通り時間が進むごとに誤差が大きくなる.
nt を大きくすると少しはましになる.
実際に上のコードで nt を大きくして再計算してみてほしい.

あとまずいのは phase space の図だ.
この系はエネルギーが保存する系だから相空間内の軌道が閉じてほしいのにそうなっていない.
シンプレクティックにやれば解消できるようだが, とにかくここではよろしくない.

この方程式でオイラー法はよろしくないことがわかった.
とりあえずルンゲ-クッタでやってみよう.

ルンゲ-クッタ法

とりあえず近似式を書く.

\begin{align}
x_{n+1}
&=
x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}),
\end{align}
\begin{align}
t_{n+1}
&=
t_{n} + h,
\end{align}
\begin{align}
k_{1}
&=
v_{n},
\end{align}
\begin{align}
k_{2}
&=
v_{n} + \frac{h}{2} k_{1},
\end{align}
\begin{align}
k_{3}
&=
v_{n} + \frac{h}{2} k_{2},
\end{align}
\begin{align}
k_{4}
&=
v_{n} + h k_{3}.
\end{align}

次が $v$ の式.

\begin{align}
v_{n+1}
&=
v_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}),
\end{align}
\begin{align}
t_{n+1}
&=
t_{n} + h,
\end{align}
\begin{align}
k_{1}
&=
– \omega^2 x_{n},
\end{align}
\begin{align}
k_{2}
&=
x_{n} – \frac{h}{2} \omega^2 k_{1},
\end{align}
\begin{align}
k_{3}
&=
x_{n} – \frac{h}{2} \omega^2 k_{2},
\end{align}
\begin{align}
k_{4}
&=
x_{n} – h \omega^2 k_{3}.
\end{align}

ルンゲ-クッタ法のプログラム

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

def harmonic_rk(nt, init = 10):
    dt = t_range / (nt - 1)
    # 初期条件設定
    x = np.zeros(nt)
    v = np.zeros(nt)
    x[0] = init[0]
    v[0] = init[1]

    def fx(t, x, v):
        return v

    def fv(t, x, v):
        return - (omega ** 2) * x

    # ベクトル計算で書き直したい
    for i in range(1, nt):
        xk1 = fx(dt * (i - 1), x[i-1], v[i-1])
        vk1 = fv(dt * (i - 1), x[i-1], v[i-1])

        xk2 = fx(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)
        vk2 = fv(dt * (i - 1/2), x[i-1] + xk1 * dt / 2, v[i-1] + vk1 * dt / 2)

        xk3 = fx(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)
        vk3 = fv(dt * (i - 1/2), x[i-1] + xk2 * dt / 2, v[i-1] + vk2 * dt / 2)

        xk4 = fx(dt * (i - 1),   x[i-1] + xk3 * dt,     v[i-1] + vk3 * dt)
        vk4 = fv(dt * (i - 1),   x[i-1] + xk3 * dt,     v[i-1] + vk3 * dt)

        x[i] = x[i-1] + dt / 6 * (xk1 + 2 * xk2 + 2 * xk3 + xk4)
        v[i] = v[i-1] + dt / 6 * (vk1 + 2 * vk2 + 2 * vk3 + vk4)

    return (x, v)


omega = 2 * np.pi
nt = 101
t_range = 2
init = (5, 0)
harm = harmonic_rk(nt, init)
t = np.linspace(0, 2, nt)


x_exact = init[0] * np.cos(- omega * t)
v_exact = - omega * init[0] * np.sin(omega * t)


plt.subplot(3, 1, 1)
plt.title('x-t graph')
plt.plot(np.linspace(0, 2, nt), harm[0])
plt.plot(t, x_exact)
plt.legend(['x approximation', 'x exact'])

plt.subplot(3, 1, 2)
plt.title('v-t graph')
plt.plot(np.linspace(0, 2, nt), harm[1])
plt.plot(t, v_exact)
plt.legend(['v approximation', 'v exact'])

plt.subplot(3, 1, 3)
plt.title('phase space')
plt.plot(harm[0], harm[1])

 # 描画
plt.tight_layout()
plt.show()
#plt.savefig(img_file)
#b64 = base64.encodestring(open(img_file, 'rb').read()).decode('utf-8')
#img_str = "" % (b64)
#print(img_str)

RESULT

ルンゲ-クッタ短評

今度の一致具合はなかなかよさそう.
あくまで見た目の感じではあるけれども.

良くなったのは一番下の図, 相空間軌道だ: ちゃんと閉じてくれた.