モデリング

モデリングとは

まずはモデリングそのものを解説しよう. モデリングに関する、わかりやすい古典的な題材を資料 モデリング古典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 スキームとよぶ.

実際、今回の対象としている微分方程式に 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} $$

となるので、実際はこの式を用いると良いだろう.

初期値は既知で $u_0 = u(0)$ とすることができるので、さらなる近似値を上の Euler スキームを用いて $u_1, u_2, \cdots$ と順番に求めていけば好きなだけ先の近似値を計算できることがわかる.

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

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

1
2
3
4
# 各種 定数の定義.Julia では a = b というのは a に b を代入するという意味.
u0 = 0.1  # 初期値
γ = 1.0
Δt = 0.1

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

Julia では関数を定義する方法がいくつかある.

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

1
2
3
4
# 古い値をもらって新しい近似値を求める Euler スキームを関数として定義する
function Euler(u)
    return u + Δt * γ * (sin(u))^1.2
end

Euler (generic function with 1 method)

単純な関数なら、数学の数式で書くような方法で次のようにしても定義できる.

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

Euler2 (generic function with 1 method)

写像の表現に似た無名関数という機能を使って次のようにしても定義できる.上手く動くなら好きな方法を選べば良い.

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

#3 (generic function with 1 method)

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

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

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( │2 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 ─ %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 ─ %1 = (Main.sin)(u) │ │ %2 = %1 ^ 1.2 │ │ %3 = Main.Δt * Main.γ * %2 │ │ %4 = u + %3 │ └── return %4 )

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

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

単に順番に計算していけばいいだけだが、あとでグラフを描きたいので、結果を配列に記録しておこう.

1
2
3
4
5
6
7
8
9
# Euler スキームを用いて、初期値から近似値を計算していく.グラフを描きやすいように、結果は配列に貯めていく.

u_sq = [u0] # 結果を貯めておく配列.最初は初期値しか入っていない.
u = u0      # その時点での最新の値.

for n in 1:100  # この for 文というもので繰り返しを行える.1:100 は数学でいう{1,2,...,100} と同じ.
    u = Euler(u)    # 古い u を右辺に入れて計算し、出てきた結果を左辺に代入している.
    push!(u_sq, u)  # push! というのは、配列に要素を追加する命令.
end

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

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

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

1
using Plots

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

1
plot(u_sq)  # plot で配列を単純にプロットしてみる

svg

tips: 上のグラフの横軸をベクトルの要素番号ではなく、きちんと時間 $t$ そのものにしたいという場合は、$x$軸用の変数も作った上で, $x$軸用の変数と $y$軸用の変数(今回は u_sq)を plot 関数に一緒に渡せば良い.具体的には以下のような感じだ.

1
2
3
4
5
t_sq = [n*Δt for n in 0:100] 

plot(t_sq, u_sq) # x軸の変数と y軸の変数を一緒に渡せば良い.
xaxis!("time t") # x軸について、xaxis! 命令を使うといろいろ変えられる.
yaxis!("u")      # y軸についても同様.

svg

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

ちなみに、上の xaxis! 等はオプションとして渡しても良い(! マークは外してね).次のような感じだ.

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

svg

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

1
plot(t_sq, u_sq, leg = false)

svg

$u$ が $\pi$ に近づいていくという話は既に授業で聞いたかもしれないね. そこで、グラフに $\pi$ の値を示す「線」を追加しておこう.定数関数を使ってもよいが、hline という命令があるので使ってみよう.次のような感じだ.

1
2
plot(t_sq, u_sq, leg = false)
hline!([pi]) # hline! には定数関数の値の配列を渡す.今回は中身が一つだけだが…

svg

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

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

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

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

1.2

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
# Euler スキームを用いて、初期値から近似値を計算していく.グラフを描きやすいように、結果は配列に貯めていく.

u_sq2 = [u0] # 結果を貯めておく配列.最初は初期値しか入っていない.
u = u0      # その時点での最新の値.

for n in 1:100  # この for 文というもので繰り返しを行える.1:100 は数学でいう{1,2,...,100} と同じ.
    println(n,": ",u,": ",sin(u))  # 何が起きているか分かるように、途中結果を画面に出すようにしておく.
    u = Euler(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(::Float64) at ./math.jl:35 [2] ^ at ./math.jl:782 [inlined] [3] Euler(::Float64) at ./In[2]:2 [4] top-level scope at ./In[26]:6

エラーを示して停まってしまう.何が起きたのか?

上のエラーメッセージを読むと、 DomainError with -0.007573793605922763 と言っていて、 x^y を複素数表現に替えるか他の方法を… とアドバイスが続いていることから、 べき計算 ^ の部分で、引数が想定外の数字 -0.007573793605922763 だった、というエラーだということがわかる.

この数字は直前の表示

7: 3.149166519606018: -0.007573793605922763

の右端の数字であることから、これは sin(u) の結果で、エラーは

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: [1] throw_exp_domainerror(::Float64) at ./math.jl:35 [2] ^(::Float64, ::Float64) at ./math.jl:782 [3] top-level scope at In[30]:1

と、全く同じエラーが出てくるので、推測が正しかっただろうことが確かめられる.

さて、エラーの原因が見えたので、何が起きたのか数字を丁寧に見よう.

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

どうやら、Euler 法の誤差によって近似解が $\pi$ 以上になってしまうという状況が発生してしまった様子だ.

そのあたりを、plot を使ってみてみよう.

1
plot(sin.(u_sq2))

svg

ふ~む、先の結果とどう違うのか.時間軸を丁寧に揃えて、同時に表示してみよう. まず、近似解の値そのものを比べてみよう.

1
2
3
4
5
6
7
t_sq2 = [ n * 1.2 for n in 0:6 ] 
plot(t_sq2, u_sq2, leg = false)  # こっちが駄目な方

t_sq = [ n * 0.1 for n in 0:72 ]
plot!(t_sq, u_sq[1:73], leg=false )  # こちらは上手くいった最初の結果.

hline!([pi])

svg

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

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

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

svg

うまく計算できている方はこの値が「正のまま」であるが、駄目な方は「負」になってしまっている.これを 1.2 乗しようとするもんだから、プログラムが動かなくなってしまっているのだ.

Report No.03:

  1. $\Delta t$ が大きいと?
       $\Delta t$ が大きいと上の Euler スキームでは近似解がうまく得られないことが上記の内容で理解できたと思う. さて、では、$\Delta t$ を大きくとりたいときはどうしたらよいだろうか.
       次回以降にさらにいろいろな手法を紹介するけれども、まずは自力で少し考えてみよう(調べてもよいが、できるだけ自分で考えてみよう). $\Delta t$ を小さくするという問題回避方法以外にどんな解決方法がありうるだろうか. そして、実際にそれを実行してうまくいくか確認してみよう.
    馬鹿げた方法と思うようなものでも良いから、アイディアを出し、実際にその方法にチャレンジしてみて、その結果をレポートしよう.

  2. 似て非なる、解きにくい問題.

\[ \frac{d u}{d t} = - \sqrt{ u } \]

という常微分方程式に基づく時間発展問題を、たとえば 時間 $t=0$ での 初期値 $u_0 = 0.25$ のもとで考えよう.

このとき、Euler スキームを構成し、$\Delta t$ をいろいろ変えて計算して結果を見比べることで以下の問にチャレンジしよう.

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