Photo by Roman Mager on Unsplash

モデリングとは

モデリングとは，対象として考えている現象に対して数学などの信頼できる理論で説明して「理論仮説 = モデル」を作り上げることである． つまり，「こうなっている」という仮説をたてることである． 通常は，その仮説を検証できるように数学を用いてモデリングする(微分方程式を用いることが多いよ)．

なお，モデルは完全である必要はなく，あくまで「この現象のこの挙動をある程度説明できる」というもので良い¹．

モデリングに関する、わかりやすい古典的な題材を資料 モデリング古典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 法とは時間発展型の微分方程式問題に対してスキームを作成する方法²で、対象とする微分方程式の解が $u(t)$ とすると、 時間刻み幅 $\Delta t$ を導入した上で、

1. 近似解として数列 $\{u_n\}_n$ を $u_n \cong u(n\Delta t)$ を想定して導入する．
   
   
2. 微分方程式中の時間微分項 $du/dt$ を陽的差分 $( u_{n+1} - u_n ) / \Delta t$ で近似し、 他の部分に出てくる $u(t)$ はそのまま $u_n$ で近似する.

としてスキームを作成するものだ． 一回やってみるとわかるが、こうして得られたスキームは、 既知量 $u_n$ を代入するだけでただちに未知量 $u_{n+1}$ が求まる ような便利な形になっている．

既知量 $u_n$ を代入するだけでただちに未知量 $u_{n+1}$ が求まる形になっているスキームを(求めたい近似解が陽的に表されているという意味で)陽的スキームとよぶ．

実際、今回の対象としている微分方程式に 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$ と順番に求めていけば好きなだけ先の近似値を計算できることがわかる．

まず各種定数などを定義しよう

数字を与えておかないといけない変数を設定しておこう．

1
2
3
4
5
6
7
8
9

# 問題に含まれる定数
γ = 1.0

# 初期値
u0 = 0.1

# 数値スキームのパラメータ
Δt = 0.1 # 特に根拠のない，とりあえずの値だ．

Euler スキームを関数として定義しよう

その前に Julia で関数を定義する方法について説明しておこう． Julia では関数を定義する方法がおおよそ ３つある．

一番手堅い方法は function 命令を使う方法で、例えば以下のようになる．ちなみに、return 命令を使うと出力値を指定できる．return を使わない場合は、最後の式の結果が出力値になる．

1
2
3
4
5
6

# Euler スキームを関数として定義する
# 入力:古い近似値, 出力:新しい近似値

function Euler(u)
    return u + Δt * γ * (sin(u))^1.2
end

２つ目の方法はシンプルな関数向けの方法で，数学の数式で書くような方法で次のようにして定義するものだ．

1

Euler2(u) = u + Δt * γ * (sin(u))^1.2

３つ目の方法は写像の表現に似た無名関数という機能を使うもので，無名関数に名前をつける形で次のようにして定義する．

1

Euler3 = u -> u + Δt * γ * (sin(u))^1.2

上手く動くならどの方法でも良い．だから好きな方法を選べば良いよ．

蛇足だが、どの定義も本当に同じか、気になる人へ

上のいずれも同じ関数を定義していることを、実際に値を入れて確認してみよう．

1
2
3

u = 1.0

Euler(u), Euler2(u), Euler3(u)

(1.0812918438981398, 1.0812918438981398, 1.0812918438981398)

それでも動作の違いが有るのではないか、と不安な人は @code_lowered マクロを用いて、次のようにそれぞれの内部コードが同一であることを確認しても良い．

1

@code_lowered Euler(u)

CodeInfo(
1 ─ %1 = Main.sin(u)
│   %2 = %1 ^ 1.2
│   %3 = Main.Δt * Main.γ * %2
│   %4 = u + %3
└──      return %4
)

1

@code_lowered Euler2(u)

CodeInfo(
1 ─ %1 = Main.sin(u)
│   %2 = %1 ^ 1.2
│   %3 = Main.Δt * Main.γ * %2
│   %4 = u + %3
└──      return %4
)

1

@code_lowered Euler3(u)

CodeInfo(
1 ─ %1 = Main.sin(u)
│   %2 = %1 ^ 1.2
│   %3 = Main.Δt * Main.γ * %2
│   %4 = u + %3
└──      return %4
)

みると、示されている内容がなんとなくわかるだろう． そして、内部的に全く同等の動作をするはずだ、ということがこれを見ると分かるだろう．

定義した Euler スキームを使って、初期値から順番に近似値を計算しよう

さて話を戻して，Euler スキームを関数として既に定義したので，あとは単に順番に計算していけばいいだけだ． ちなみに，いつものようにグラフを描きたいので結果を配列に記録しておこう．

1
2
3
4
5
6
7
8
9

# Euler スキームを用いて、初期値から近似値を計算していく．

u = u0      # 初期値を値として設定．
u_sq = [ u ] # 結果を配列に貯めていく．

for n in 1:100      # とりあえず 100ステップで試してみよう．
    u = Euler(u)    # 古い u -> 新しい u
    push!(u_sq, u)  # push! で配列に要素を追加する．
end

結果のグラフを見てみよう

上の計算が終わると、配列 u_sq に $u_0, u_1, \cdots, u_{100}$ の要素が入っているはずである． そこでこの数列をプロットしてみよう．

いつものように、Plots パッケージを下記のようにして呼び出しておいて…

1

using Plots

データを貯めてある配列を直接 plot しよう．

1

plot(u_sq)

前にも少し示したが，グラフの横軸をきちんと時間 $t$ そのものにしたい場合は、$x$軸用の変数も作った上で, $x$軸用の変数と $y$軸用の変数(今回は u_sq)を plot 関数に一緒に渡せば良い．具体的には以下のような感じになる．

1
2
3
4
5

Rx = 0.0:Δt:100*Δt # 範囲というやつだね

plot(Rx, u_sq)
xaxis!("time t")
yaxis!("u")

前にも書いたが, 命令の最後に ! がついているものは、「追加で」という意味を示す．

ちなみに plot 命令のオプションとして xaxis 等を渡しても良い( "!" は外してね)．次のような感じだ．

1

plot(Rx, u_sq, xaxis = "time t", yaxis = "u" )

さらに、legend (図に書き込まれる「凡例」)が邪魔なら、legend オプションを適当にいじれば良い． 例えば次のようにすれば消せる．

1

plot(Rx, u_sq, leg = false)

さて，この場合グラフからも想像できるように，$u$ が下から $\pi$ に近づいていく. 数式をみて納得しておこう．

そこで、グラフに $\pi$ の値を示す「線」を追加しておこう． 定数関数を使ってもよいが以前も出てきた hline という命令でやってみよう．次のような感じだ．

1
2
3
4

plot(Rx, u_sq, leg = false)

# hline! には定数関数の値の配列を渡す．今回は中身が一つだけだが…
hline!([pi]) 

Euler法は失敗しやすい．見てみよう．

Euler 法は数式の変形も簡単だし、プログラムも簡単だ．実行速度も速い． しかしその分，(まだ説明していないが)結果の精度は悪いし、計算過程は不安定だ． ここでいう不安定というのは、数値解が発散したりすることを言う．

このことを、実際に見てみよう．やり方は簡単で、時間発展を追いかける際の $\Delta t$ を少し大きめにするだけでいい…

1

Δt = 1.2 # 以前はこの値が 0.1 だった．少し欲張って大きくしてみる．

さて、さきほどと全く同じように計算してみよう．ただし、過程がわかりやすいように、途中結果をいちいち表示させよう．

1
2
3
4
5
6
7
8

u = u0      # 初期値を値として設定．
u_sq2 = [ u ] # 結果を配列に貯めておく．

for n in 1:100
    println(n,": ",u,": ",sin(u))  # 途中経過を画面に出す．
    u = Euler(u)     # 古い u -> 新しい u
    push!(u_sq2, u)  # push! で配列に要素を追加する．
end

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(x::Float64)
   @ Base.Math .\math.jl:37
 [2] ^
   @ .\math.jl:901 [inlined]
 [3] Euler(u::Float64)
   @ Main .\In[2]:5
 [4] top-level scope
   @ .\In[21]:6
 [5] eval
   @ .\boot.jl:360 [inlined]
 [6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base .\loading.jl:1116

エラーを示して停まってしまう．何が起きた？

エラーが起きて停まった． プログラムを書いて動かす場面ではエラーはよく出てくるものであるから，これと付き合えないと話が進まない．

今回はせっかくなので，エラーメッセージの読み方を知って，状況を把握できるようになっておこう．

最初に，上のエラーメッセージの stacktrace という部分を見よう． ここには，エラーが起きた箇所から始まって，その動作を呼出した場所へと逆上っての情報が記されている．

真面目にいくならば，[1] から順に読んでいって，まずは top-level scope という文字が出てくる箇所を探そう． top-level scope というのはプログラムを動かしている一番表面ということで, この場合は jupyter に書き込んだプログラムそのものなのだ．

実際，この場合，

[4] top-level scope
   @ .\In[21]:6

というものが見つかる． そしてこの @ 以降を読む． これが，このプログラムのエラー箇所の「一番表面的な箇所」になる.

この場合 In[21]:6 となっていて， In[21] は jupyter の入力の In [21] という箇所³だ，ということを意味し，そして最後の :6 はその 6行目だ，ということを表している．

実際，上の入力の 6行目をみると u = Euler(u) となっているので，自分で定義した関数 Euler でエラーが出た，ということがわかる．

次に，この場所から stacktrace を上へ向かって読み進めて状況を細かく把握していく． この場合，

[3] Euler(u::Float64)
   @ Main .\In[2]:5

となっていて，関数 Euler の 5行目，つまり return u + Δt * γ * (sin(u))^1.2 でエラーになったということがわかる．

次にさらに stacktrace を上へ向かい，

[2] ^
   @ .\math.jl:901 [inlined]

を見ることでエラーの起きた行の ^ という箇所，つまり，sin(u)^1.2 の 1.2乗の計算でエラーになったということが判明する．

さらに stacktrace を上に行くと

[1] throw_exp_domainerror(x::Float64)
   @ Base.Math .\math.jl:37

とあり，「exp_domainerror というエラーで⁴，処理をした」と言っている．

そしてこのエラーについてはこのエラーメッセージ中に解説が書いてあるのでそれを読めば良い(stacktrace の追いかけはここで終了だ)．

そしてそれは DomainError with -0.007573793605922763: で始まる 3行のエラー説明部分である． 読んでみると，

Exponentiation yielding a complex result requires a complex argument.
意訳: 結果が複素数になる指数計算は引数が複素数であるべき

と言っている(さらにちょっとした例も載っているね)．

さて，stacktrace の追いかけと，このメッセージの両方から考えると， sin(u)^1.2 の ^1.2 の結果が複素数になってしまった，だからエラーだ，ということだ． となると ^1.2 の引数の sin(u) が怪しい． 具体的には，sin(u)が負の値になり，エラーの原因になったのでは？ と怪しまれる⁵．

実際，プログラムが出力した直前の表示

7: 3.149166519606018: -0.007573793605922763

の右端の数字が sin(u)の値で確かにマイナスになっている．

しかもこの数字 -0.007573793605922763 はエラーメッセージにも現れている． よって、エラーは

1

(-0.007573793605922763)^1.2

を計算しようとして出たものと考えてよいはずだ．

確認のため，実際に上の１行計算を実行させてみよう．そうすると、

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:
 …略…

と出力され，先のものと全く同じエラーメッセージが現れる． この推測が正しかったのだな．

さて、エラーの原因が sin の値がマイナスになったことだとはっきりとわかったので、 次に，なぜマイナスになったのか，数字を丁寧に追いかけてみよう．

すると、7回目の出力の $u$ が $\pi$ をわずかに超えている様子がみてとれる． しかしそもそも、初期値が小さな正の実数ならば,もとの微分方程式では解は $\pi$ 未満にとどまるためこれはあり得ないシチュエーションである．

どうやら、Euler 法の誤差によって近似解が $\pi$ 以上になってしまうという状況が発生し，そのために sin(u) が負になってしまった様子だ．

そのあたりを、plot を使って確かめていこう．

まず，$\Delta t$ が小さいときの結果とどう違うのか．時間軸を丁寧に揃えて、同時に表示してみよう． まず、近似解の値そのものを比べてみよう．

1
2
3
4
5
6
7
8
9

default(leg = false) # leg = false を共通設定に．

t_sq2 = [ n * 1.2 for n in 0:6 ] # plot にわたすのは範囲でなくて配列でも良い．
plot(t_sq2, u_sq2)  # Δt = 1.2 の方

t_sq = [ n * 0.1 for n in 0:72 ]
plot!(t_sq, u_sq[1:73]) # Δt = 0.1 の方

hline!([pi])

見ると分かるが、「失敗したほう, つまり粗い方」は解が大きくなっていくスピードにブレーキがかかるのが大変に遅れる、という違いになっている．そのため、突破するはずのない値(この問題の場合は $\pi$ のこと)を近似解が突破してしまう、という格好だ．

その結果、$\sin(u)$ がありえない値である「負」になってしまうのだ．それもグラフで比較してみよう．

1
2
3

plot(t_sq2, sin.(u_sq2))
plot!(t_sq, sin.(u_sq[1:73]) )
hline!([0.0])

うまく計算できている方はこの値が「正のまま」であるが、駄目な方は「負」になってしまっている．これを 1.2 乗しようとするもんだから、プログラムが動かなくなってしまっている，ということがこの図でも確認でき，これでエラーの正体がきっちり把握できた，ということになるだろう．

じゃあ Δt はいくつなら動くのか？

このエラーの経験から，「ならば Δt はいくつにするべきなのか？」という疑問が生ずると思う． 実はこれは「数値安定性を確保するためのパラメータの決定」という問題⁶で，一般には結構難しい．

ただまあ今回はスキームが大変シンプルなので，実は結構簡単に調べがつく． やってみよう．

まず，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}} \]

となる．

よって，この右辺の最小値を調べ，$Δt$ をその値より小さくしておけば安全，ということになる．

そこでこの右辺を \[ Up(u) = \frac{\pi - u}{\gamma\,\,\sin(u)^{1.2}} \label{up:up-u} \] と定義して，まずはこのグラフを描いてみよう．

1
2
3

Up(u) =  (π - u )/ ( γ * sin(u)^(1.2) )

plot( 0:0.01:π, Up ) # u_n の範囲は 0 より大きく π より小さいので．

グラフは左端が一番大きく，よーく見ると右端も少し上がっている． そこでもう少しわかりやすいように範囲を狭めてグラフを描いてみよう．

1

plot( 1.0:0.01:π, Up )

このグラフに依ると $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 にて，
• 題名を 2021-numerical-analysis-report-05 として,
• 教官宛(アドレスは web の "TOP" を見よう)に，
• 自分の学籍番号と名前を必ず書き込んで，
• 内容はテキストでも良いし，pdf などの電子ファイルを添付しても良いので，

下記の内容を実行して，結果や解析，感想等をレポートとして提出しよう．

1. 上の $(\ref{eq:u-min})$ 式の $u_{\mbox{min}}$ の値を自分でも求めてみよう．なお，近似精度はあまり気にしなくて良い．
   
      やり方はいろいろある． 一番オーソドックスなのは式 $(\ref{up:up-u})$ を $u$ で微分した式を手で求めておいて，それがゼロになるところを求める方法だろう．これなら前回までのゼロ点求解の知識が使える．
   他には，たとえば x in 2.3:0.01:2.5 の範囲で数値的に最小値を探す，なんて強引な方法もあるだろう． いろいろ工夫してみるとよいだろう．
   
   
2. $\Delta t$ を大きくしたい場合は？
   
      今回の Euler スキームの場合は $\Delta t$ が大きいと近似解がうまく得られないことがわかった． では、$\Delta t$ を大きくとりたいときはどうしたらよいだろうか．
   次回以降にさらにいろいろな手法を紹介するけれども、まずは自力で少し考えてみよう(調べてもよいが、できるだけ自分で考えてみよう)． $\Delta t$ を小さくするという問題回避方法以外にどんな解決方法がありうるだろうか． そして、実際にそれを実行してうまくいくか確認してみよう．
   馬鹿げた方法と思うようなものでも良いから、アイディアを出し、実際にその方法にチャレンジしてみて、その結果をレポートしよう．
   
   
3. 似て非なる、解きにくい問題に取り組んでみる. \[ \frac{d u}{d t} = - \sqrt{ u } \] という常微分方程式に基づく時間発展問題を、たとえば 時間 $t=0$ での 初期値 $u_0 = 0.25$ のもとで考えよう．
   
   このとき、Euler スキームを構成し、$\Delta t$ をいろいろ変えて計算して結果を見比べることで以下の問にチャレンジしよう．

• 厳密解があるとして、どんな形であると推測されるか．推測してみよう．
• Euler 法で実質的に $t \in [0, T)$ で近似解が求めることができるとしよう．その場合，$T$ はいくつぐらいか．そしてそれはどうしてその値になるのか
• (チャレンジ問題) 厳密解を式で書けるか. 手で問題を解いてみよう．