Photo by Roman Mager on Unsplash

モデリングとは

  モデリング = 仮説をたてること

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

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

モデリングに関する、わかりやすい古典的な題材を資料 モデリング古典3題 (pdf)で紹介しておこう．

単純な常微分方程式を Euler 法で計算してみよう

対象の微分方程式

対象として, 上の資料に載っている人口問題を考えよう． そして例えば，$u = u(t)$ ($t$: 時間) に対する常微分方程式 \begin{equation} \frac{du}{dt} = \gamma\, \left( \sin(u) \right)^{1.2}, \,\, \gamma > 0, t \in (0,\infty) \end{equation} をモデルとする仮説をたてた，としよう． ただし、$u$ の初期値 $u(0)$ については $0 < u(0) < 1$ としておく．

早速，このモデルのシミュレーション(つまり数値計算)を試みよう． まず，時間刻み幅 $\Delta t$ を一定にして、微分を差分で近似することで近似計算式を作ってみる．

そのために，$u_n \cong u(n\Delta t)$ という近似値 $u_n$ を用いて話を進めよう．

なお，もとの問題に対して実際に計算が可能な数値(近似)計算式のことを数値スキームもしくはスキームとよぶ．以降この単語が時々出てくるので覚えておこう．

さて、「もとの問題を近似してれば良いよね」と大雑把に考えるとスキームはの作り方は無数にある．

  Euler法 = 微分の近似以外は「知っている値」だけ使う方法

その中で、ここでは最も計算やプログラミングが簡単な 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 法を適用してみると、

\begin{equation} \frac{u_{n+1} - u_n}{\Delta t} = \gamma\, \left( \sin(u_n)\right)^{1.2} \end{equation}

となる. 左辺に新しい値による項を集め，右辺に古い値による項を集める形で整理するとさらにわかりやすく

\begin{equation} u_{n+1} = u_n + \Delta t \,\, \gamma \,\, \left( \sin(u_n)\right)^{1.2} \label{eq:euler-scheme} \end{equation}

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

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

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

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

1# 問題に含まれる定数
2γ = 1.0
3
4# 初期値
5u0 = 0.1
6
7# 数値スキームのパラメータ
8Δt = 0.1 # 特に根拠のない，とりあえずの値だ．

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

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

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

1# Euler スキームを関数として定義する
2# 入力:古い近似値, 出力:新しい近似値
3
4function Euler(u)
5    return u + Δt * γ * (sin(u))^1.2
6end

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

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

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

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

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

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

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

1u = 1.0
2
3Euler(u), Euler2(u), Euler3(u)

(1.0812918438981398, 1.0812918438981398, 1.0812918438981398)

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

1@code_lowered Euler(u)

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

1@code_lowered Euler2(u)

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

1@code_lowered Euler3(u)

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

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

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

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

1# Euler スキームを用いて、初期値から近似値を計算していく．
2
3u = u0      # 初期値を値として設定．
4u_sq = [ u ] # 結果を配列に貯めていく．
5
6for n in 1:100      # とりあえず 100ステップで試してみよう．
7    u = Euler(u)    # 古い u -> 新しい u
8    push!(u_sq, u)  # push! で配列に要素を追加する．
9end

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

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

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

1using Plots

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

1plot(u_sq)

  おめでとう！ 君はおそらく初めての自力でのシミュレーションに成功した！
仮説を数値で検証するスキルを手に入れたのだ．

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

1Rx = 0.0:Δt:100*Δt # 範囲というやつだね
2
3plot(Rx, u_sq)
4xaxis!("time t")
5yaxis!("u")

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

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

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

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

1plot(Rx, u_sq, leg = false)

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

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

1plot(Rx, u_sq, leg = false)
2
3# hline! には定数関数の値の配列を渡す．今回は中身が一つだけだが…
4hline!([π])  # π は \pi と入力してから Tabキーを押すと入力できる．

ちなみにここまでの情報を全てのせて，かつ，データは線ではなくて点で表すという場合は以下のようなコードと図になるね．

1plot(Rx, u_sq, marker = :circ,  xaxis = "time t", yaxis = "u", leg = false)
2
3# hline! には定数関数の値の配列を渡す．今回は中身が一つだけだが…
4hline!([π])  # π は \pi と入力してから Tabキーを押すと入力できる．

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

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

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

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

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

1u = u0      # 初期値を値として設定．
2u_sq2 = [ u ] # 結果を配列に貯めておく．
3
4for n in 1:100
5    println(n,": ",u,": ",sin(u))  # 途中経過を画面に出す．
6    u = Euler(u)     # 古い u -> 新しい u
7    push!(u_sq2, u)  # push! で配列に要素を追加する．
8end

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:41
 [2] ^
   @ .\math.jl:1206 [inlined]
 [3] Euler(u::Float64)
   @ Main .\In[2]:2
 [4] top-level scope
   @ .\In[17]:6

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

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

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

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

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

実際，この場合，

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

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

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

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

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

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

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

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

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

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

さらに stacktrace を上に行くと

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

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

そしてこのエラーについてはこのエラーメッセージ中に解説が書いてあるのでそれを読めば良い. なお，エラーの発生箇所にたどり着いたので stacktrace の追いかけはここで終了だ．

さてエラーの内容に話を戻すと， DomainError with -0.007573793605922763: で始まる 3行のエラー説明部分である． 読んでみると，

Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

意訳: 結果が複素数になる指数計算は引数が複素数であるべき.
 (複素計算であるならば) x^y を (x+0im)^y か Complex(x)^y といった表記に直せ．

と言っている.

さて，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$ が小さいときの結果とどう違うのか．時間軸を丁寧に揃えて、同時に表示してみよう． まず、近似解の値そのものを比べてみよう．

1default(leg = false) # leg = false を共通設定に．
2
3t_sq2 = [ n * 1.2 for n in 0:6 ] # plot にわたすのは範囲でなくて配列でも良い．
4plot(t_sq2, u_sq2)  # Δt = 1.2 の方
5
6t_sq = [ n * 0.1 for n in 0:72 ]
7plot!(t_sq, u_sq[1:73]) # Δt = 0.1 の方
8
9hline!([π])

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

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

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

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

じゃあ Δt をいくつにすれば安全にプログラムが動くのか？

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

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

まず，Eulerスキームの計算が直前までうまくいっている(つまり，$u_n < \pi$ ということ)状況下で今回のエラーを回避するための条件

\begin{equation} u_{n+1} < \pi \end{equation}

について今回の Euler スキーム $\eqref{eq:euler-scheme}$ を眺めて考える. するとそれは

\begin{equation} u_n + \Delta t \,\, \gamma \,\, \left( \sin(u_n)\right)^{1.2} < \pi \end{equation}

ということであるので，これを $\Delta t$ についての式に変形すると

\begin{equation} \Delta t < \frac{\pi - u_n}{\gamma\,\, \sin(u_n)^{1.2}} \end{equation}

となる．

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

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

1Up(u) =  (π - u )/ ( γ * sin(u)^(1.2) )
2
3plot( 0:0.01:π, Up ) # u_n の範囲は 0 より大きく π より小さいので．

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

1plot( 1.0:0.01:π, Up )

このグラフに依ると $u$ が 2.4 の付近？ に最小値があるように見受けられる． 実際，真面目に最小値を満たす $u$ を求めると⁷

\begin{equation} u_{\mbox{min}} = 2.446142705263389 \label{eq:u-min} \end{equation}

と求められ，そしてそのときの $Up(u)$ の値は

\begin{equation} Up(u_{\mbox{min}}) = 1.1864652653081016 \end{equation}

となる．つまり，

\begin{equation} \Delta t < 1.1864652653081016 \end{equation}

ならば今回のエラーを起こさずに計算が可能だ，ということになる．

実際，Δt = 1.186465 として計算しても確かに計算がうまくいくことなどが実際に確認できる．やってみよう．

Euler法の速度と精度を測定してみよう

次回，他の解法と比較してみるのだが，同じ問題でも解法によって速度や精度が異なる． 今回，まずはちょっと練習として Euler法の計算速度や精度を測定してみよう．

まずは速度だ

準備として、初期値と計算回数を入力すると最後の解を出力する形で関数を作ろう．以下のような感じだ．

1Δt = 0.1 # 念の為，値をもとに戻しておこう．
2
3function Euler_all(init_u, N)
4  u = init_u  # 最初の値
5  for n in 1:N
6    u = Euler(u)  # Euler スキームで新しい値を求める．
7  end
8  return u
9end

念のために、きちんと動くことを確認しよう．

1Euler_all(u0,100)

3.119420135870497

では、これを使って計算時間を比較しよう． それには、BenchmarkTools というパッケージに入っている、@benchmark というマクロが便利だ．

BenchmarkTools インストール (VIORAS 環境では不要)

自分の PC に自前で julia をインストールしているような場合はこの package は自分でインストールする必要がある(一回だけやれば良い)．

  VIORAS環境ではこの package のインストールは不要なので，下記のインストール作業を実行しないように！

個人の PC に Julia をインストールしているような場合は，以下のようにしてインストールしよう．

1using Pkg
2Pkg.add("BenchmarkTools")

BenchmarkTools を使って計算時間などを測定する！

さて，早速だがこのマクロを使って計算速度を測定、比較してみよう． まずパッケージの利用宣言を以下のようにして、

1using BenchmarkTools

あとは測定したい関数を @benchmark マクロに以下のように渡すだけだ．

1@benchmark Euler_all(u0, 100)

すると，次のように計算時間に関する簡単な統計情報を教えてくれる．

BenchmarkTools.Trial: 10000 samples with 3 evaluations.
 Range (min … max):   9.067 μs … 92.033 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     11.033 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   11.175 μs ±  1.564 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆▄          ▆▆▆█▆█▆▅▅▃▂▂▂▁▂▂▁▂ ▁▂ ▁▂                        ▃
  ██▅▃▃▆▆▅▁▃▁▆█████████████████████████▇█▆▆▅▃▅▁▄▅▅▄▃▄▅▆▅▆▇▆▇▆ █
  9.07 μs      Histogram: log(frequency) by time      16.2 μs <

 Memory estimate: 6.27 KiB, allocs estimate: 401.

結果を見ると、10000回のサンプルで, 平均が 11.2 μsぐらいの計算時間だということや，その他諸々がわかる． 真ん中のグラフは計算時間の頻度グラフだ． 通常は悪くても 16μs ほどで済むということなども読み取れる． なお，Jupyter の画面だとこのグラフに「青い線」「緑の線」が縦に入っていて，おそらくそれぞれ median値, mean値を意味する箇所だろうと思われる．

この速度が速いのか遅いのかは，次回に他の解法と比較して調べてみよう．今回はまずは速度測定の練習，という感じだ．

Benchmark 結果の中の GC は garbage collection (ゴミ処理) の頭文字で，使わなくなったメモリの整理等の作業を意味している． 状況によっては GC が計算時間の 99%以上を食ってしまうケースが有り，GC が発生するとかなり時間を食う． PC で使えるメモリが少ないほど GC が頻繁に発生するので，メモリが少ないとつらいね．

精度も見てみよう！

いきなりだが，time $t = 10$ でのこの問題の近似解値: $$u_{approx}(t=10) = 3.1184213196530965$$ を厳密解の代わりにしよう⁸．

そして，近似解とこの値との違いを誤差とみなして評価してみよう．

具体的には，異なる値の $\Delta t$ に対する誤差をグラフにプロットして様子を見るのだ． これでいろいろわかるはずだ．

まず，様々な $\Delta t$ に対して計算して $t = 10$ のときの誤差がどうなるのか，記録を取ろう． 少し長いが，下記のようなプログラムになるだろう．

 1# 全体条件
 2T = 10.0 # この時間の…
 3u_exact = 3.1184213196530965 # 厳密解もどき．
 4
 5ratio  = 1.5 # 計算回数を毎回 1.5倍する <=> Δt を毎回 2/3 に．
 6N_iter = 15  # この試行を何回やるか．
 7
 8# 初期化
 9N = 100
10Δt = T/N
11u0 = 0.1
12
13# 最初に一回計算しておく．
14e_euler = Euler_all(u0,N) - u_exact
15
16println("iter, N, Δt, e_euler")  # 画面表示
17println("0, $N, $Δt, $e_euler")
18
19# 記録をとりはじめる．
20sq_dt = [ Δt ] 
21sq_euler = [ e_euler ]
22
23# Δt を小さくしていったときの近似解をそれぞれ求める．
24for i in 1:N_iter
25  N = round(Int, ratio * N) # N を ratio 倍して，
26  Δt = T / N                # Δt は 1/ratio 倍になるはず．
27  push!( sq_dt, Δt )        # 記録をとる
28  
29  e_euler = Euler_all(u0,N) - u_exact     # 誤差を評価
30  push!(sq_euler, e_euler)        # 記録をとる 
31
32  println("$i, $N, $Δt, $e_euler")  # 画面表示
33end

すると次のような結果が得られる．

iter, N, Δt, e_euler
0, 100, 0.1, 0.000998816217400389
1, 150, 0.06666666666666667, 0.0006564740322243523
2, 225, 0.044444444444444446, 0.0004334550057212283
3, 338, 0.029585798816568046, 0.0002866709631739184
4, 507, 0.01972386587771203, 0.00019028442774038368
5, 760, 0.013157894736842105, 0.0001265708703082069
6, 1140, 0.008771929824561403, 8.421618377019158e-5
7, 1710, 0.005847953216374269, 5.607101881643928e-5
8, 2565, 0.003898635477582846, 3.734817716161132e-5
9, 3848, 0.002598752598752599, 2.4881098828366532e-5
10, 5772, 0.0017325017325017325, 1.65809781838e-5
11, 8658, 0.001155001155001155, 1.1051131361128341e-5
12, 12987, 0.00077000077000077, 7.366152340981813e-6
13, 19480, 0.000513347022587269, 4.910330449359179e-6
14, 29220, 0.00034223134839151266, 3.2733030281839604e-6
15, 43830, 0.00022815423226100844, 2.1820906250802352e-6

もう少し状況を把握するために，グラフにしてみよう．

 1default( marker = :auto, xaxis = ("Δt", :log10), 
 2  yaxis = ("error", :log10) )
 3# 共通の設定は default 宣言で一回だけすれば良い．
 4
 5plot(sq_dt, abs.(sq_euler), label = "Euler")
 6
 7# グラフを見やすくする為の補助線も描いてみる
 8default(marker = :none)
 9a1 = 2.0e-2
10p1(t) = a1 * t 
11plot!(sq_dt, p1.(sq_dt), label = "y = Cx")

この両対数グラフでのグラフの傾きは 1 であるので(目盛りで見てみよう)，Euler 法では $\Delta t$と誤差は比例している ことがわかる⁹．

レポート

下記要領でレポートを出してみよう．

• e-mail にて，
• 題名を 2024-numerical-analysis-report-05 として,
• 教官宛(アドレスは web の "TOP" を見よう)に，
• 自分の学籍番号と名前を必ず書き込んで，
• 内容はテキストでも良いし，pdf などの電子ファイルを添付しても良いので，

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

  注意
　　近年はセキュリティ上の懸念から，実行形式のプログラムなどをメールに添付するとそのメールそのものの受信を受信側サーバが拒絶したりする． そういうことを避けるため，レポートをメールで提出するときは添付ファイルにそういった懸念のあるファイルが無いようにしよう．

まあ要するに，レポートは 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. 似て非なる、解きにくい問題に取り組んでみる. \begin{equation} \frac{d u}{d t} = - \sqrt{ u } \end{equation} という常微分方程式に基づく時間発展問題を、たとえば 時間 $t=0$ での 初期値 $u_0 = 0.25$ のもとで考えよう．
   
   このとき、Euler スキームを構成し、$\Delta t$ をいろいろ変えて計算して結果を見比べることで以下の問にチャレンジしよう．

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