線形多段解法 + 予測子修正子法を試してみよう¶

題材として、カオス現象を引き起こす有名な方程式の一つである Lorenz 方程式を扱ってみよう．¶

式としては以下の様な常微分方程式である． $$ \left\{ \begin{array}{rcl} \displaystyle \frac{d u}{dt} & = & \sigma (-u + v), \\ \displaystyle \frac{d v}{dt} & = & u (-w + r) - v, \\ \displaystyle \frac{d w}{dt} & = & u v - b w . \end{array} \right. $$

ただし、パラメータ(ザルツマンのパラメータ)は、 $ \sigma = 10, $ $ r = 28, $ $ b = 8/3 . $

In [1]:
const σ, r, b = 10, 28, 8/3  # まずパラメータを設定してしまおう. 定数の場合、const を宣言しておくと参照が速い！
Out[1]:
(10, 28, 2.6666666666666665)

方程式の右辺を定義してしまおう．

In [2]:
function f(U)
    (u, v, w) = U  # こうやって成分毎に受け取ることも出来る．
    
    return [
        σ * (-u + v)
        u * (- w + r) - v
        u*v - b*w
        ]
    
    end
Out[2]:
f (generic function with 1 method)

初期値の幾つかは Runge-Kutta 法によるので、その時間発展も、定義してしまおう．

In [3]:
Δt = 0.01

function RK(U)
    f1 = f(U)
    f2 = f(U + Δt/2 * f1)
    f3 = f(U + Δt/2 * f2)
    f4 = f(U + Δt * f3)
    return U + Δt * (f1 + 2*f2 + 2*f3 + f4)/6
end
Out[3]:
RK (generic function with 1 method)

Runge-Kutta 法を使って、$U_2,U_1,U_0, f_2, f_1, f_0$ を返す関数を作ろう

In [26]:
function make_init(U)
    V = copy(U)
    sq_V = copy(U)
    sq_f = f(U)
    
    for n in 1:2
        V = RK(V)
        sq_V = hcat(sq_V, V)
        sq_f = hcat(sq_f, f(V))
    end
    
    return sq_V, sq_f
end
Out[26]:
make_init (generic function with 1 method)

初期値を適当に用意する

In [4]:
U0 = [1.0, 1, 1]
Out[4]:
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

試しに、初期値がどうなっているか見てみよう

In [29]:
make_init(U0)[1]  # これが $U_0, U_1, U_2$ のはず
Out[29]:
3×3 Array{Float64,2}:
 1.0  1.01257   1.04882 
 1.0  1.25992   1.524   
 1.0  0.984891  0.973114
In [31]:
make_init(U0)[2]  # これが $f_0, f_1, f_2$ のはず．
Out[31]:
3×3 Array{Float64,2}:
  0.0       2.47351   4.75173 
 26.0      26.0947   26.8224  
 -1.66667  -1.35062  -0.996567

さて、線形多段解法 + 予測子修正子法のプログラムだ．

線形計算なので、線形代数の記号を使うと簡単にかける ことに着目しよう．

In [54]:
function LM(U, vF) # vF = Fn, Fn-1, Fn-2 のつもり．
    
    # P
    up = U + Δt* ( vF * [23/12, -4/3, 5/12] )
    # E
    fp = f(up)
    # C
    up = U + Δt* ( hcat(fp, vF) * [3/8, 19/24, -5/24, 1/24] )
    # E
    fp = f(up)
    
    return up, hcat(fp, vF[:,1:2]) # 入力と同じ型のデータで更新したものを返す
    end
Out[54]:
LM (generic function with 1 method)

早速、計算してみよう

In [55]:
sq_U, sq_f = make_init(U0) # 初期値を作成．

U = sq_U[:,3]
vF = copy(sq_f)

for i in 1:3000
    U, vF = LM(U,vF)
    sq_U = hcat(sq_U, U)
end

試しに結果の数字をちょろっと見てみる．

In [56]:
sq_U
Out[56]:
3×3003 Array{Float64,2}:
 1.0  1.01257   1.04882   1.07379   1.16296   …   2.87227   2.31572   1.8142 
 1.0  1.25992   1.524     1.78498   2.07339      -2.97887  -2.9685   -2.93744
 1.0  0.984891  0.973114  0.959705  0.957358     29.1586   28.3153   27.5102

あまり変なことにはなっていなさそうなので、プロットしてみよう．

In [7]:
using Plots
gr()
Out[7]:
Plots.GRBackend()
In [57]:
plot(sq_U', fmt = :png) # 描画点が多く、web browser に負荷が大きいので、画像形式を png にして負荷を下げよう
Out[57]:

大丈夫そうだ．では、相図にして見てみよう．

In [58]:
X, Y, Z = sq_U[1,:], sq_U[2,:], sq_U[3,:]
plot(X, Y, Z, fmt = :png)
Out[58]:

おお、大丈夫そうだな．

RK とどちらが速いかみてみよう¶

そこで、初期値をいれたら解を出力する形で両方の関数を作ろう．

In [61]:
# まずは Runge-Kutta 法．

function RKout(init_U, N) # 初期値と最終ステップ数だけを与えれば良いとしよう．
    U = copy(init_U)
 sq_U = copy(init_U)

for i in 1:N-1 # 正確には、これで最終ステップが N になる
    U = RK(U)
    sq_U = hcat(sq_U, U)
end
    
    return sq_U
end
Out[61]:
RKout (generic function with 1 method)
In [62]:
# 次に線形多段階法

function LMout(init_U, N) # 初期値と最終ステップ数だけを与えれば良いとしよう．
    sq_U, sq_f = make_init(init_U) # 初期値を作成．Runge-Kutta が2回動作する．
    U = sq_U[:,3]
    vF = copy(sq_f)

    for i in 1:N-3  # 正確には、これで最終ステップが N になる
      U, vF = LM(U,vF)
      sq_U = hcat(sq_U, U)
    end

    return sq_U
end
Out[62]:
LMout (generic function with 1 method)

試しに、両方動かしてみよう．データのサイズが一緒になることを確認しよう．

In [64]:
RKout(U0,3000)
Out[64]:
3×3000 Array{Float64,2}:
 1.0  1.01257   1.04882   1.10721   1.18687   …  11.3317   10.8724   10.3808 
 1.0  1.25992   1.524     1.79831   2.08854       6.93398   6.10074   5.33638
 1.0  0.984891  0.973114  0.965159  0.961737     35.0101   34.8026   34.486
In [65]:
LMout(U0,3000)
Out[65]:
3×3000 Array{Float64,2}:
 1.0  1.01257   1.04882   1.07379   1.16296   …   4.89749   4.16174   3.48692
 1.0  1.25992   1.524     1.78498   2.07339      -2.76237  -2.89064  -2.95742
 1.0  0.984891  0.973114  0.959705  0.957358     31.934    30.9673   30.0422

では、両方の速度を測ってみよう．それには、BenchmarkTools というパッケージに入っている、@benchmark というマクロが便利だ．

In [68]:
using BenchmarkTools 
# 使う準備をする．インストールしてない人は、Pkg.add("BenchmarkTools") としてインストールだ．
In [66]:
@benchmark RKout(U0,3000)  # まずは Runge-Kutta 法の速度を見よう．
Out[66]:
BenchmarkTools.Trial: 
  memory estimate:  109.94 MiB
  allocs estimate:  112265
  --------------
  minimum time:     35.006 ms (24.46% GC)
  median time:      46.338 ms (26.33% GC)
  mean time:        45.794 ms (25.77% GC)
  maximum time:     62.930 ms (26.28% GC)
  --------------
  samples:          110
  evals/sample:     1
In [67]:
@benchmark LMout(U0,3000)  # 次に、線形多段解法の速度をみよう．
Out[67]:
BenchmarkTools.Trial: 
  memory estimate:  108.48 MiB
  allocs estimate:  79319
  --------------
  minimum time:     36.275 ms (25.46% GC)
  median time:      44.255 ms (26.69% GC)
  mean time:        44.545 ms (26.47% GC)
  maximum time:     59.616 ms (27.84% GC)
  --------------
  samples:          113
  evals/sample:     1

ふむ、理論的には、微分値 $f(t,u)$ の評価回数が半分の線形多段解法が有理なはずなのだが、こんな単純な $f$ だと差がでないみたいだな…

レポート課題¶

線形多段解法で、Lorenz 方程式の初期値鋭敏性について一通り実験してみよ．