モデリング
モデリングとは
まずはモデリングそのものを解説しよう. モデリングに関する、わかりやすい古典的な題材を資料 モデリング古典3題 (pdf)で紹介する.
単純な常微分方程式を Euler 法で計算してみよう
対象の微分方程式
対象として, 授業時に示した、人口 $u = u(t)$ ($t$: 時間) に対する常微分方程式 $$ \frac{du}{dt} = \gamma\, \left( \sin(u) \right)^{1.2}, \gamma > 0, t \in (0,\infty) $$ を考えよう. ただし、$u$ の初期値 $u(0)$ については $0 < u(0) < 1$ としておく.
時間刻み幅 $\Delta t$ を一定にして、微分を差分で近似することで近似計算式(スキームと呼ぶ)を作ることを考えよう. 以降、$u_n \cong u(n\Delta t)$ と想定した近似値 $u_n$ を用いて話を進めよう.
スキームの作り方は無数にあるが、ここでは最も計算やプログラミングが簡単な Euler 法 を使ってみよう.
これは、 微分方程式中の時間微分項 $du/dt$ を陽的差分 $( u_{n+1} - u_n ) / \Delta t$ で近似し、 他の部分に出てくる $u(t)$ はそのまま $u_n$ で近似することで、 既知量 $u_n$ を代入するだけでただちに未知量 $u_{n+1}$ が求まる ように近似式を構成する方法を言う.
実際、今回の対象としている微分方程式に Euler 法を適用してみると、
$$ \frac{u_{n+1} - u_n}{\Delta t} = \gamma\, \left( \sin(u_n)\right)^{1.2} $$
となる(これが Euler スキームと呼ばれる). 整理するとさらにわかりやすく
$$ u_{n+1} = u_n + \Delta t \gamma\, \left( \sin(u_n)\right)^{1.2} $$
となるので、実際はこの式を用いると良いだろう.
初期値は既知で $u_0 = u(0)$ とすることができるので、さらなる近似値を上の Euler スキームを用いて $u_1, u_2, \cdots$ と順番に求めていけば好きなだけ先の近似値を計算できることがわかる.
まず各種定数などを定義しよう
数字を与えておかないといけない変数を設定しておこう.
|
|
0.1
Euler スキームを関数として定義しよう
Julia では関数を定義する方法がいくつかある.
以前示したように、一番手堅い方法は function 命令を使う方法で、例えば以下のようになる.ちなみに、return 命令を使うと出力値を指定できる.return を使わない場合は、最後の式の結果が出力値になる.
|
|
Euler (generic function with 1 method)
単純な関数なら、数学の数式で書くような方法で次のようにしても定義できる.
|
|
Euler2 (generic function with 1 method)
写像の表現に似た無名関数という機能を使って次のようにしても定義できる.上手く動くなら好きな方法を選べば良い.
|
|
(::#1) (generic function with 1 method)
蛇足だが、どの定義も本当に同じか、気になる人へ
上のいずれも同じ関数を定義していることを、実際に値を入れて確認してみよう.
|
|
(1.0812918438981398, 1.0812918438981398, 1.0812918438981398)
それでも動作の違いが有るのではないか、と不安な人は @code_lowered マクロを用いて、次のようにそれぞれの内部コードが同一であることを確認しても良い.
|
|
CodeInfo(:(begin nothing return u + Main.Δt * Main.γ * (Main.sin)(u) ^ 1.2 end))
|
|
CodeInfo(:(begin nothing return u + Main.Δt * Main.γ * (Main.sin)(u) ^ 1.2 end))
|
|
CodeInfo(:(begin nothing return u + Main.Δt * Main.γ * (Main.sin)(u) ^ 1.2 end))
内部的に全く同等の動作をするはずだ、ということがこれを見ると分かるだろう.
定義した Euler スキームを使って、初期値から順番に近似値を計算しよう
単に順番に計算していけばいいだけだが、あとでグラフを描きたいので、結果を配列に記録しておこう.
|
|
結果のグラフを見てみよう
上の計算が終わると、配列 u_sq に $u_0, u_1, \cdots, n_{100}$ の要素が入っているはずである. そこでこの数列をプロットしてみよう.
いつものように、Plots パッケージとそのサブパッケージ GR を下記のようにして呼び出しておいて…
|
|
Plots.GRBackend()
|
|
tips: 上のグラフの横軸をベクトルの要素番号ではなく、きちんと時間 $t$ そのものにしたいという場合は、$x$軸用の変数も作った上で, $x$軸用の変数と $y$軸用の変数(今回は u_sq)を plot 関数に一緒に渡せば良い.具体的には以下のような感じだ.
|
|
ちなみに、上の xaxis! 等はオプションとして渡しても良い(! マークは外してね).次のような感じだ.
|
|
さらに、legend (図に書き込まれる「凡例」)が邪魔なら、legend オプションを適当にいじれば良い. 例えば次のようにすれば消せる.
|
|
$u$ が $\pi$ に近づいていくという話は既に授業で聞いたかもしれないね. そこで、グラフに $\pi$ の値を示す「線」を追加しておこう.次のような感じだ.
|
|
Euler法は失敗しやすい.見てみよう.
Euler 法は数式の変形も簡単だし、プログラムも簡単だ.実行速度も速い.しかしその分、結果の精度は悪いし、計算過程は不安定だ.ここでいう不安定というのは、数値解が発散したりすることを言う.
このことを、実際に見てみよう.やり方は簡単で、時間発展を追いかける際の $\Delta t$ を少し大きめにするだけでいい…
|
|
1.2
さて、さきほどと全く同じように計算してみよう.ただし、過程がわかりやすいように、途中結果をいちいち表示させよう.
|
|
1: 0.1: 0.09983341664682815 2: 0.17556355249608085: 0.17466305501555435 3: 0.3234139541098293: 0.3178053673300265 4: 0.6266457368366771: 0.5864311118638564 5: 1.2591182402790926: 0.9518203144684844 6: 2.390078148041807: 0.6827461247045289 7: 3.149166519606018: -0.007573793605922763 DomainError: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar. Stacktrace: [1] nan_dom_err at ./math.jl:300 [inlined] [2] ^(::Float64, ::Float64) at ./math.jl:699 [3] macro expansion at ./In[17]:8 [inlined] [4] anonymous at ./
:?
エラーを示して停まってしまう.何が起きたのか?
上のエラーメッセージを読むと、(負の数)^1.2 の部分で出ていて、そこで文句を言っている.
そしてさらに数字を丁寧に見よう.
すると、7回目の出力の $u$ が $\pi$ をわずかに超えている様子がみてとれる. しかしそもそも、初期値が小さな正の実数ならば,もとの微分方程式では解は $\pi$ 未満にとどまるためこれはあり得ないシチュエーションである.
どうやら、Euler 法の誤差によって近似解が $\pi$ 以上になってしまうという状況が発生してしまった様子だ.
そのあたりを、plot を使ってみてみよう.
|
|
ふ~む、先の結果とどう違うのか.時間軸を丁寧に揃えて、同時に表示してみよう. まず、近似解の値そのものを比べてみよう.
|
|
見ると分かるが、「グラフの $y1$, つまり粗い方」は解が大きくなっていくスピードにブレーキがかかるのが大変に遅れる、という違いになっている.そのため、突破するはずのない値(この問題の場合は $\pi$ のこと)を近似解が突破してしまう、という格好だ.
その結果、$\sin(u)$ がありえない値である「負」になってしまうのだ.それもグラフで比較してみよう.
|
|
うまく計算できている方はこの値が「正のまま」であるが、駄目な方は「負」になってしまっている.これを 1.2 乗しようとするもんだから、プログラムが動かなくなってしまっているのだ.
Report No.03 近似解がうまくない.自分でアイディアを.
次回以降にさらにいろいろな手法を紹介するけれども、まずは自力で少し考えてみよう(調べてもよいが、できるだけ自分で考えてみよう). $\Delta t$ を小さくするという方法以外にどんな解決方法がありうるだろうか. そして、実際にそれを実行できるだろうか.
馬鹿げた方法と思うようなものでも良いから、アイディアを出し、実際にその方法にチャレンジしてみて、その結果をレポートしよう.