05. モデリングと Euler法
Photo by Roman Mager on Unsplash
モデリングとは
モデリングとは,対象として考えている現象に対して数学などの信頼できる理論で説明して「理論仮設 = モデル」を作り上げることである.
つまり,「こうなっている」という仮設をたてることである.
通常は,その仮設を検証できるように数学を用いてモデリングする(微分方程式を用いることが多いよ).
なお,モデルは完全である必要はなく,あくまで「この現象のこの挙動をある程度説明できる」というもので良い1.
モデリングに関する、わかりやすい古典的な題材を資料 モデリング古典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 法 を使ってみよう.
この Euler 法とは時間発展型の微分方程式問題に対してスキームを作成する方法2で、対象とする微分方程式の解が $u(t)$ とすると、 時間刻み幅 $\Delta t$ を導入した上で、
- 近似解として数列 $\{u_n\}_n$ を $u_n \cong u(n\Delta t)$ を想定して導入する.
- 微分方程式中の時間微分項 $du/dt$ を陽的差分 $( u_{n+1} - u_n ) / \Delta t$ で近似し、 他の部分に出てくる $u(t)$ はそのまま $u_n$ で近似する.
としてスキームを作成するものだ. 一回やってみるとわかるが、こうして得られたスキームは、 既知量 $u_n$ を代入するだけでただちに未知量 $u_{n+1}$ が求まる ような便利な形になっている.
既知量 $u_n$ を代入するだけでただちに未知量 $u_{n+1}$ が求まる形になっているスキームを (陽的)Euler スキームとよぶ.
実際、今回の対象としている微分方程式に Euler 法を適用してみると、
$$ \frac{u_{n+1} - u_n}{\Delta t} = \gamma\, \left( \sin(u_n)\right)^{1.2} $$
となる. 左辺に新しい値による項を集め,右辺に古い値による項を集める形で整理するとさらにわかりやすく
$$ u_{n+1} = u_n + \Delta t \,\, \gamma \,\, \left( \sin(u_n)\right)^{1.2} \label{eq:euler-scheme} $$
となるので、実際はこの式を用いると良いだろう.
初期値は既知で $u_0 = u(0)$ とすることができるので、それ以降の近似値を上の Euler スキームを用いて $u_1, u_2, \cdots$ と順番に求めていけば好きなだけ先の近似値を計算できることがわかる.
まず各種定数などを定義しよう
数字を与えておかないといけない変数を設定しておこう.
|
|
Euler スキームを関数として定義しよう
Julia では関数を定義する方法がいくつかある.
以前示したように、一番手堅い方法は function 命令を使う方法で、例えば以下のようになる.ちなみに、return 命令を使うと出力値を指定できる.return を使わない場合は、最後の式の結果が出力値になる.
|
|
単純な関数なら、数学の数式で書くような方法で次のようにしても定義できる.
|
|
写像の表現に似た無名関数という機能を使って,その無名関数に名前をつける形で次のようにしても定義できる.
|
|
上手く動くならどの方法でも良い.だから好きな方法を選べば良いよ.
蛇足だが、どの定義も本当に同じか、気になる人へ
上のいずれも同じ関数を定義していることを、実際に値を入れて確認してみよう.
|
|
(1.0812918438981398, 1.0812918438981398, 1.0812918438981398)
それでも動作の違いが有るのではないか、と不安な人は @code_lowered マクロを用いて、次のようにそれぞれの内部コードが同一であることを確認しても良い.
|
|
CodeInfo(
1 ─ %1 = Main.sin(u)
│ %2 = %1 ^ 1.2
│ %3 = Main.Δt * Main.γ * %2
│ %4 = u + %3
└── return %4
)
|
|
CodeInfo(
1 ─ %1 = Main.sin(u)
│ %2 = %1 ^ 1.2
│ %3 = Main.Δt * Main.γ * %2
│ %4 = u + %3
└── return %4
)
|
|
CodeInfo(
1 ─ %1 = Main.sin(u)
│ %2 = %1 ^ 1.2
│ %3 = Main.Δt * Main.γ * %2
│ %4 = u + %3
└── return %4
)
みると、示されている内容がなんとなくわかるだろう. そして、内部的に全く同等の動作をするはずだ、ということがこれを見ると分かるだろう.
定義した Euler スキームを使って、初期値から順番に近似値を計算しよう
単に順番に計算していけばいいだけだが、あとでグラフを描きたいので、結果を配列に記録しておこう.
|
|
結果のグラフを見てみよう
上の計算が終わると、配列 u_sq に $u_0, u_1, \cdots, n_{100}$ の要素が入っているはずである. そこでこの数列をプロットしてみよう.
いつものように、Plots パッケージを下記のようにして呼び出しておいて…
|
|
データを貯めてある配列を直接 plot しよう.
|
|
上のグラフの横軸をベクトルの要素番号ではなく、きちんと時間 $t$ そのものにしたいという場合は、$x$軸用の変数も作った上で, $x$軸用の変数と $y$軸用の変数(今回は u_sq)を plot 関数に一緒に渡せば良い.具体的には以下のような感じだ.
|
|
(前にも書いたが) 命令の最後に ! がついているものは、「追加で」という意味を示す.
ちなみに plot 命令のオプションとして xaxis! 等を渡しても良い( "!" は外してね).次のような感じだ.
|
|
さらに、legend (図に書き込まれる「凡例」)が邪魔なら、legend オプションを適当にいじれば良い. 例えば次のようにすれば消せる.
|
|
さて,この場合グラフからも想像できるように,$u$ が $\pi$ に近づいていく(数式をみて納得しておこう).
そこで、グラフに $\pi$ の値を示す「線」を追加しておこう.
定数関数を使ってもよいが以前も出てきた hline
という命令でやってみよう.次のような感じだ.
|
|
Euler法は失敗しやすい.見てみよう.
Euler 法は数式の変形も簡単だし、プログラムも簡単だ.実行速度も速い.しかしその分、結果の精度は悪いし、計算過程は不安定だ.ここでいう不安定というのは、数値解が発散したりすることを言う.
このことを、実際に見てみよう.やり方は簡単で、時間発展を追いかける際の $\Delta t$ を少し大きめにするだけでいい…
|
|
さて、さきほどと全く同じように計算してみよう.ただし、過程がわかりやすいように、途中結果をいちいち表示させよう.
|
|
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 with -0.007573793605922763:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
[1] throw_exp_domainerror(::Float64) at .\math.jl:37
[2] ^ at .\math.jl:888 [inlined]
[3] Euler(::Float64) at .\In[2]:5
[4] top-level scope at .\In[18]:9
[5] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
[6] execute_code(::String, ::String) at c:\julia\PKG\v1.5\packages\IJulia\a1SNk\src\execute_request.jl:27
[7] execute_request(::ZMQ.Socket, ::IJulia.Msg) at c:\julia\PKG\v1.5\packages\IJulia\a1SNk\src\execute_request.jl:86
[8] #invokelatest#1 at .\essentials.jl:710 [inlined]
[9] invokelatest at .\essentials.jl:709 [inlined]
[10] eventloop(::ZMQ.Socket) at c:\julia\PKG\v1.5\packages\IJulia\a1SNk\src\eventloop.jl:8
[11] (::IJulia.var"#15#18")() at .\task.jl:356
エラーを示して停まってしまう.何が起きた?
エラーが起きて停まった. プログラムを書いて動かす場面ではエラーはよく出てくるものであるから,これと付き合えないと話が進まない.
今回はせっかくなので,エラーメッセージの読み方を知って,状況を把握できるようになっておこう.
最初に,上のエラーメッセージの stacktrace という部分を見よう. ここには,エラーが起きた箇所から始まって,その動作を呼出した場所へと逆上っての情報が記されている.
真面目にいくならば,[1] から順に読んでいって,まずは top-level scope という文字が出てくる箇所を探そう. top-level scope というのはプログラムを動かしている一番表面ということで, この場合は jupyter に書き込んだプログラムそのものなのだ.
実際,この場合,
[4] top-level scope at .\In[18]:9
というものが見つかる.
そしてこの at
以降を読む.
これが,このプログラムのエラー箇所の「一番表面的な箇所」になる.
この場合 In[18]:9
となっていて, In[18]
は jupyter の入力の In [18] という箇所3だ,ということを意味し,そして最後の :9
はその 9行目だ,ということを表している.
実際,上の入力の 9行目をみると u = Euler(u)
となっているので,自分で定義した関数 Euler
でエラーが出た,ということがわかる.
次に,この場所から stacktrace を上へ向かって読み進めて状況を細かく把握していく. この場合,
[3] Euler(::Float64) at .\In[2]:5
となっていて,関数 Euler
の 5行目,つまり return u + Δt * γ * (sin(u))^1.2
でエラーになったということがわかる.
次にさらに stacktrace を上へ向かい,
[2] ^ at .\math.jl:888 [inlined]
を見ることでエラーの起きた行の ^
という箇所,つまり,sin(u)^1.2
の 1.2乗の計算でエラーになったということが判明する.
さらに stacktrace を上に行くと
[1] throw_exp_domainerror(::Float64) at .\math.jl:37
とあり,「exp_domainerror4 というエラー処理をした」と言っている.
そしてこのエラーについてはこのエラーメッセージ中に解説が書いてあるのでそれを読めば良い(stacktrace の追いかけはここで終了だ).
そしてそれは
DomainError with -0.007573793605922763:
で始まる 3行のエラー説明部分である.
読んでみると,
Exponentiation yielding a complex result requires a complex argument.
意訳: 結果が複素数になる指数計算は引数が複素数であるべき
と言っている(さらにちょっとした例も載っているね).
さて,stacktrace の追いかけとこのメッセージから考えると,
sin(u)^1.2
の
^1.2
の結果が複素数になってしまった,だからエラーだ,ということだ.
となるとその引数の
sin(u)
が怪しい.
具体的には, sin(u) = -0.007573793605922763
となり,負の値になってしまったのがエラーの原因では? と怪しまれる5.
この数字はプログラムが出力した直前の表示
7: 3.149166519606018: -0.007573793605922763
の右端の数字であることから、まさにこれは sin(u)
の結果であることが確認できる.
よって、エラーは
|
|
を計算しようとして出たものと考えてよいはずだ.
確認のため,実際にこれを実行させてみよう.そうすると、
DomainError with -0.007573793605922763:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
[1] throw_exp_domainerror(::Float64) at .\math.jl:37
[2] ^(::Float64, ::Float64) at .\math.jl:888
[3] top-level scope at In[22]:1
…略…
と出力される.これは先のものと全く同じエラーなので,この考えが正しかったことが確かめられる.
さて、エラーの原因がはっきりと見えたので、何が起きたのか,数字を丁寧に追いかけてみよう.
すると、7回目の出力の $u$ が $\pi$ をわずかに超えている様子がみてとれる. しかしそもそも、初期値が小さな正の実数ならば,もとの微分方程式では解は $\pi$ 未満にとどまるためこれはあり得ないシチュエーションである.
どうやら、Euler 法の誤差によって近似解が $\pi$ 以上になってしまうという状況が発生し,そのために sin(u)
が負になってしまった様子だ.
そのあたりを、plot を使って確かめていこう.
まず,$\Delta t$ が小さいときの結果とどう違うのか.時間軸を丁寧に揃えて、同時に表示してみよう. まず、近似解の値そのものを比べてみよう.
|
|
見ると分かるが、「失敗したほう, つまり粗い方」は解が大きくなっていくスピードにブレーキがかかるのが大変に遅れる、という違いになっている.そのため、突破するはずのない値(この問題の場合は $\pi$ のこと)を近似解が突破してしまう、という格好だ.
その結果、$\sin(u)$ がありえない値である「負」になってしまうのだ.それもグラフで比較してみよう.
|
|
うまく計算できている方はこの値が「正のまま」であるが、駄目な方は「負」になってしまっている.これを 1.2 乗しようとするもんだから、プログラムが動かなくなってしまっている,ということがこの図でも確認でき,これでエラーの正体がきっちり把握できた,ということになるだろう.
じゃあ Δt はいくつなら動くのか?
このエラーの経験から,「ならば Δt はいくつにするべきなのか?」という疑問が生ずると思う. 実はこれは「数値安定性を確保するためのパラメータの決定」という問題6で,一般には結構難しい.
ただまあ今回はスキームが大変シンプルなので,実は結構簡単に調べがつく. やってみよう.
まず,Eulerスキームの計算が直前までうまくいっている(つまり,$u_n < \pi$ ということ)状況下で今回のエラーを回避するための条件
\[ u_{n+1} < \pi \]
について今回の Euler スキーム $\eqref{eq:euler-scheme}$ を眺めて考える. するとそれは
\[ u_n + \Delta t \,\, \gamma \,\, \left( \sin(u_n)\right)^{1.2} < \pi \]
ということであるので,これを $\Delta t$ についての式に変形すると
\[ \Delta t < \frac{\pi - u_n}{\gamma\,\, \sin(u_n)^{1.2}} \]
となる. よって,この右辺の最小値を調べれば良い.
そこで \[ Up(u) = \frac{\pi - u}{\gamma\,\,\sin(u)^{1.2}} \label{up:up-u} \] と定義して,まずはこのグラフを描いてみよう.
|
|
グラフは左端が一番大きく,よーく見ると右端も少し上がっている. そこでもう少しわかりやすいように範囲を狭めてグラフを描いてみよう.
|
|
このグラフに依ると $u$ が 2.4 の付近? に最小値があるように見受けられる. 実際,真面目に最小値を満たす $u$ を求めると
\[ u_{\mbox{min}} = 2.446142705263389 \label{eq:u-min} \]
と求められ,そしてそのときの $Up(u)$ の値は
\[ Up(u_{\mbox{min}}) = 1.1864652653081016 \]
となる.つまり,
\[ \Delta t < 1.1864652653081016 \]
ならば今回のエラーを起こさずに計算が可能だ,ということになる.
実際,Δt = 1.186465
として計算しても確かに計算がうまくいくことなどが実際に確認できる.やってみよう.
レポート
下記要領でレポートを出してみよう.
- e-mail にて,
- 題名を 2020-numerical-analysis-report-05 として,
- 教官宛(アドレスは web の "TOP" を見よう)に,
- 自分の学籍番号と名前を必ず書き込んで,
- 内容はテキストでも良いし,pdf などの電子ファイルを添付しても良いので,
下記の内容を実行して,結果や解析,感想等をレポートとして提出しよう.
- 上の $(\ref{eq:u-min})$ 式の $u_{\mbox{min}}$ の値を自分でも求めてみよう.なお,近似精度はあまり気にしなくて良い.
やり方はいろいろある. 一番オーソドックスなのは式 $(\ref{up:up-u})$ を $u$ で微分した式を手で求めておいて,それがゼロになるところを求める方法だろう.これなら前回までのゼロ点求解の知識が使える.
他には,たとえばx in 2.3:0.01:2.5
の範囲で数値的に最小値を探す,なんて強引な方法もあるだろう. いろいろ工夫してみるとよいだろう. - $\Delta t$ を大きくしたい場合は?
今回の Euler スキームの場合は $\Delta t$ が大きいと近似解がうまく得られないことがわかった. では、$\Delta t$ を大きくとりたいときはどうしたらよいだろうか.
次回以降にさらにいろいろな手法を紹介するけれども、まずは自力で少し考えてみよう(調べてもよいが、できるだけ自分で考えてみよう). $\Delta t$ を小さくするという問題回避方法以外にどんな解決方法がありうるだろうか. そして、実際にそれを実行してうまくいくか確認してみよう.
馬鹿げた方法と思うようなものでも良いから、アイディアを出し、実際にその方法にチャレンジしてみて、その結果をレポートしよう. - 似て非なる、解きにくい問題に取り組んでみる.
\[
\frac{d u}{d t} = - \sqrt{ u }
\]
という常微分方程式に基づく時間発展問題を、たとえば
時間 $t=0$ での
初期値
$u_0 = 0.25$
のもとで考えよう.
このとき、Euler スキームを構成し、$\Delta t$ をいろいろ変えて計算して結果を見比べることで以下の問にチャレンジしよう.
- 厳密解があるとして、どんな形であると推測されるか.推測してみよう.
- Euler 法では実質的に $t \in [0, T)$ で近似解が求めることができるとすると、$T$ はいくつぐらいか.そしてそれはどうしてその値になるのか
- (チャレンジ問題) 厳密解を式で書けるか. 手で問題を解いてみよう.
-
より完全に近いモデルの方が良いかというとそうでもない.というのも,より細かい挙動を説明できるモデルはそれだけ細かくて複雑になることが多く,本質をつかみにくくなる.また,複雑な分だけ解析も難しいし,パラメータも増えるし,シミュレーションも大変だ.モデルは「知りたい挙動が上手く説明できる範囲で」簡潔な方が良いと言っていいだろう. ↩︎
-
Euler 法とはスキームを作る手法のことなのだが,得られた数値スキームである Euler スキームのことを Euler法と呼ぶ人もいる. ↩︎
-
これは教員が動かしたとき,In [18] にこのプログラムを書き込んだからだ. ↩︎
-
一応,
domainerror
についてはマニュアルの Core.DomainError に解説があり,
The argument val to a function or constructor is outside the valid domain.
と書いてある.わかりやすい例も載っている.見ておくと良いね. ↩︎ -
これは ^1.2 という計算がどうやって行われるかを考えればわかる. $x^{1.2} = ( e^{\log(x)} )^{1.2} = e^{1.2\log(x)}$ として計算されるので,計算 $\log(x)$ ができないときはエラーになる. そしてそれは $x \leq 0$ のときで,まさに今回の話になる. ↩︎
-
今回のようなシンプルな問題ならば簡単なのだが,一般にはパラメータ(今回は Δt)とスキームの数値安定性との関係を数学的に解析してその結果に基づいて決めるということになるので結構難しいのだ. ↩︎