カテゴリー: メルマガ

こちらに 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

ルンゲ-クッタ短評

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

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

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

ここまでの復習

ここまでは微分方程式が自然現象のシミュレーションに使えることを紹介し,
実際の微分方程式をいくつか見てきました.
実際のシミュレーションでは本当に微分方程式を厳密に解いているわけじゃなくて,
近似計算をしているんだ,
そしてそれは四則演算なんだというのを大雑把に紹介してきました.

近似といってもけっこうきちんと計算できますよ,
ということを見るためにもとにかくいくつかの実例を紹介することを優先にしてきました.
我慢強いあなたも「もういい加減きちんと数学の紹介して!」
と思ってらっしゃるはずです.

ベクトルからはじめます

というわけで数学の紹介も簡単にやっていきましょう.
微分方程式なんだからまずは微分じゃないの?
と思う方もいらっしゃるでしょう.
でも今回の話の流れでは微分の話をする前にベクトルをやった方がスムーズなので,
ベクトルからいきます.

ベクトルは矢印

高校でベクトルは向きと大きさがある矢印のことと言われます.
それはそれで大事な見方ですが,
ここでは始点と終点がある矢印と思ってみてください.
始点と終点をつないだたくさんのベクトルで作った折れ線を方程式の近似解とみなすのが今回の話のキモだからです.

図のように曲線を適当な間隔でわけましょう.

わけた点に $A_0, A_1, \dots, A_n$ と順に名前をつけていきます.
そして $A_0$ と $A_1$, $A_1$ と $A_2$ というふうに始点と終点を入れかえながら隣の点どうしを結んでいきます.
$A_0$ から $A_1$ に向かって進むんだ,
という気持ちを表したいのでそれを $\overrightarrow{A_0 A_1}$ と書くことにします.
記号の上の右向き矢印で点 $A_0$ から点 $A_1$ に向かっている気分もはっきり書いています.

ここで何度も始点と終点をつないだ矢印を継ぎ足しています.
この継ぎ足しという図形の操作がベクトルの足し算です.
例を 1 つだけ書いておくと $\overrightarrow{A_0 A_1} + \overrightarrow{A_1 A_2} = \overrightarrow{A_0 A_2}$ です.

微分方程式のシミュレーションから見たベクトル

高校だとベクトルに関してもっといろいろなことが出てきます.
定数倍したり内積を取ったり長さを調べたりなどなど.
もちろんベクトルを使っていろいろやるならやっておかないと困ります.
でも微分方程式から数学を眺めてみる立場ではまずここさえ乗り越えればどうにかなります.
他のことは他のことをやるときに絡めてやってみてください.

折れ線をつないでいって曲線を近似する様子をもっときちんと見てみましょう.
具体的な曲線としては円を取り,
それを等分して折れ線近似した様子をプログラムで描きました.
次のページを見てください.
最後に近似がよくなっていくアニメーションもつけています.

• http://tinyurl.com/zhoptac

基本的なプログラムはこの記事の後半にも載せています。
アニメーションは Jupyter との連携でやっているので,
この記事ではうまく動かせません.
上のリンクのページから確認してください.

今回はここまでです.
お疲れ様でした.

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

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

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

プログラミングパート

円を折れ線近似する関数

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

div_num = 9
def draw_circle_by_polygonal_line(div_num):
    thetas = np.linspace(0, 2 * np.pi, div_num)
    xs = np.cos(thetas)
    ys = np.sin(thetas)

    plt.plot(xs, ys)
    plt.gca().set_aspect('equal', adjustable='box')

ベクトルでつないだ折れ線で円を近似する

円周上の点を $n$ 分割し折れ線で分点をつないでいった.
一旦個別に図示している.

for i in range(4, 12):
    pnum = i - 3
    plt.subplot(2, 4, pnum)
    plt.xlim([-1.2, 1.2])
    plt.ylim([-1.2, 1.2])
    plt.title("div = " + str(i-1))
    draw_circle_by_polygonal_line(i)
    draw_circle_by_polygonal_line(1000)

plt.tight_layout()

RESULT

ベクトルをつないだ折れ線で円を近似する

9 分割までを重ねて 1 つの図にしてみた.

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

for div_num in range(10):
    draw_circle_by_polygonal_line(div_num)

draw_circle_by_polygonal_line(1000)
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])
plt.show()

RESULT

アニメーションにしてみる

分割を細かくしていくことで精度が上がる様子をアニメーションで見てみよう.

下のプログラムを Jupyter 上で実行すると,
図の両端にある「+」「-」ボタンでアニメーションのスピードを変えられる.
適当に変えて遊んでみよう.

注意
次のプログラムは Jupyter 上で実行することを前提にしている.
他の実行環境で動くかどうかは確認していない.
ここを見ると Jupyter のアニメーションが見られる.

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

def circle_data(div_num):
    thetas = np.linspace(0, 2 * np.pi, div_num)
    xs = np.cos(thetas)
    ys = np.sin(thetas)
    return (xs, ys)


circles = []
for i in range(4, 31):
    circles.append(circle_data(i))


#fig = plt.figure(figsize=(10, 10));
fig = plt.figure();


ax = plt.axes(xlim=(-1.2, 1.2), ylim=(-1.2, 1.2));
line, = ax.plot([],[],lw=2);


plt.gca().set_aspect('equal', adjustable='box')

def animate(data):
    x = data[0]
    y = data[1]
    line.set_data(x, y)
    return line

anim = animation.FuncAnimation(fig, animate, frames=circles, interval=100)
display_animation(anim, default_mode='reflect')