数列とは何か?プログラム解説 中高数学駆け込み寺

この記事は6分で読めます

中高の数学の復習から専門的な数学・物理までいろいろな情報を発信しています.
中高数学に関しては中高数学駆け込み寺,
大学数学に関しては現代数学観光ツアーという無料の通信講座があります.
ご興味のある方はぜひお気軽にご登録ください!


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}

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


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

import base64
img_file = "tmp/image.tmp.png"

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}), \
t_{n+1}
&=
t_{n} + h, \
k_{1}
&=
v_{n}, \
k_{2}
&=
v_{n} + \frac{h}{2} k_{1}, \
k_{3}
&=
v_{n} + \frac{h}{2} k_{2}, \
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}), \
t_{n+1}
&=
t_{n} + h, \
k_{1}
&=
– \omega^2 x_{n}, \
k_{2}
&=
x_{n} – \frac{h}{2} \omega^2 k_{1}, \
k_{3}
&=
x_{n} – \frac{h}{2} \omega^2 k_{2}, \
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

ルンゲ-クッタ短評

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

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


中高の数学の復習から専門的な数学・物理までいろいろな情報を発信しています.
中高数学に関しては中高数学駆け込み寺,
大学数学に関しては現代数学観光ツアーという無料の通信講座があります.
ご興味のある方はぜひお気軽にご登録ください!

  • このエントリーをはてなブックマークに追加
  • LINEで送る

関連記事

  • コメント (0)

  • トラックバックは利用できません。

  1. この記事へのコメントはありません。

このサイトについて

はじめまして。相転移Pです。数学・物理の情報を中心にアカデミックな話題を発信しています。このサイトを見て興味があればぜひご連絡ください。 mail: phasetr@gmail.com LINE: oxg2753d
  • このエントリーをはてなブックマークに追加
  • LINEで送る

YouTube チャンネル登録

講義など動画を使った形式の方が良いコンテンツは動画にしています。ぜひチャンネル登録を!

メルマガ登録

メルマガ登録ページからご登録ください。 数学・物理の専門的な情報と大学受験向けのメルマガの 2 種類があります。

役に立つ・面白い記事があればクリックを!

記事の編集ページから「おすすめ記事」を複数選択してください。