はじめに¶
大気や海洋の運動は、ニュートンの第2法則に基づいた微分方程式で記述されます。
$$
\frac{d{\bf u}}{dt}=\frac{\bf F}{\rho}
$$
(${\bf u}$は速度、$t$は時間、$\rho$は密度、${\bf F}$は単位質量の流体に加わる力)
物体の運動方程式を解いて運動を求めるように、この方程式を解けば将来の空気・水の様子が予測できるわけです。
しかし、話はそう単純ではありません。予測すべき速度$u$は空間の関数であり、周囲の状況がその点の時間変化に影響する偏微分方程式だからです。
一部のごく単純な偏微分方程式はきれいに解が求まりますが、ほとんどの偏微分方程式は手に負えないほど複雑です。
それでも未来の天気を予測するために、1950年代以降コンピュータの力を借りて微分方程式を解く試みが続けられてきました。
その結実が今日の天気予報や気候予測なのです。
コンピュータは有限個の数値しか扱えません。数値シミュレーションとは、コンピュータが解けるよう微分方程式を近似表現し、元々の方程式を模擬する手法です。 ここでは数値シミュレーションの基本的手法を学び、体験しましょう。今回は偏微分方程式よりいくらか単純な常微分方程式の数値シミュレーションを、実際の現象を踏まえつつ学んでいきましょう。
慣性振動の数値シミュレーション¶
KEOデータの例¶
第2回・第3回で扱ったKEOには、水中の流れを音波で計測する音波式ドップラー流速計 (acoustic Doppler current profiler; ADCP) が設置されています。 ADCPは上向きまたは下向きに音波を飛ばし、反射波のドップラー効果からさまざまな深さでの流速分布を測定する機器です。その中から、一部のデータを見てみましょう。
使うデータは2018年7月27日~8月14日の1時間ごとのデータで、16m深の東西流速u（東向きが正）および南北流速v（北向きが正）です。
https://yasushifujiwara.github.io/dataanalysistutorial/data_KEO/KEO_2018typhoon.npz
import numpy as np
import matplotlib.pyplot as plt
# .npzファイルはnumpy配列をバイナリ形式で複数保存する形式です。
# np.savezで保存できます
data = np.load(r"/Users/yfuji/Desktop/lesson/data/data_KEO/KEO_2018typhoon.npz")
time = data["time"]
u = data["u"]
v = data["v"]
plt.figure(figsize=(14,4))
plt.plot(time, u, label="u")
plt.plot(time, v, label="v")
plt.legend()
plt.show()
この中で、7/28と8/7に台風が通過しています。台風の通過後、流速がある周期で振動しているのが見て取れると思います。
その周期は（この図からは見づらいですが）1日程度で、位相関係を考えるとuがvの1/4周期遅れ（uがsin, vがcos）なので、流れの方向が時計回りに回転していることがわかります。
これは慣性振動と呼ばれる現象で、低気圧の通過などによって広範囲で一斉に水平の流れが起こった時に、コリオリ力の作用で流れが曲げられて円運動をするために起こります。
今回は、この現象をシミュレーションしてみましょう。
現象のシミュレーションを行う（コンピュータ内で模擬する）ためには、2つの段階を踏む必要があります。
- 数理モデル化：関連する量を変数として選び、それらの関係を数式として表すことで、注目する現象を数式として表現する。
- 数値モデル化：コンピュータは有限桁の数を有限回の加減乗除算で扱うことしかできないため、1.の方程式をそうした形式で近似するアルゴリズムを構築し、近似解を求める。
ここでは、2.の段階でどのようなことを行っているのか、近似によってどのような誤差が生じるのかを学びます。
そのためには、1.の結果として出てくる方程式の正しいふるまいを知ることが重要です。
まず次節では1.の段階を踏み、解の「あるべき姿」を確認します。
方程式（モデル化）と真の解¶
海水の運動方程式は複雑ですが、流れは水平方向だけで、一様（空間的に一定）に起こり、圧力勾配力や粘性の影響も小さいと仮定して無視することで、次の2本の式にまで単純化できます。
$$ \frac{du}{dt} = fv\\ \frac{dv}{dt} = -fu $$ここでコリオリパラメータ$f=2\Omega \sin\phi$（$\Omega=2\pi/86164$ [rad/s]は地球の自転角速度、$\phi$は緯度）はコリオリ力の大きさに関わる係数です。
元々の運動方程式は$u(x,y,z,t)$や$v(x,y,z,t)$に関する、空間微分を含む偏微分方程式でしたが、ここでは流れが一様と仮定することで$u(t)$と$v(t)$に関する非常に単純な常微分方程式になっています。この一般解は次のようになります。
$$
u(t) = C_1 \cos ft + C_2 \sin ft\\
v(t) = -C_1 \sin ft + C_2 \cos ft
$$
時刻$t=0$で初期条件$u=u_0, v=v_0$ならば、$C_1=u_0, C_2=v_0$となるため次のようになります。
$$
u(t) = u_0 \cos ft + v_0 \sin ft\\
v(t) = -u_0 \sin ft + v_0 \cos ft
$$
この解の様子を、グラフを描いて見ておきましょう。
latitude = 32.3 # KEOの緯度
omega = 2*np.pi / 86164 # 地球の自転角速度
f = 2*omega * np.sin(np.pi / 180 * latitude) # コリオリパラメタの計算。度→ラジアンの変換に注意
time_sec = np.linspace(0, 5*24*60*60., 241) # 5日間まで、1800秒刻み
time_day = time_sec / 86400 # 日単位の時間、図示用
u0 = 1. # 初期条件、何でもいい
v0 = 0.
u = u0 * np.cos(f*time_sec) + v0 * np.sin(f*time_sec)
v = -u0 * np.sin(f*time_sec) + v0 * np.cos(f*time_sec)
plt.figure(figsize=(14,4))
plt.plot(time_day, u, label="u")
plt.plot(time_day, v, label="v")
plt.legend()
plt.show()
観測された流速変動と同じ位相差の振動が生じていることがわかります。
この現象でもう一つ重要な性質は、エネルギーが保存することです。$du/dt=fv$の両辺に$u$, $dv/dt=-fu$の両辺に$v$を掛けて辺々足すと、
$\frac{d}{dt}\left[\frac{1}{2}(u^2+v^2)\right]=0$、すなわち$\frac{1}{2}(u^2+v^2)=\textrm{const.}$が得られます。慣性振動する水の質量をかけると運動エネルギーになります。
現実に観測される慣性振動は、乱流が摩擦として働いたり、波（内部重力波）としてエネルギーが遠方に運び去られることで次第に振幅が減衰していきますが、ここではそれらの過程を無視して単純化しているため、エネルギーは一定になります。
近似方程式を解く¶
常微分方程式から漸化式へ¶
上で言及した通り、コンピュータは原則有限桁の数値しか扱えず、有限回の加減乗除算しかできません。そのため、極限操作を含む微積分は行えません。 そこで、離散化が必要になります。連続的な関数を扱う微分方程式を、飛び飛びの点での値を求める漸化式に変形させる、数値モデル化の中心的な行程です。
例えば微分であれば無限に小さい幅での変化を表す代わりに有限だが十分小さい幅での変化で近似することになります。これによって微分は単なる引き算、つまり差分で近似表現されます。
$$
\frac{du}{dt}=\lim_{\Delta t\rightarrow 0}\frac{u(t+\Delta t)-u(t)}{\Delta t}\approx\frac{u(t+\Delta t)-u(t)}{\Delta t}\\
\frac{dv}{dt}=\lim_{\Delta t\rightarrow 0}\frac{v(t+\Delta t)-v(t)}{\Delta t}\approx\frac{v(t+\Delta t)-v(t)}{\Delta t}
$$
ここで最右辺の$\Delta t$はある小さい幅で、状況に合わせて適切に選びます（その指針は今回の後半で扱います）。
これによって、今回の方程式は次のように表されます。
$$
u(t+\Delta t)=u(t) + (f\Delta t) v(t)\\
v(t+\Delta t)=v(t) - (f\Delta t) u(t)
$$
すると、開始時刻から$\Delta t$ごとの$u$, $v$が飛び飛びに得られるでしょう。
開始時刻を$t=0$とし、$t_n = n\Delta t$と表します。$t=t_n$での$u(t_n), v(t_n)$を$u_n, v_n$と表すと、上の式は次のようになります。
$$
u_{n+1}=u_n + (f\Delta t) v_n\\
v_{n+1}=v_n - (f\Delta t) u_n
$$
これは、漸化式ですね。一般項を求めなくても、代入法で時刻$n$の状態がわかれば時刻$n+1$の状態を求めることができます。その結果、数列$(u_1, u_2, u_3, ...), (v_1, v_2, v_3, ...)$が好きなだけ求められます。
上の漸化式をいきなり解く前に、より簡単な漸化式を代入法で求めていくプログラムを作ってみましょう。
練習5-1¶
次の漸化式を満たす数列の、最初の21項（$a_{20}$まで）を求めなさい。
- $a_0 = 0,\ a_{n+1}=2a_n+1$
- $a_0 = 0, a_1 = 1,\ a_{n+2}=a_{n+1}+a_{n}$（フィボナッチ数列、ヒント：前2項を記憶しておく必要がある）
実際に方程式を解く¶
それでは、実際に連立漸化式$u_{n+1}=u_{n}+f\Delta t v_n$, $v_{n+1}=v_n-f\Delta t u_n$を解いていきましょう。
実装の前に、どの時刻まで計算するのか・$\Delta t$としてどんな値を選ぶのか・結果をどれだけ記録するのかというのかということを決める必要があります。
計算期間（積分期間とも言います）は必要に応じ決めればよいのですが、ここでは5日間とすることにしましょう。
連立漸化式に含まれる定数のうち、$f$は上で「真の解」を求めた際に計算された通りですが、$\Delta t$の値は選択の余地があります。
これは非常に重要な問題で、後程じっくり分析しますが、ここではひとまず10分=600秒として進めてみましょう。
1ステップが10分なので、5日間の計算には 5日/10分=720ステップが必要となります。今回の慣性振動の式のように隣接2項式の場合は、最低でも直前の結果（第n+1項を求めるときは第n項）を記憶しておけば漸化式を代入法で解くことができます。
実際にはu, vの時間変化の様子も残したいので、データを保存していく必要があります。
今回は毎ステップのu, vの値を配列の中に記録してゆきます。
# 下準備
latitude = 32.3 # KEOの緯度
omega = 2*np.pi / 86164 # 地球の自転角速度
f = 2*omega * np.sin(np.pi / 180 * latitude) # コリオリパラメタの計算。度→ラジアンの変換に注意
num_step = 720 # 全ステップ数
dt = 600 # delta t (時間間隔)、単位は[s]
time_sec = dt * np.arange(0, num_step) # 0, 1*600, 2*600, ..., 719*600秒（5日）
time_day = time_sec / 86400 # 日単位の時間、図示用
u0 = 1. # uの初期値, [m/s]
v0 = 0. # vの初期値, [m/s]
u = np.zeros(num_step) # uの記録、入れ物を先に作っておく
v = np.zeros(num_step) # vの記録、入れ物を先に作っておく
u[0] = u0 # uの初期値
v[0] = v0 # vの初期値
# 真の解も計算しておく
u_true = u0 * np.cos(f*time_sec) + v0 * np.sin(f*time_sec)
v_true = -u0 * np.sin(f*time_sec) + v0 * np.cos(f*time_sec)
# 実際にforループで計算してゆく
for n in range(1, num_step):
u[n] = u[n-1] + f * dt * v[n-1]
v[n] = v[n-1] - f * dt * u[n-1]
これで、5日間の流速変化が求まったはずです。以下で、結果を見てみましょう。
plt.figure(figsize=(14,4))
plt.plot(time_day, u, "C0", label="u")
plt.plot(time_day, v, "C1", label="v")
plt.legend()
plt.show()
この結果を見て、どのように思いますか。慣性振動をシミュレーションできているのでしょうか、それともできていないと捉えられるでしょうか。 今回は「真の解」がわかっているので、何が良くて何が悪いかを量的に議論できます。
練習5-2¶
上記の時系列グラフに、先に計算してある真の解u_true, v_trueを描き加えなさい。
u_trueは線色"C0", 線種"--", v_trueは線色"C1", 線種"--"として描きなさい。
エネルギーの時系列¶
上の真の解との比較により、「振幅が増加してしまう問題がある」という結論が得られたのではと思います。
このことを別の観点で確かめましょう。この常微分方程式は、エネルギー$\frac{1}{2}(u^2+v^2)$が保存される性質を持つのでした。
エネルギーを計算し、時系列を描いてみましょう。
KE = 0.5 * (u**2 + v**2)
KE_true = 0.5 * (u_true**2 + v_true**2)
plt.figure(figsize=(14,4))
plt.plot(time_day, KE, "C0", label="KE")
plt.plot(time_day, KE_true, "C1", label="KE_true")
plt.legend()
plt.show()
本当は一定であるべきなのに、エネルギーが増加してしまっている（しかも加速的に）ということがわかります。
この誤差をどう評価し、誤差が計算法によってどう変化するのかは後に議論しますが、ここでは流速の変化を描画する別の方法を試してみましょう。
ホドグラフ¶
これまでは横軸に時間、縦軸に時間変化する量を描いてきました。
これに対して、水平速度ベクトルを一気に描きたいときには横軸に東西流速$u$, 縦軸に南北流速$v$を描くことがあります。これをホドグラフと呼びます。
このようにすることで、上から見下ろしたときの速度ベクトルの先端が得られます。
一つのベクトルは1つの点になりますが、たくさんのベクトルの先端を繋ぐと線を描けます。よく描かれる図として流れの鉛直分布$(u(z), v(z))$や流れの時間変化$(u(t), v(t))$があります。
ホドグラフを描くときには、以下のようなコツがあります。
- どの点がどの時刻かわからないので、ある決まった時刻の点に印をつける（
"o"を使う） - 図の縦横比を1:1にする（
plt.gca().aspect("equal")）
plt.figure()
plt.plot(u, v, "-", color="C0", label="simulation")
plt.plot(u[::144], v[::144], "o", color="C0") # 24時間ごとに点を描く
plt.plot(u_true, v_true, "-", color="C1", label="true")
plt.plot(u_true[::144], v_true[::144], "o", color="C1") # 24時間ごとに点を描く
plt.gca().set_aspect('equal')
plt.legend()
plt.xlabel("u [m/s]")
plt.ylabel("v [m/s]")
plt.show()
真の解が同じ半径の円周上を回り続けているのに対し、計算結果は螺旋を描いて外に広がっているのがわかります。
計算誤差¶
上で見た通り、今回行ったシミュレーションは振動自体はそれらしく再現するものの、エネルギーが増大してしまうという問題があることがわかりました。
これは、本来極限を用いて定義される微分を有限の$\Delta t$で近似したことにより生じる誤差です。
コンピュータが有限桁の数値しか扱えず有限回の計算しか行えないことから、こうした誤差は避けられないと言えるでしょう。大事なのは結果がそうした誤差を含んでいることを認識したうえで、誤差に左右されない核心的結論を見出すことです。
ここでは計算誤差をどう評価するか・誤差はどうしたら減らせるのかを学びます。
異なる時間刻みでの実験¶
本来$\Delta t\rightarrow 0$であるべきところを有限の$\Delta t$で近似しているので誤差が生じているのであれば、$\Delta t$を小さくすることで誤差が小さくなる、と期待できます。また、逆に$\Delta t$を大きくしていくとどんなことが起こるのでしょうか。
ここでは、様々な$\Delta t$を用いて同様の計算を行ってみましょう。
誤差の評価¶
比較にあたり、誤差の大きさをどのように評価するか定義しておく必要があります。
今回は真の解を知っているので、ある決まった時刻（例えば1日後）での計算値$(u,v)$と真値$(u_\textrm{true},v_\textrm{true})$の違いを誤差と定義しましょう。
ただし、単なる差では誤差の大小が初期条件によって決まる速度の大小に依存してしまうので、差を典型的な速度で割った相対誤差を用いることにします。
相対誤差の基準も様々に選べますが、ここでは初期流速$(u_0,v_0)$を用いることにします。
練習5-3¶
上で計算したu, u_true, v, v_trueを用いて、誤差の時系列（各時刻の誤差）を計算しなさい。
誤差の時間変化をグラフとして描きなさい。
関数を使おう¶
今回は様々な時間刻み幅で同じ計算を繰り返すので、関数を用いてコードを使いまわすのがコーディングを減らすことにつながります。
「1日のステップ数を引数として与え、1日分の計算を行い、全ステップのtime, u, vを記録した配列を返す」関数を作りましょう。
実験によってステップ数（配列の大きさ）が異なるので、timeもそれぞれについて生成します。
後に詳しく調べますが、今回の計算法は「前進オイラー法」というので、それにちなんだ関数名を付けます。
def calc_forward(step_oneday):
# step_onedayは1日（86400秒）のステップ数
dt = 86400 / step_oneday # delta t (時間間隔)、単位は[s]
time = dt * np.arange(0, step_oneday+1) # 時間の配列
# 今回は、step_oneday+1要素の配列にしてちょうど1日後の結果までを含むようにします
u = np.zeros(step_oneday+1) # uの記録、入れ物を先に作っておく
v = np.zeros(step_oneday+1) # vの記録、入れ物を先に作っておく
u[0] = 1.0 # uの初期値, [m/s]
v[0] = 0.0 # vの初期値, [m/s]
for n in range(1, step_oneday+1):
u[n] = u[n-1] + f * dt * v[n-1]
v[n] = v[n-1] - f * dt * u[n-1]
return time, u, v
先ほどは$\Delta t$として10分を用いました。今回は、5分・10分・20分の設定で実験を3通り行いましょう。
5分ならば1日は288ステップ、10分ならば144ステップ、20分ならば72ステップです。
第1回では説明しませんでしたが、関数が複数の値を返すこともでき、これを受け取るにはtime, u, v = calc_forward(144)のようにします。
細かく言えば、関数は3つの返り値をまとめたタプル(time,u,v)を返していて、その要素を個別に受け取っています。詳しくは"tuple unpacking"で調べてみてください。
time5, u5, v5 = calc_forward(288)
time10, u10, v10 = calc_forward(144)
time20, u20, v20 = calc_forward(72)
plt.figure(figsize=(8,4))
plt.plot(time5/3600, u5, "-", color="C0", label="u, 5min")
plt.plot(time5/3600, v5, "--", color="C0", label="v, 5min")
plt.plot(time10/3600, u10, "-", color="C1", label="u, 10min")
plt.plot(time10/3600, v10, "--", color="C1", label="v, 10min")
plt.plot(time20/3600, u20, "-", color="C2", label="u, 20min")
plt.plot(time20/3600, v20, "--", color="C2", label="v, 20min")
plt.legend(loc="upper center")
plt.xlabel("time [hour]")
plt.show()
見ての通り、振幅の増加は5分→10分→20分の順に大きくなっていることがわかります。予想通り、時間ステップを細かく刻むほど誤差は小さくなると言えそうです。
誤差の$\Delta t$依存性¶
それでは、具体的にそれぞれの誤差の大きさを求めましょう。
1日経過時点、つまり最終ステップでの相対誤差を求め、$\Delta t$に対してどう変化しているかを考えます。
# 誤差を求める
def error(u, v):
# まず真値を求める
t_1d = 24 * 60 * 60
utrue_1d = u0 * np.cos(f*t_1d) + v0 * np.sin(f*t_1d)
vtrue_1d = -u0 * np.sin(f*t_1d) + v0 * np.cos(f*t_1d)
# 配列u, vの最後の要素（1日後の値）を用いて相対誤差を計算する
err = np.sqrt((u[-1]-utrue_1d)**2 + (v[-1]-vtrue_1d)**2) / np.sqrt(u0**2 + v0**2)
return err
err5 = error(u5, v5)
err10 = error(u10, v10)
err20 = error(u20, v20)
print(" 5min: ", err5)
print("10min: ", err10)
print("20min: ", err20)
$\Delta t$が5分のとき1日後の誤差は8.2%, 10分のとき17.0%, 20分のとき36.8%となっています。誤差はおおよそ$\Delta t$に比例していると言えるでしょう。そのため、例えば1日後の誤差を1%程度に抑えたければ時間ステップを5分の1/8、40秒程度にすればよいと推測できます。
このことから、この数値モデルは時間について1次の精度をもつとか1次精度といいます。もし誤差が$\Delta t^2$に比例していれば2次精度といい、3次精度以上も同様です。2次以上の精度がどのように実現されるかは、後に紹介します。
練習5-4¶
上で見た通り、この計算法は時間について1次精度で、$\Delta t$が5分のとき1日後の誤差が8%になります。この情報をもとに、次を推定しなさい。
- $\Delta t$が1分の場合の1日後の誤差
- $\Delta t$が5分の場合の5日後の誤差
(1)について、上で定義した関数calc_forwardを用いて確かめなさい。
数値スキーム¶
今回は微分方程式$\frac{du}{dt}=fv$を$\frac{u_{n+1}-u_{n}}{\Delta t}=fv_n$（もう一方の式についても同様）と近似しました。
しかし、微分方程式を差分方程式で近似する方法（数値スキームと言います）はこれだけではありません。
数値スキームにはそれぞれ誤差の大きさや計算量の多さといった特徴があり、それを知ることで$\Delta t$の選択など計算の設定を適切に行うことができます。
ここでは、前進オイラー法を含めたいくつかの数値スキームを紹介します。
前進オイラー法¶
前進オイラー法では、時間微分項を$(u_{n+1}-u_n)/\Delta t$と近似し、それ以外の項を第nステップの値で表現します。その結果、次のような式になります。 $$ \frac{u_{n+1}-u_{n}}{\Delta t}=fv_n,\\ \frac{v_{n+1}-v_{n}}{\Delta t}=-fu_n $$ この式を整理すると、次のようになります。以下、簡単のために$\gamma\equiv f\Delta t$と表します。 $$ u_{n+1}=u_n + \gamma v_{n},\\ v_{n+1}=v_n - \gamma u_{n} $$ 前進オイラー法は$u_{n+1}=$（nステップ目の量）と表せるので、単純な代入法で漸化式計算を進めていけるという利点があります。この特徴を持つスキームを陽解法と言います。
これまで調べた通り、振動のエネルギーが増幅する誤差があり、その誤差は$\Delta t$に比例して減少する、という特徴があります。
後退オイラー法¶
後退オイラー法では、時間微分項を$(u_{n+1}-u_n)/\Delta t$と近似し、それ以外の項を第n+1ステップの値で表現します。その結果、次のような式になります。
$$
\frac{u_{n+1}-u_{n}}{\Delta t}=fv_{n+1},\\
\frac{v_{n+1}-v_{n}}{\Delta t}=-fu_{n+1}
$$
この式を整理すると、次のようになります。
$$
u_{n+1} - \gamma v_{n+1}=u_n,\\
\gamma u_{n+1}+v_{n+1}=v_n
$$
nステップ目の値からn+1ステップ目の値を求める際に、単なる代入法では求められず、連立方程式を解く必要があることがわかります。こうしたスキームを陰解法と言います。
一般的にはこの連立方程式は簡単には解けず、反復法などの方法で手間をかけて解くことになります。しかし、今回に限っては両辺に次の行列を左からかけることで、代入法の形式に変形することができます。
$$
\frac{1}{1+\gamma^2}
\begin{pmatrix}
1 & \gamma \\ -\gamma & 1
\end{pmatrix}
$$
変形の結果、次の式が得られます。
$$
u_{n+1}=\frac{1}{1+\gamma^2}(u_n+\gamma v_n)\\
v_{n+1}=\frac{1}{1+\gamma^2}(-\gamma u_n+v_n)
$$
これを用いて、上と同様に5分・10分・20分の$\Delta t$で計算を行ってみましょう。
# 1日分の計算を行う関数を定義する
def calc_backward(step_oneday):
# step_onedayは1日（86400秒）のステップ数
dt = 86400 / step_oneday # delta t (時間間隔)、単位は[s]
gamma = f * dt
time = dt * np.arange(0, step_oneday+1) # 時間の配列
# 今回は、step_oneday+1要素の配列にしてちょうど1日後の結果までを含むようにします
u = np.zeros(step_oneday+1) # uの記録、入れ物を先に作っておく
v = np.zeros(step_oneday+1) # vの記録、入れ物を先に作っておく
u[0] = 1.0 # uの初期値, [m/s]
v[0] = 0.0 # vの初期値, [m/s]
for n in range(1, step_oneday+1):
u[n] = (u[n-1] + gamma * v[n-1]) / (1 + gamma**2) # gammaの定義と、
v[n] = (v[n-1] - gamma * u[n-1]) / (1 + gamma**2) # この2行だけ改変しました
return time, u, v
# 計算を実行する
time5, u5, v5 = calc_backward(288)
time10, u10, v10 = calc_backward(144)
time20, u20, v20 = calc_backward(72)
# 図を描く
plt.figure(figsize=(8,4))
plt.plot(time5/3600, u5, "-", color="C0", label="u, 5min")
plt.plot(time5/3600, v5, "--", color="C0", label="v, 5min")
plt.plot(time10/3600, u10, "-", color="C1", label="u, 10min")
plt.plot(time10/3600, v10, "--", color="C1", label="v, 10min")
plt.plot(time20/3600, u20, "-", color="C2", label="u, 20min")
plt.plot(time20/3600, v20, "--", color="C2", label="v, 20min")
plt.legend(loc="upper center")
plt.xlabel("time [hour]")
plt.show()
# 誤差を求める
err5 = error(u5, v5)
err10 = error(u10, v10)
err20 = error(u20, v20)
print(" 5min: ", err5)
print("10min: ", err10)
print("20min: ", err20)
前進オイラー法は振幅が増大してしまう誤差を持っていましたが、後退オイラー法は振幅が減衰してしまう誤差を持つことがわかります。
また、誤差の大きさは$\Delta t$に比例していることから、後退オイラー法も1次精度であるとわかります。
台形法¶
前進オイラー法がエネルギー増大、後退オイラー法がエネルギー減衰の誤差を持つのなら、この2つを組み合わせればちょうど良くなるのでは、と予想できます。
時間微分項を$(u_{n+1}-u_n)/\Delta t$と近似し、それ以外の項を第nステップの値と第n+1ステップの値の平均で表現したスキームを台形法といいます。
$$
\frac{u_{n+1}-u_{n}}{\Delta t}=f\frac{v_n+v_{n+1}}{2},\\
\frac{v_{n+1}-v_{n}}{\Delta t}=-f\frac{u_n+u_{n+1}}{2}
$$
$\gamma\equiv f\Delta t$を用いてこれを整理すると、
$$
u_{n+1}-\frac{\gamma}{2}v_{n+1}=u_n+\frac{\gamma}{2}v_n,\\
\frac{\gamma}{2}u_{n+1}+v_{n+1}=-\frac{\gamma}{2}u_n+v_n,
$$
となりますので、これも陰解法の一種です。
後退オイラー法のときと同様に式変形を行うと、次の式が得られます。
$$
u_{n+1}=\frac{1}{1+\gamma^2/4}\left[\left(1-\frac{\gamma^2}{4}\right)u_n+\gamma v_n\right]\\
v_{n+1}=\frac{1}{1+\gamma^2/4}\left[-\gamma u_n+ \left(1-\frac{\gamma^2}{4}\right)v_n\right]
$$
練習5-5¶
上の台形法の式を用いて、5分・10分・20分の$\Delta t$について時系列のグラフを描き、誤差を評価しなさい。
関数名はcalc_trapzとしなさい。（台形法=trapezoidal method）
ルンゲ・クッタ法¶
台形法は2次精度で誤差が小さい一方、陰解法であるという困難がありました。
陽解法でありながら高次精度を達成する方法として、ルンゲ・クッタ法という方法がよく用いられるので紹介します。
ルンゲ・クッタ法では前進オイラー法のような陽解法を複数回繰り返して適切な重みをつけてそれらを平均し、次のステップの値を高次精度で求めてゆきます。
中点法¶
2次精度のルンゲ・クッタ法の一つとして、中点法は次のような式で計算が行われます。まず、通常の前進オイラー法で中間の時刻の速度 $$ u^* = u_n + (f\Delta t/2) v_n\\ v^* = v_n - (f\Delta t/2) u_n $$ を求め、その時刻で求めた時間変化率を用いてnステップからn+1ステップを積分します。 $$ u_{n+1} = u_n + (f\Delta t) v^*\\ v_{n+1} = v_n - (f\Delta t) u^* $$
中点法の誤差を上と同様に評価してみましょう。
# 1日分の計算を行う関数を定義する
def calc_midpoint(step_oneday):
# step_onedayは1日（86400秒）のステップ数
dt = 86400 / step_oneday # delta t (時間間隔)、単位は[s]
gamma = f * dt
time = dt * np.arange(0, step_oneday+1) # 時間の配列
# 今回は、step_oneday+1要素の配列にしてちょうど1日後の結果までを含むようにします
u = np.zeros(step_oneday+1) # uの記録、入れ物を先に作っておく
v = np.zeros(step_oneday+1) # vの記録、入れ物を先に作っておく
u[0] = 1.0 # uの初期値, [m/s]
v[0] = 0.0 # vの初期値, [m/s]
for n in range(1, step_oneday+1):
um = u[n-1] + 0.5*gamma * v[n-1]
vm = v[n-1] - 0.5*gamma * u[n-1]
u[n] = u[n-1] + gamma * vm
v[n] = v[n-1] - gamma * um
return time, u, v
# 計算を実行する
time5, u5, v5 = calc_midpoint(288)
time10, u10, v10 = calc_midpoint(144)
time20, u20, v20 = calc_midpoint(72)
# 図を描く
plt.figure(figsize=(8,4))
plt.plot(time5/3600, u5, "-", color="C0", label="u, 5min")
plt.plot(time5/3600, v5, "--", color="C0", label="v, 5min")
plt.plot(time10/3600, u10, "-", color="C1", label="u, 10min")
plt.plot(time10/3600, v10, "--", color="C1", label="v, 10min")
plt.plot(time20/3600, u20, "-", color="C2", label="u, 20min")
plt.plot(time20/3600, v20, "--", color="C2", label="v, 20min")
plt.legend(loc="upper center")
plt.xlabel("time [hour]")
plt.show()
# 誤差を求める
err5 = error(u5, v5)
err10 = error(u10, v10)
err20 = error(u20, v20)
print(" 5min: ", err5)
print("10min: ", err10)
print("20min: ", err20)
1ステップ進めるのに必要な計算量が前進オイラー法に比べて倍になっていますが、誤差が大幅に低減されてることがわかります。
$\Delta t$が20分でも誤差が1%未満ですが、前進オイラー法で同じだけの誤差を実現するには$\Delta t$が40秒である必要でしたので、ステップ数は40秒/20分=1/30, 計算量は1/15まで低減されていることがわかります。
4次精度ルンゲ・クッタ法¶
ルンゲ・クッタ法にはさらに高次精度の方法があり、特に4次精度の公式は常微分方程式の計算によく用いられます。ベクトル量$u$の時間発展が $$ \frac{du}{dt}= f(u, t) $$ と表されるとき、$u_n$から${u}_{n+1}$は次の手続きで求められます。 $$ \begin{align*} h^{(1)}&=f(u_n, t_n)\\ h^{(2)}&=f(u_n+h^{(1)}\Delta t/2, t_n+\Delta t/2)\\ h^{(3)}&=f(u_n+h^{(2)}\Delta t/2, t_n+\Delta t/2)\\ h^{(4)}&=f(u_n+h^{(3)}\Delta t, t_n+\Delta t)\\ u_{n+1}&=u_n+\frac{\Delta t}{6}(h^{(1)}+2h^{(2)}+2h^{(3)}+h^{(4)}) \end{align*} $$
高次精度解法の威力をお見せするため、3時間・1時間・20分の$\Delta t$での誤差を求めてみましょう。
# 1日分の計算を行う関数を定義する
def calc_rk4(step_oneday):
# step_onedayは1日（86400秒）のステップ数
dt = 86400 / step_oneday # delta t (時間間隔)、単位は[s]
time = dt * np.arange(0, step_oneday+1) # 時間の配列
# 今回は、step_oneday+1要素の配列にしてちょうど1日後の結果までを含むようにします
u = np.zeros(step_oneday+1) # uの記録、入れ物を先に作っておく
v = np.zeros(step_oneday+1) # vの記録、入れ物を先に作っておく
u[0] = 1.0 # uの初期値, [m/s]
v[0] = 0.0 # vの初期値, [m/s]
for n in range(1, step_oneday+1):
hu1 = f*v[n-1]
hv1 = -f*u[n-1]
hu2 = f*(v[n-1] + 0.5*dt*hv1)
hv2 = -f*(u[n-1] + 0.5*dt*hu1)
hu3 = f*(v[n-1] + 0.5*dt*hv2)
hv3 = -f*(u[n-1] + 0.5*dt*hu2)
hu4 = f*(v[n-1] + dt*hv3)
hv4 = -f*(u[n-1] + dt*hu3)
u[n] = u[n-1] + dt/6 * (hu1 + 2*hu2 + 2*hu3 + hu4)
v[n] = v[n-1] + dt/6 * (hv1 + 2*hv2 + 2*hv3 + hv4)
return time, u, v
# 計算を実行する
time20, u20, v20 = calc_rk4(72)
time60, u60, v60 = calc_rk4(24)
time180, u180, v180 = calc_rk4(8)
# 図を描く
plt.figure(figsize=(8,4))
plt.plot(time20/3600, u20, "-", color="C2", label="u, 20min")
plt.plot(time20/3600, v20, "--", color="C2", label="v, 20min")
plt.plot(time60/3600, u60, "-", color="C3", label="u, 1hour")
plt.plot(time60/3600, v60, "--", color="C3", label="v, 1hour")
plt.plot(time180/3600, u180, "-", color="C4", label="u, 3hour")
plt.plot(time180/3600, v180, "--", color="C4", label="v, 3hour")
plt.legend(loc="upper center")
plt.xlabel("time [hour]")
plt.show()
# 誤差を求める
err20 = error(u20, v20)
err60 = error(u60, v60)
err180 = error(u180, v180)
print("20min: ", err20)
print("1hour: ", err60)
print("3hour: ", err180)
$\Delta t$が20分の計算では誤差が0.0004%, 3時間でも2.7%ということで、非常にざっくりした計算でも高精度な結果が得られているとわかります。このように、計算法によって誤差と計算量が様々に変化するので、計算法の特性を知っておくことは有用なことです。
適切な数値シミュレーションのために (1)¶
ここまで紹介した常微分方程式の数値計算法の基本は、自分が数値計算プログラムを作成する場合はもちろん、既存のプログラムを利用するときにも重要な知識です。
今回の内容から見出せるいくつかの教訓をまとめておきましょう。
適切な数理モデリングの重要性¶
実際の海水の運動方程式は非常に複雑で、その代わりに様々な現象を解として含みます。
ここではADCPのデータをもとに水平一様性やコリオリ力以外の力を省略するなどの近似を施し、現象の本質を抜き出した非常に単純な常微分方程式に帰着した結果、簡単な式で現象の再現をできました。
このように、数理モデリングの段階で注目する現象に必要十分な要素を抜き出しておくことで、複雑な現象の本質をとらえた実験を最小限の労力で実施できます。これを行うには、現象の本質を見抜く観察眼が重要です。
数値スキームの選択¶
今回多くの数値スキームを紹介しました。実用においては、どのようにスキームを選ぶべきなのでしょうか。
1つの対比軸として、単純なスキーム（前進オイラー法など）と複雑なスキーム（ルンゲ・クッタ4次など）があります。 これまで見た通り、一般的に単純なスキームは計算量が少ないという利点がある一方、誤差が大きいという欠点があります。 逆に、複雑なスキームは誤差が小さく、しかも$\Delta t$を少し小さくしただけで誤差を大きく減らせる一方、計算量やメモリ使用量が多くなるという欠点があります。
また、別稿で解説する予定の安定性という観点も重要です。不安定な数値スキームはエネルギーが徐々に増大し、やがて数列が発散してしまいます。
興味のある現象ならば高次精度スキームを使ったり差分幅を小さくとることで増大する誤差を抑えればよいのですが、現象の再現精度が重要でない場合には安定なスキームを用いることが得策になってきます。
例えば、1000年分の海洋シミュレーションを行って数年規模で変動する風成循環の気候値を計算したい場合に慣性振動の再現性は重要ではありません。
その際、コリオリ力を前進オイラー法などの不安定なスキームで積分していては振動がどんどん増大し、やがて非現実的な流れになってしまうので、$\Delta t$を非常に細かくとる必要が出てくるでしょう。
興味のある海洋現象の精度向上にはほとんど寄与しないのに、いたずらに$\Delta t$を増やすことは非効率的です。
そこで、後退オイラー法のような安定なスキームを採用することで、大きな$\Delta t$でも慣性振動が増大せずに計算を続けられます。
循環の持続時間より十分短い$\Delta t$を取っていれば興味のある現象はきちんと表現できる一方、周期の短い慣性振動は勝手に減衰していってくれます。
差分幅の選択¶
差分幅$\Delta t$を小さくとれば誤差は小さくなりますが、その分同じ期間（例えば1日）を計算するのに必要な計算量が増えてしまいます。$\Delta t$はどの程度小さくとれば十分と言えるのでしょうか。
一般的な指針として、計算誤差が誤差許容範囲よりも十分小さく、シミュレーションから得られる結論が$\Delta t$に左右されない程度に小さい値を選択すべきです。
前進オイラー法での計算の際、5分刻みでは1日後に流速が約8%増大してしまうことがわかりました。
例えば5日後の流速を予想しようとしており、初期条件が20%程度の不確実性を含んでいるとわかっているならば、計算誤差の影響は少なくとも20%の数分の1, できれば10分の1以下に押さえたいところです。
もし前進オイラー法で5日後の誤差を2%にすることを目指すのならば、1日あたり$(1.02^{1/5}-1)\times 100\approx0.4\%$の誤差になるように、$\Delta t=(0.4/8)\times 5\textrm{min}=15\textrm{sec}$の刻み幅を選ぶのが目安となるでしょう。
もちろん高次精度スキームならより大きな$\Delta t$で間に合います。
今回はあえて「答えを知っている」問題を考え、誤差の$\Delta t$依存性を真の解との比較で調べました。しかし、実際の問題のほとんどでは真の解を得られません。
そんな時に有効な指針として、現在の実験結果を$\Delta t$を小さくとった計算と比較し、結論が変わらないことを確かめるということが挙げられます。
適切な離散化が行われた数値スキームは時間刻みを0に近づけることで元の微分方程式に収束する（一貫性という）ので、例えば$\Delta t$を半分に取った計算と結論が変わらなければ、これ以上$\Delta t$を小さくとっても結論は変わらないと考えられるでしょう。
とはいえ、すべての判断について$\Delta t$を半分にした計算で検証を行うことはできませんので、設定を新しく作った時などの折々に検証するのが良いでしょう。