2021年9月4日洛星高校土曜講座 Part1

最終更新：2021/09/14

土曜講座の内容（微分方程式のシミュレーション）を公開します. 多くを板書で説明したので要約となっていますが, （特に数式については）できるだけ詳しく書くようにします.

思ったより長くなったので, 2つにわけます. Part2は近日公開.

目的

この講座のテーマに微分方程式のシミュレーションを選んだ理由を簡単に説明します.

計算科学は第3の科学

計算科学は, 理論, 実験に次ぐ第三の科学とも言われています. 理論, 実験は高校生にとって想像しやすいですし, そもそも高校のカリキュラムで習っている数学や理科の内容であるともいえるでしょう. 一方, 計算科学を高校で習うことは（少なくとも私が高校生のときは）ありませんでした. 土曜講座の趣旨の一つは将来の可能性を皆に模索してもらうことだと聞いたので, 今回は皆にとって馴染みが薄いと思われる計算科学をとりあげてみました.

ちなみに, 最近はデータサイエンスが第4の科学ともいわれています. 最初はこちらを取り上げようかと思いましたが, 統計や線形代数の説明に時間を取りそうだったのでやめました. しかしながら, 今回講義した微分法手式の数値計算法などは, データサイエンス・機械学習でも不可欠な技術です.

シミュレーションを知ってもらう

最近は, コロナ対策の文脈で飛沫がどのように飛散するかのシミュレーションが話題になりました. 津波のシミュレーションなども馴染み深いものだと思います. これらも, 物理的な理論（数式）が先にありそれを計算機に解かしています. この, 「数式を計算機に解かす」というプロセスのイメージを持ってもらいたいと思いました.

大学で微分方程式を避けて通るのは難しい

少なくとも理系学部に進学すると, 常微分方程式・偏微分方程式には年次が低い時点でほぼ確実に出会うことになると思います. しかしながら, 見た目のイカツさから理解を諦めてしまったり, 手計算で解く過程を丸覚えするだけで, 自然現象との繋がりに気づかず終わっていしまう可能性があります. そこで, この講座では実際に大学レベルの偏微分方程式を取り上げ, シミュレーションを通して数学的理論ではなく物理的意味を強調することで, 数式への心理的障壁を取り除こうともいました.

数学的準備

ほんの少し数学的準備をおいておきます.

1変数関数の微分

関数 $u(x)$ の導関数が以下で定義されることは知っていると思います.

\begin{align}\large{ u^\prime (x) = \od{u}{x}{} = \lim_{\Delta x \rightarrow 0} \frac{u(x+\Delta x) - u(x)}{\Delta x} }\end{align}

このとき, 細かいことはおいといて

\begin{align}\large{ u^\prime (a) = \lim_{\Delta x \rightarrow 0} \frac{u(a+\Delta x) - u(a)}{\Delta x} }\end{align}

が, 関数 $f(x)$ の $x=a$ における接線の傾き, もしくは $f$ の変化率であるということを理解しておればOK. 例えば, $t$ に関する関数である速度 $v(t)$ の導関数（変化率） $v^\prime (t)$ は加速度でした.

2変数関数の偏微分

このノートでは, 高校レベルを超えた偏微分が出てきますが, 以下のように簡単に理解しておればOK.

2変数関数 $u(x, t)$ に対して偏導関数が以下のように定義されます.

\begin{align} \large \pd{u}{x}{} &\large= \lim_{\Delta x \rightarrow 0} \frac{u(x+\Delta x, t) - u(x, t)}{\Delta x}, \label{pdx} \\[5pt] \large \pd{u}{t}{} &\large= \lim_{\Delta t \rightarrow 0} \frac{u(x, t+\Delta t) - u(x, t)}{\Delta t}. \label{pdt} \end{align}

それぞれ, $u(x, t)$ の $x$ に関する変化率, $u(x, t)$ の $t$ に関する変化率と認識していれば充分です.

厳密には, それぞれ「$x$ を固定したときの変化率」, 「$t$ を固定したときの変化率」というようにい但し書きつきの変化率です. 図を見たほうがわかりやすいと思うので, 興味がある人は「偏微分傾き」などで画像検索してみて下さい.

簡単な例: 空気抵抗のある自由落下

モデル

もう早速計算してみましょう. 質量 $m$ の物体が自由落下する際の運動方程式が $ma = mg$ であることは既知だと思います. しかし, 現実には空気抵抗があるので運動方程式は以下のようになります.

\begin{align}\large{ m \od{v}{t}{} = mg - k v. \label{eom} }\end{align}

ここで, $v$ が時間 $t$ の関数 $v(t)$ であることに注意して下さい. また, $k$ は空気抵抗係数と呼ばれます. この式が表していることは「自由落下なのでもちろん速度はどんどん大きくなるが, それに比例して空気抵抗力 $kv$ も大きくなる」ということです.

シミュレーション結果

結果から見てみましょう. 初速度 $0\;(v(0) = 0)$ で物体を落したとき, どのように速度 $v$ が時間発展していくのかを上の式を用いて計算してみました. なお, $m = 0.3$, $k = 0.24$, $g = 9.81$ としています.

見ての通り, 最初はどんどん速さが大きくなっていきますが, だんだん空気抵抗の影響が大きくなり, 最終的には一定の速度 $v_\infty$ に落ち着いています. 式\eqref{eom}に基づいて議論すると, 一定の速度に到達したとき左辺の加速度は $0$ であるので, $v_\infty$ は $mg = k v_\infty$ を解いて得られます. つまり, $v_\infty = 0.3 \times 9.81 / 0.24 = 12.26\; \mathrm{[m/s]}$. なお, $v_\infty$ は, 終端速度と呼ばれます.

どうやって解いたのか : 離散化

ところで, この結果はなめらかな曲線に見えますが, 計算機で解いている以上本当に連続的な解を求めることはできません. 実際には, 以下のようにデータ（解）は点的なもので, それぞれの点の間を直線で結んでなめらかに見せているのです.

本当に滑らかな曲線を図示するには, 無限個の点を用意しなければならず, そのためには無限回の計算や無限のメモリが必要となるので不可能なのです. そこで, 計算機（コンピュータ）で方程式を解くときには, 飛び飛びの値として計算しなければなりません. 連続的な式\eqref{eom}を離散的な式, コンピュータで解ける式に変換する作業を離散化といいます.

さっそく離散化してみましょう. そのために, 速度に以下のように番号を振ります. このとき, 各点の時間間隔を時間刻み $\Delta t$ といいます. 今の場合, $\Delta t = 0.1 $ です.

これは, 上の図の原点付近を拡大したものですが, まだ $v_i$ の値はわかっていないと思って下さい. ただし, 初速度は $v_1 = 0$ とわかっています. なので, $v_1$ から $v_2$ を求め, $v_2$ から $v_3$ を求め・・・と逐次的に計算すればよいのです. 一般化すると, $v_{i}$ から $v_{v+1}$ を求めるための式を作ればよいということです. もちろんその際, 式\eqref{eom} を使います. 少し変形して以下のように書きましょう.

\begin{align}\large{ \od{v}{t}{} = g - \frac{k}{m} v. \label{eom2} }\end{align}

左辺は, 速度の変化率でした. 則ち, $v_i$ 付近で左辺は,

\begin{align}\large{ \od{v}{t}{} \simeq \frac{v_{i+1}-v_{i}}{\Delta t} \label{zensin} }\end{align}

と近似することができます. 以下の図も参考にしてみて下さい.

これをみて,

\begin{align}\large{ \od{v}{t}{} \simeq \frac{v_{i}-v_{i-1}}{\Delta t} \label{koutai} }\end{align}

でもいいのじゃないか？と思うかもしれません. 実際それは正解で, そのように近似することもあります.

これを使うと, 式\eqref{eom2}は, 以下のようにかけて

\begin{align}\large{ \frac{v_{i+1}-v_{i}}{\Delta t} = g - \frac{k}{m} v_i \label{eomdisc} }\end{align}

はれて$v_{i+1}$ と $v_i$ の関係式

\begin{align}\large{ v_{i+1} = v_{i} + \Delta t \left( g - \frac{k}{m} v_i \right) \label{eomdisc2} }\end{align}

が得られました. あとは, $v_1 = 0$ とわかっているのだから

\begin{align} \large v_{2} &\large= v_{1} + \Delta t \left( g - (k/m) v_1 \right) \\[5pt] \large v_{3} &\large= v_{2} + \Delta t \left( g - (k/m) v_2 \right) \notag \\[5pt] \large v_{4} &\large= v_{3} + \Delta t \left( g - (k/m) v_3 \right) \notag \\ &\qquad\qquad \vdots \notag \end{align}

と解いていけば, 電卓でもできます.

コード

最初に示した計算結果は, 実際に式\eqref{eomdisc2}をMATLABというプログラム言語を用いて解いたものです. 以下にコードを示します. プログラムに馴染みが無くても, なんとなくはわかると思います.

MATLAB

% 入力
DT     = 0.1;   % 時間刻み [sec]
TMAX   = 10;    % 計算時間 [sec] 
G      = 9.81;  % 重力加速後 [m*s^(-2)]
M      = 0.3;   % 質量 [kg]
K      = 0.24;  % 空気抵抗係数 [kg/m]
V0     = 0;     % 速度の初期値

% 計算
T = 0:DT:TMAX;         % 計算範囲 ( 時間 ) 
v = zeros(size(T));    % 事前割り当て 
L = length(T);
v(1) = V0;

for i = 1:L-1
    v(i+1) = v(i) + DT * (G - (K/M) * v(i));
end

ここで, DTを小さい値にすると, より正確な結果が得られますが, 同時に計算時間も増加します. 極端な話, DT = 0 とすれば, 微分係数の定義そのものとなりますが, 無限の計算時間がかかります.

2変数の例: 熱方程式

モデル

さて, 次は2変数の例を解いてみましょう. このモデルは, $u(x, t)$ という $x$, $t$ に関する2変数関数として表される量 $u$ の挙動を示し, 以下の式で表されます.

\begin{align}\large{ \pd{u}{t}{} = c \pd{u}{x}{2} \label{heateq} }\end{align}

先程は，$u$ が時間 $t$ に関する1変数関数 $u(t)$ でした. そして離散化し, $i$ 番目の時間における $u$ の値, つまり$u(i\Delta t)$, を $u_i$ と表したのでした.

今回は, $u$ が $x$ と $t$ との2変数（わかりにくければ, $x$ が位置を, $t$ が時間をあらわすと思えば良い）あるので, それにあわせて $x$ と $t$ とを離散化する必要がある. $i$ 番目の $x$, $j$ 番目の $t$ における $u$ の値, つまり $u(i\Delta x, j\Delta t)$, を $u_{ij}$ と表そう. そして, まず左辺を離散化してみます.

\begin{align}\large{ \pd{u}{t}{} \simeq \frac{u_{i, j+1} - u_{i, j}}{\Delta t} \label{dudt} }\end{align}

ここで, $i$ つまり $x$ 座標の値が動いていないことに注意して下さい. これは, 偏微分の定義\eqref{pdt}に含まれているのでした.

次に, 右辺を離散化してみます. 少し複雑ですが結論からいうと

\begin{align}\large{ c \pd{u}{x}{2} \simeq c \frac{u_{i+1, j} - 2u_{i, j} + u_{i-1, j}}{{\Delta x}^2} }\end{align}

となります. 以下に導出を書きます. 簡単です.

導出（クリック）

\begin{align} \pd{u}{x}{2} \simeq \frac{u_{i+1, j} - 2u_{i, j} + u_{i-1, j}}{{\Delta x}^2} \end{align}

を示そう. まず,

\begin{align} \pd{u}{x}{2} = \pd{}{x}{} \pd{u}{x}{} = \pd{U}{x}{}, \quad U = \pd{u}{x}{} \end{align}

とおきます. すると, 式\eqref{koutai}を参考に,

\begin{align} \pd{u}{x}{2} = \pd{U}{x}{} \simeq \frac{U_{i, j} - U_{i-1, j}}{\Delta x}, \label{eq2} \\[5pt] \end{align}

であることがわかります. さらに, \eqref{zensin}より,

\begin{align} & U_{i,j} = \pd{u_{ij}}{x}{} \simeq \frac{u_{i+1, j} - u_{i, j}}{\Delta x}, \\[5pt] & U_{i-1,j} \simeq \frac{u_{i, j} - u_{i-1, j}}{\Delta x}, \notag \end{align}

です. これらを, \eqref{eq2}に代入することで, 求めたい関係式が得られます.

$\blacksquare$

これで, 熱方程式\eqref{heateq}を離散化することができました.

\begin{align}\large{ \frac{u_{i, j+1} - u_{i, j}}{\Delta t} = c \frac{u_{i+1, j} - 2u_{i, j} + u_{i-1, j}}{{\Delta x}^2}. }\end{align}

これを変形して,

\begin{align}\large{ u_{i, j+1} = u_{i, j} + c\frac{\Delta t}{{\Delta x}^2} (u_{i+1, j} - 2u_{i, j} + u_{i-1, j}) \label{heatdisc} }\end{align}

とします. ここで, 境界条件を決めないといけませんが, 少し複雑なので $u_{ij}$ を以下のような格子として考えてみます.

すると, 式\eqref{heatdisc}は赤丸●の値を, 3つの黒丸●から求めていることが分かると思います. この方法で全ての格子点の値を求めようとすると, 赤線で示される境界の値を与える必要があります. いま,

$u_{1, 1}$ ~ $u_{M, 1}$ を時間的に一番最初の状態を示すので, 初期条件,
$u_{1, 1}$ ~ $u_{1, N}$ を空間的に端っこの状態を示すので, 境界条件,
$u_{M, 1}$ ~ $u_{M, N}$ も空間的に端っこの状態を示すので, 境界条件,

と呼ぶことにします. 今回, 境界条件は30, 初期条件は以下のグラフで与えられる値とします.

ここで, どういう意図でこの初期・境界条件にしたのかを説明します. タイトルに有る通り, このシミュレーションは熱方程式と呼ばれる偏微分方程式の数値解析であり, 温度がどう変化していくかを表しています. 今回想定したのは, 以下のようにフライパンを半分に切った断面の温度分布です.

ちょっと異様ですが, 初期条件は直径1mあるフライパンに火をかけ終わった瞬間の状態です. 火が掛かっていた真ん中は100度の温度となっており, （これもまた非現実的ですが）フライパンの端っこは室温である30度になっています. 式\eqref{heatdisc}を解くことで, 時間が立つとこのフライパンの温度がどう変化していくかがわかります.

結果

まず, 各パラメーターを以下のように設定します.

$L = 1$（空間方向の最大値 [m], フライパンの例でいうと直径）
$\Delta t = 0.002$（時間刻み）
$N = 5000$（時間刻みの数, 即ち $N \times \Delta t = 10$ 秒後まで計算）
$\Delta x = 0.02$（空間方向の刻み幅 [m]）
$c = 0.07$（定数 $c$）

結果を見てみましょう.

どうですか? 予想通りでしたか? 時間が立つにつれて, 中央部の温度が下がり, 最終的には室温である30度に全体が収束していることがわかります. これは, 直感や経験に沿った結果であると思います.

では, なぜ式\eqref{heateq}でこのような温度変化が表現できるのでしょうか? それを理解するために, 式\eqref{heateq}の右辺の $u$ の $x$ に関する二階微分の値もプロットしてみましょう.

高校2年ならすんなり理解できると思うけれど, 一応説明すると, $x$ に関する二階微分とは, 俗に言う「変化の割合」です. 初期状態（サムネイルの状態）の通り, グラフのてっぺんは急上昇から旧減少に転じているので二階微分は大きな負の値となっています. 逆に, その左側ではほとんど一定（30度）から急上昇, 旧減少からほとんど一定に転じており, 二階微分は大きな正の値となっています. 御存知の通り, 二回微分が $0$ となるところが変曲点であり, 普段君等が描いている増減表を思い出してもらうと合点がいくでしょう.

これを踏まえると, 式\eqref{heateq}（再掲します）

\begin{align}\large{ \pd{u}{t}{} = c \pd{u}{x}{2} \notag }\end{align}

は, $c$ が正であるとき,

温度の時間変化（左辺）は温度の空間的な2階微分（変化の割合）に比例するといっていることがわかります.

実際, 上の動画で温度が最も高いところは大きな負の二階微分値をもつので, 急激に温度が冷めていきます. しかしながら, 時間が立つにつれ温度分布が滑らかになってゆき, 二階微分の値が小さくなるので温度の変化スピードも小さくなっていくことがわかります.

余談

上のセクションでは, 結果を時系列でみせたのでいかにも温度の時間変化に見えました. しかしながら, 行った計算は式\eqref{heateq}で表される2変数の偏微分方程式を初期値・境界条件を満たすように解いたにすぎません. 得られた解を一気にみると以下のようになります.

時系列の動画で見ていたのは, $x$ 軸に平行にこの表面を輪切りにしたものを連続的に表示していたのです.

函数の不定積分を求めると, 積分定数がでてきます. つまり, 積分定数で表されるだけの解の不定性があるのでした. もし, 積分定数を求めて解を確定させたければ, ある点における函数の値のような条件を与えなければなりません.

同様に, 偏微分方程式の解にも, 一般的に不定性があります. なぜなら, \eqref{heateq} をみればわかるように, やっていることは積分を行って $u$ を求めているからです. そして, 解を一つに決めるために（今回の場合上の図の曲面を決めること）, 境界・初期条件を与える必要があるのです.

コード

MATLABのコードを載せます. あまり詳しく説明していません（でもこれだけのコードで計算できることは驚いてもらいたいです）. また, 二階微分の値を描くコードや動画を作るコードは作っているうちにぐちゃぐちゃになったので載せていません. ですが, 希望があれば載せますので連絡ください.

MATLABコード（クリック）