対象として､ 授業時に示した、人口 $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 法 を使ってみよう． これは､ 微分方程式中の時間微分項 $du/dt$ を陽的差分 $(u_{n+1} - u_n)/\Delta t$ で近似し、他の部分に出てくる $u(t)$ はそのまま $u_n$ で近似するという方法である．
今回の対象としている微分方程式に 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$ と順番に求めていけば好きなだけ先の近似値を計算できることがわかる．
数字を与えておかないといけない変数を設定しておこう．
# 各種 定数の定義．Julia では a = b というのは a に b を代入するという意味．
u0 = 0.1 # 初期値
γ = 1.0
Δt = 0.1
Julia では関数を定義する方法がいくつかある．
以前示したように、一番手堅い方法は function 命令を使う方法で、例えば以下のようになる．ちなみに、return 命令を使うと出力値を指定できる．return を使わない場合は、最後の式の結果が出力値になる．
# 古い値をもらって新しい近似値を求める Euler スキームを関数として定義する
function Euler(u)
return u + Δt * γ * (sin(u))^1.2
end
単純な関数なら、数学の数式で書くような方法で次のようにしても定義できる．
Euler2(u) = u + Δt * γ * (sin(u))^1.2
写像の表現に似た無名関数という機能を使って次のようにしても定義できる．上手く動くなら好きな方法を選べば良い．
Euler3 = u -> u + Δt * γ * (sin(u))^1.2
単に順番に計算していけばいいだけだが、あとでグラフを描きたいので、結果を配列に記録しておこう．
# 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 がインストールされていないという場合は、それぞれ
Pkg.add("Plots")
Pkg.add("GR")
を一回だけやっておこう．
using Plots # グラフ描画フレームワークの利用を宣言する. juliabox だと既にインストールされているかも？
gr() # 具体的には GR というツールを背後で利用する．こちらもおそらくインストール済みかも？
plot(u_sq) # plot で配列を単純にプロットしてみる
上のグラフは横軸が $t$ そのものになっていないので困る、直したいという場合は、$x$軸用の変数も作った上で, $x$軸用の変数と $y$軸用の変数(今回は u_sq)を plot 関数に一緒に渡せば良い．具体的には以下のような感じだ．
t_sq = [n*Δt for n in 0:100] # julia だとこんな風にして配列を一度に生成できる．便利だね．
plot(t_sq, u_sq) # x軸の変数と y軸の変数を一緒に渡せば良い．
xaxis!("time t") # x軸について、xaxis! 命令を使うといろいろ変えられる．
yaxis!("u") # y軸についても同様．
さらに、legend (図に書き込まれる「凡例」)が邪魔なら、次のようにすれば消せる．
plot(t_sq, u_sq, leg = false)
$u$ が $\pi$ に近づいていくという話は既に授業で聞いたかもしれないね． そこで、グラフに $\pi$ の値を示す「線」を追加しておこう．次のような感じだ．
plot(t_sq, u_sq, leg = false)
hline!([pi]) # hline! には定数関数の値の配列を渡す．今回は中身が一つだけだが…
Euler 法は数式の変形も簡単だし、プログラムも簡単だ．実行速度も速い．しかしその分、結果の精度は悪いし、計算過程は不安定だ．ここでいう不安定というのは、数値解が発散したりすることを言う．
このことを、実際に見てみよう．やり方は簡単で、時間発展を追いかける際の $\Delta t$ を少し大きめにするだけでいい…
Δt = 1.2 # 以前はこの値が 0.1 だった．少し欲張って大きくしてみる．
さて、さきほどと全く同じように計算してみよう…
# Euler スキームを用いて、初期値から近似値を計算していく．グラフを描きやすいように、結果は配列に貯めていく．
u_sq = [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_sq, u) # push! というのは、配列に要素を追加する命令．
end
上のエラーメッセージを読むと、(負の数)^1.2 の部分で出ていて、そこで文句を言っている．
そしてさらに数字を丁寧に見よう．
すると、7回目の出力の $u$ が $\pi$ をわずかに超えている様子がみてとれる． しかしそもそも、初期値が小さな正の実数ならば,もとの微分方程式では解は $\pi$ 未満にとどまるためこれはあり得ないシチュエーションである．
どうやら、Euler 法の誤差によって近似解が $\pi$ 以上になってしまうという状況が発生してしまった様子だ．
次回以降に紹介するけれども、まずは自力で少し考えてみよう．