たのしい！数値計算
この記事は、奈良高専 Advent Calendar 2022の24日目の記事です。前回の方(Jin)のブログはこちらです。
はじめに
皆さん初めまして。ぴかです。
最近BLEACHの千年血戦篇を見ています。めちゃくちゃおもろいですね。
↑バンビエッタちゃん。かわいい。
さて、今回は数値計算をテーマになんか書いていこうと思います。なんで数値計算かというと、クリスマスだからです*1。クリスマスは数値計算をしろ。
わたくしの分際でクリスマスイブの記事を書くのは、たいへんおこがましいと思いつつも、なんかごちゃごちゃと書いていこうかなといった次第であります。
目次？
たいそうな文章ではないため、目次とは言えない何かです。うーん、おしながきといったところでしょうか。
数値計算ってなに～
数値計算ってなんなんでしょうね。数値を使って計算するなら算数じゃね？って思われるかもしれません。デジタル大辞泉によると、
コンピューターを使って、代数的または解析的に解くことが困難な問題、または手計算では不可能な計算量が膨大な問題を解くこと。
とのことです。
すなわち、手計算するのがしんどかったり、解析的（大雑把に言うと式変形）に解くのが難しいものを計算する方法が数値計算法というわけですね。
数値計算の対象
数値計算だったり数値解析の対象となるものはめっちゃ多いです。
数値補間、行列式、微積分、微分方程式（常微分も偏微分も）とかその他たくさんあります。
その応用例は、流体解析とか電磁界解析、ウイルスの感染*2をシミュレーションしたりとかです。身近なものだと天気予報もあります。
実際に数値計算してみる
実際に数値計算してみます。むずかしいモデルを数値計算するには、その分プログラム書くのも大変だったりします。
そのため、ここでは古典力学上の簡単な計算例を紹介します。
自由落下問題
質量をもった物体（ここでは質点とします）がある高さから落下していくモデルを例に数値計算します。
↑スカイダイビングする人
まず、自由落下について理論式から考えていきます。質量を$ m $、速度を$ v(t) $、重力加速度を$ g $とすれば、Newtonの運動方程式は、
$$ m\frac{dv(t)}{dt} = mg \tag{1} $$
となります。
この微分方程式を解けば、物体が移動する距離や物体の速度が求まるわけです。
ところで微分ってどうやって数値計算するんでしょうか？？
微分はTaylor展開により表すことができます。
$ t = \Delta t $まわりのTaylor展開により$ \frac{dv(t)}{dt} $は、
$$ \frac{dv(t)}{dt} \approx \frac{v(t + \Delta t )-v(t)}{\Delta t} \tag{2} $$
となります。（これは、$ \Delta t \approx 0 $とすれば微分の定義$ \frac{dv(t)}{dt} = \lim_{\Delta t \rightarrow 0} \frac{v(t + \Delta t )-v(t)}{\Delta t} $と一致し、微分を近似表現できていると言えます。）
ここで、速度の微分を$ \frac{dv(t)}{dt} = f(t,v) $とすると、式(2)は、
$$ v(t+\Delta t) = v(t) + f(t,v) \cdot \Delta t \tag{3} $$
こんな感じになります。同様に、移動距離$ x(t) $についても、
\begin{align} x(t+\Delta t) &= x(t) + \frac{dx(t)}{dt} \cdot \Delta t \\ &= x(t) + v(t) \cdot \Delta t \tag{4} \end{align}
と表すことができます。
これら式(3)、(4)により微分方程式の数値解を計算する方法をEuler法といいます。
このEuler法によって式(1)で表される自由落下を数値計算します。なおPythonで数値計算しました。Pythonは個人的に使いやすいなっていうのと、colaboratoryとかweb上で実行できたりするので便利です（VScodeのほうが使いやすいけど）。
コードはこんな感じ（きたないかもだけど許して。あとスマホだとこのコード表示が見にくくなるんやが、どうしたらいいんでしょうか。）。
#簡単な力学系のシミュ(1次元質点運動問題) #ここでは自由落下を扱う(空気抵抗を一切含んでいない F = mg の式) import numpy as np import matplotlib.pyplot as plt #const g = 9.8 #計算エンジン------------------------------------------------------------ #変数 t = 0.0 #時刻t d = 0.01 #微分刻み幅d #格納配列(csv出力とかグラフとかにするため) tlist = [] ylist = [] vlist = [] #入力部 y = float(input("初期高度y0を入力してほしいにゃ:")) #初期位置を入力 v = float(input("初速度v0を入力してほしいにゃ:")) #初速を入力 print("{:.5f} {:.5f} {:.5f}" .format(t,y,v)) #入力時の数値出力 #自由落下計算 while y >= 0: #高度が0になるまでの時間 t += d #時刻更新 v += g*d #速度更新 y -= v*d #位置更新 print("{:.5f} {:.5f} {:.5f}" .format(t,y,v)) #計算課程から終了まで数値出力 tlist.append(t) ylist.append(y) vlist.append(v) #グラフ出力します figy = plt.figure(1) plt.plot(tlist,ylist) plt.xlabel('Time [s]') plt.ylabel('Height [m]') plt.show() figv = plt.figure(2) plt.plot(tlist, vlist) plt.xlabel('Time [s]') plt.ylabel('Velocity [m/s]') plt.show()
初期位置を50 m、初速を0 m/sとして実行します。なお計算条件は、時刻の刻み幅を0.01 s、質点の高度が0になる時刻で計算終了としました。
質点の高度変化をFig3に、速度変化をFig4に示します。
わ～、自由落下運動の数値解を求めれた～......
実はこれ、空気抵抗を考慮していないモデルなので、初期高度が高ければ高いほど落下時における質点の速度は速くなります。
例えば雨を降らす乱層雲は地上から1 km～2 kmぐらいの高さですが、このとき空気抵抗が無いと雨粒は140 ～198 m/sの速度で降ることになります（雨が10 kgの量降るとしたら）。
時速換算すると、504～713 km/h程度です（この速度で雨が降ればヤバイですよね）。
ということで空気抵抗を考慮した自由落下を考えましょう。
閑話休題
空気抵抗を考慮する前にちょっとだけ。
最近Polyphiaっていうバンドの曲をよく聞いてるんですが、まぁえぐいんですね。インストバンドなんですが、ギターがうますぎる。ギターはTim HensonとScott Lepageの2人なんですが、2人ともえげつうまい。
Playing Godとかえげついです。クラシックギターでピロピロ弾いてるし。ぜひ聞いてみてね。
空気抵抗を考慮した自由落下
物理の問題とかで速度に比例する空気抵抗とかよくありますよね*3。アレを考えます。
速度$ v $に比例する空気抵抗を考慮した場合、運動方程式は、
$$ m\frac{dv(t)}{dt} = mg - kv(t) \tag{5} $$
と表すことができます。ここで$ k $は空気抵抗の比例係数です。
この比例定数は粘性抵抗と言われるぽいです。
式(5)について初速を0とし、速度$ v(t) $を導くと、
$$ v(t) = \frac{mg}{k}\Bigr ( 1 - \exp\bigr [ -\frac{kt}{m} \bigl ] \Bigl ) \tag{6} $$
となります。
時間$ t $が十分長いとき、速度は、
$$ \lim_{\Delta t \rightarrow \infty} v(t) = \frac{mg}{k} \tag{7} $$
となります。
この速度を終端速度といい、質点は一定の速度に達します。
さて、この空気抵抗を考慮した自由落下を数値計算しましょう。
コードはこんな感じ。
#簡単な力学系のシミュ(1次元質点運動問題) #速度に比例する空気抵抗(粘性抵抗)がかかる自由落下を扱う(F = mg - kv の式) import numpy as np import matplotlib.pyplot as plt #const g = 9.8 #計算エンジン------------------------------------------------------------ #変数 t = 0.0 #時刻t d = 0.01 #微分刻み幅d #格納配列(csv出力とかグラフとかにするため) tlist = [] ylist = [] vlist = [] #入力部 m = float(input("質量mを入力してほしいにゃ:")) #物体の質量を入力 k = float(input("空気抵抗kを入力してほしいにゃ:")) #空気抵抗を入力 y = float(input("初期高度y0を入力してほしいにゃ:")) #初期位置を入力 v = float(input("初速度v0を入力してほしいにゃ:")) #初速を入力 print("{:.5f} {:.5f} {:.5f}" .format(t,y,v)) #入力時の数値出力 #自由落下計算 while y >= 0: #高度が0になるまでの時間 t += d #時刻更新 v += (g-k/m*v)*d #速度更新 y -= v*d #位置更新 print("{:.5f} {:.5f} {:.5f}" .format(t,y,v)) #計算課程から終了まで数値出力 tlist.append(t) ylist.append(y) vlist.append(v) #グラフ出力します figy = plt.figure(1) plt.plot(tlist,ylist) plt.xlabel('Time [s]') plt.ylabel('Height [m]') plt.show() figv = plt.figure(2) plt.plot(tlist, vlist) plt.xlabel('Time [s]') plt.ylabel('Velocity [m/s]') plt.show()
質量を10 kg、比例定数を10 kg/s、初期位置を50 m、初速を0 m/s、として実行します。なお計算条件は、時刻の刻み幅を0.01 s、質点の高度が0になる時刻で計算終了としました。
質点の高度変化をFig7に、速度変化をFig8に示します。
ちゃんと速度が一定のとこ（およそ9.8 m/s）に漸近していますね。
質点の高度変化とか速度変化をグラフだけで見ていても退屈になってきたので、アニメーションで質点の動きを見てみましょう。
動き始めよりも地面に近づくときのほうが速い落下になってますね。
おわりに
自由落下という簡単な物理現象を数値計算してシミュレーションしました。わりとおもろいですよね。
微分方程式で表される現象は古典力学の運動だけでなく、電磁界を表すMaxwell方程式や流体のNavier-Stokes方程式であったりと、たくさんあります。（私は電磁界解析を研究でやってたりやってなかったり）
気になった方はぜひ数値計算してみてください。
おまけ
アニメーションを含むコードです。
# 簡単な力学系のシミュ(1次元質点運動問題) # 速度に比例する空気抵抗(粘性抵抗)がかかる自由落下を扱う(F = mg - kv の式) import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation # const g = 9.8 # 計算エンジン------------------------------------------------------------ # 変数 d = 0.01 # 微分刻み幅d t = 0.0 # 時刻 # 格納配列(csv出力とかグラフとかにするため) tlist = [] ylist = [] vlist = [] td = np.arange(0, 100, 0.01) # アニメ時刻 # 入力部 m = float(input("質量mを入力してほしいにゃ:")) # 物体の質量を入力 k = float(input("空気抵抗kを入力してほしいにゃ:")) # 空気抵抗を入力 y = float(input("初期高度y0を入力してほしいにゃ:")) # 初期位置を入力 v = float(input("初速度v0を入力してほしいにゃ:")) # 初速を入力 print("{:.5f} {:.5f} {:.5f}" .format(t, y, v)) # 入力時の数値出力 y0 = y # 自由落下計算 while y >= 0: # 高度が0になるまでの時間 t += d # 時刻更新 v += (g-k/m*v)*d # 速度更新 y -= v*d # 位置更新 print("{:.5f} {:.5f} {:.5f}" .format(t, y, v)) # 計算課程から終了まで数値出力 tlist.append(t) ylist.append(y) vlist.append(v) t_posi0 = t # 地面に到着する時間 # 計算 ind_0 = np.argmin(np.abs(td-t_posi0)) ylist = ylist[:ind_0] # データをグラフに(ここでは高度変化) data = np.concatenate([ylist]) plt.plot(data) plt.ylabel("Height [m]") plt.xlabel("Time [s]") plt.show() data.shape # アニメーション fig, ax = plt.subplots(figsize=(6, 4)) ball, = ax.plot(0, data[0]+5, 'C2o', ms=15) ax.hlines(0, -0.25, 0.25) ax.set_xticks([]) ax.set_ylim(-10, data.max()+20) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.spines['bottom'].set_visible(False) ax.set_title("$Height$="+str(y0)+"[m]") def update(num): ball.set_ydata(data[num]+5) ani = animation.FuncAnimation(fig, update, 1000, interval=10) ani.save('ball.gif', dpi=100)
一昨日、全然寝れてなくて執筆・投稿が少し遅くなりました。ごめんなさい🎅
ちなみに「たのしい！数値計算」ってタイトルですが、アニメーションとか見てるときはたのしいけど、それ以外なんもたのしくなかったです、ねむすぎて。
