Runge-Kutta 法を試してみよう¶

題材として、カオス現象を引き起こす有名な方程式の一つである 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)

初期値を適当に用意する

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

早速、計算してみよう

In [5]:
U = copy(U0)
sq_U = copy(U0)

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

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

In [6]:
sq_U
Out[6]:
3×3001 Array{Float64,2}:
 1.0  1.01257   1.04882   1.10721   1.18687   …  10.8724   10.3808    9.86686
 1.0  1.25992   1.524     1.79831   2.08854       6.10074   5.33638   4.64963
 1.0  0.984891  0.973114  0.965159  0.961737     34.8026   34.486    34.0768

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

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

ふ～む、少しわかりにくい挙動だな．では、相図にして見てみよう．

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

これは有名なバタフライ図というやつだな． 大きく分けて二つの状態を行き来していることが相図から見て取れる．

問題は、その状態間を行き来する条件が微妙そうなことだ． それが初期値の微妙な差の影響をどれだけ受けるか、ちょっと試してみよう．

初期値鋭敏性を見てみよう¶

そこで、初期値をいれたら解を出力する形で関数を作ろう． ちなみに、上の数値解をみると v(t) をみるとわかりやすそうなので、出力は v(t) だけでいいかな．

In [10]:
function v_out(init_U, N) # 初期値とステップ数だけを与えれば良いとしよう．
    U = copy(init_U)
 sq_U = copy(init_U)

for i in 1:N
    U = RK(U)
    sq_U = hcat(sq_U, U)
end
    
    return sq_U[2,:] # v の数値近似解のみを出力、でいいかな．
end
Out[10]:
v_out (generic function with 1 method)

ではまず、先の結果をあらためて．

In [11]:
result0 = v_out(U0,3000) # 初期値を U0, ステップ数 3000 での、v の数値解全体．
Out[11]:
3001-element Array{Float64,1}:
  1.0    
  1.25992
  1.524  
  1.79831
  2.08854
  2.40015
  2.73855
  3.10915
  3.51757
  3.96962
  4.47141
  5.02941
  5.65044
  ⋮      
 13.8324 
 13.1635 
 12.3935 
 11.5428 
 10.6353 
  9.69613
  8.7508 
  7.82309
  6.93398
  6.10074
  5.33638
  4.64963

初期値のある成分を 1% だけ変えてみて、影響をみてみよう¶

次に、初期値を少しだけ変えてみよう．

In [12]:
U1 = U0 + [0,0,0.01]    #  初期値の第3成分を 1% だけ変えてみる．
Out[12]:
3-element Array{Float64,1}:
 1.0 
 1.0 
 1.01
In [13]:
result1 = v_out(U1, 3000)  # そして、この初期値のもとでやはり 3000ステップほど計算．
Out[13]:
3001-element Array{Float64,1}:
  1.0    
  1.25982
  1.5238 
  1.798  
  2.08812
  2.3996 
  2.73785
  3.10829
  3.51652
  3.96834
  4.46988
  5.0276 
  5.64829
  ⋮      
  9.4971 
 10.3631 
 11.2916 
 12.2801 
 13.3234 
 14.4125 
 15.5338 
 16.6678 
 17.7881 
 18.8608 
 19.8438 
 20.6874

result0 と result1 は数字を見る限り、あとの方ではかなり異なるようだ． そうしたことを、プロットを通じて見てみよう．

In [14]:
plot( [result0 result1] , fmt = :png)  # matrix 形式で渡せば、同時にプロットしてくれる．
Out[14]:

どうやら、1600ステップあたりからズレがはっきりと出て、その後の挙動はお互いに似ても似つかない ことが見て取れる．

この二つの出力の関係性のなさ、を、次のプロットを通して感じ取ってみよう．

In [15]:
plot( result0, result1, fmt = :png ) #  横軸を result0, 縦軸を result 1 としている．
Out[15]:

上のプロットから横軸と縦軸の関係性が「読み取れれば」お互いの相関は強いし、「読み取れなければ」相関は薄い、ということになるのは直感的に分かるだろう．さて、上のプロットはどちらを示しているだろうか．

初期値の差を、0.1 % にしてみよう¶

In [16]:
result2 = v_out(U0 + [0,0,0.001], 3000) # やりかたは同じ．
Out[16]:
3001-element Array{Float64,1}:
 1.0    
 1.25991
 1.52398
 1.79828
 2.0885 
 2.4001 
 2.73848
 3.10907
 3.51746
 3.96949
 4.47126
 5.02923
 5.65022
 ⋮      
 2.79645
 3.05142
 3.21372
 3.30466
 3.34256
 3.34286
 3.31827
 3.27909
 3.23355
 3.18811
 3.1478 
 3.11646

先に関係性をみてみよう

In [17]:
plot( result0, result2, fmt = :png)
Out[17]:

どうやらこれでも、少し先へ行くと出力の関係性が失われるようだ． 念のため、時系列プロットもしてみよう．

In [18]:
plot( [result0, result2] , fmt = :png )
Out[18]:

どうやら、1750 ステップあたりで関係性が切れているようだ．

初期値の差を 0.01% にしてみよう¶

In [19]:
result3 = v_out( U0 + [0,0,0.0001], 3000)
Out[19]:
3001-element Array{Float64,1}:
  1.0     
  1.25992 
  1.524   
  1.79831 
  2.08854 
  2.40015 
  2.73854 
  3.10915 
  3.51756 
  3.9696  
  4.4714  
  5.0294  
  5.65042 
  ⋮       
 -0.1374  
 -0.139188
 -0.168239
 -0.218623
 -0.285468
 -0.364854
 -0.453708
 -0.549689
 -0.651086
 -0.756726
 -0.865882
 -0.978206
In [20]:
plot(result0, result3, fmt = :png)
Out[20]:

ふ～む、やはり関係性が失われるようだ．

In [21]:
plot( [result0, result3], fmt = :png )
Out[21]:

時系列データのプロットからもそれはみてとれる．この場合、2200ステップぐらいで関係性を喪失するようだ．

思い切って、初期値の差を 0.00001 % にしてみる…¶

In [22]:
result4 = v_out( U0 + [0,0,0.0000001], 3000)
Out[22]:
3001-element Array{Float64,1}:
   1.0    
   1.25992
   1.524  
   1.79831
   2.08854
   2.40015
   2.73855
   3.10915
   3.51757
   3.96961
   4.47141
   5.02941
   5.65044
   ⋮      
 -17.8923 
 -17.1442 
 -16.153  
 -14.9448 
 -13.5595 
 -12.0474 
 -10.4651 
  -8.87045
  -7.31716
  -5.85114
  -4.50757
  -3.30991
In [23]:
plot( result0, result4, fmt = :png )
Out[23]:

やはり関係性は失われる．図がさっぱりしているのは、関係性が失われる時刻が遅いからだろう．

In [24]:
plot( [result0, result4], fmt = :png )
Out[24]:

やはり、関係性が失われる時刻がだいぶ遅い．とはいえ、その後はもう全然関係ないと言って良い感じだ．

レポート課題¶

上の例で、初期値の差の小ささと、解の関係性が失われるまでの時間には、正の相関がありそうだ． つまり、初期値の差が小さくなればなるほど、解の関係性が失われるまでの時間は長くなる、ということだ．

そこで、この二つの量を大雑把にで良いので計測し、グラフを描いて、どういう関係があるのか、その関係式を推測してみよう．