モデリング/微分方程式の出現 & Euler 法

現象・問題とモデリング

まずはモデリングそのものを解説しよう． モデリングに関する、わかりやすい古典的な題材を資料 授業資料: モデリング古典3題 (pdf) で紹介する．

Euler 法

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$ としておく．

Euler 法とはなにか、そしてその適用

時間刻み幅 $\Delta t$ を一定にして、微分を差分で近似することで近似計算式(スキームと呼ぶ)を作ることを考えよう． 以降、$u_n \cong u(n\Delta t)$ と想定した近似値 $u_n$ を用いて話を進めよう．

もちろん、近似式であるスキームの作り方は無数にあるが、ここでは上に述べた、最も計算やプログラミングが簡単な Euler 法 を使って近似式を作ってみよう．

この Euler 法であるが、定義はかなり簡単で、 常微分方程式中の時間微分項 $du/dt$ を陽的差分 $( u_{n+1} - u_n ) / \Delta t$ で近似し、他の部分の $u(t)$ はそのまま $u_n$ で近似することで、

既知量 $u_n$ を代入するだけでただちに未知量 $u_{n+1}$ が求まるように近似式を構成する方法

を言う…としておこう． そして、こうしてできた近似式を Eulerスキーム と言ったりする．

まあ、「今の近似値を代入するだけで次の時点での近似値が求まる形に、もとの微分方程式を近似する」と思っておけば十分だ．

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

ではプログラムへ

まず各種定数などを定義しよう． 具体的には、数字を与えておかないといけない変数を設定する、ということになる．

1
2
3
4

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

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

今の近似値を入力として入れると、次の時点での近似値が出力として出てくる「関数」として Euler スキームをプログラムしよう．

そして、Julia では関数を定義する方法がいくつかある．

一番手堅い方法は function 命令を使う方法で、例えば以下のようになる． 見てすぐわかるだろうから、文法的には難しいところはないはずだ．

あえて言うならば、return 命令を使うと出力値を指定できるし、return を使わない場合は、最後の式の結果が出力値になる、ということぐらいだろうか．

1
2
3
4

# 古い値をもらって新しい近似値を求める 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

この 3つの定義の仕方のうち、(上手く動くならば) 好きな方法で定義すれば良い．

なお、この「無名関数」の文法が必要になるシーンも有る．そうした場合はここを思い出そう．

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

上のいずれも同じ関数を定義していることを確認したいのであれば、そうだなあ、まず、実際に値を入れて確認してみよう．

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(:(begin nothing return u + Main.Δt * Main.γ * (Main.sin)(u) ^ 1.2 end))

1

@code_lowered Euler2(u)

CodeInfo(:(begin nothing return u + Main.Δt * Main.γ * (Main.sin)(u) ^ 1.2 end))

1

@code_lowered Euler3(u)

CodeInfo(:(begin nothing return u + Main.Δt * Main.γ * (Main.sin)(u) ^ 1.2 end))

この３つの結果を見ると一字一句変わらないので、内部的に全く同等の動作をするはずだ、ということが分かるだろう．

定義した 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 パッケージとそのサブパッケージ GR を

1
2

using Plots
gr()

と一回だけ呼び出しておいて plot しよう.

1

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

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軸についても同様．

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

1

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

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

1

plot(t_sq, u_sq, leg = false)

$u$ が $\pi$ に近づいていくという話は既に授業で聞いたかもしれないね． そこで、グラフに $\pi$ の値を示す「線」を追加しておこう．次のような感じだ．

1
2

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

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

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

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

1

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

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

 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

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

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

上のエラーメッセージを読むと、(負の数)^1.2 の部分で出ていて、そこで文句を言っている．

そしてさらに数字を丁寧に見よう．

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

どうやら、Euler 法の誤差によって近似解が $\pi$ 以上になってしまうという「ありえない」状況が発生してしまったために、数学上の問題が発生したのだ．

そのあたりを、plot を使ってもう少し詳しくみてみよう．

ちょっと面倒なので、先に、デフォルトで「マーカー有り、凡例無し」にしておこう．

1

default(marker = :auto, legend = false)

さて、うまくいくときといかないときどう違うのか, それぞれの数値近似解を同時に表示してみよう(もちろん時間軸を丁寧に揃えて)．

1
2
3
4
5

t_sq2 = [ n * 1.2 for n in 0:6 ] 
plot(t_sq2, u_sq2)  # こっちが駄目な方

t_sq = [ n * 0.1 for n in 0:72 ]
plot!(t_sq, u_sq[1:73] )  # plot! で重ね描きができる．こちらは上手くいった先のやつ．

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

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

1
2

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

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

本来は中間の時点での値をとるべきなのに古い時点の微分値を近似に一方的に使うことから、一般に言って、Euler 法は微分の効果が厳密な解よりも「しばらくの間強めに出る」か「しばらくの間弱めに出る」ということになりがちであり、その結果、近似解がおかしなことになることが多い．

Report No.03 似て非なる、解きにくい問題

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

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

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

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