スポンサードリンク
モチベーション
先日、非線形力学系ローレンツ方程式のアトラクタ描画をルンゲクッタ法で計算しました。x,y,zに関する時系列データも簡単に計算できたので、描画してみたところ以下のようなグラフになりました。以下はx-tグラフです。
上記のグラフがカオス的な振る舞いを示しているというのは前回の記事に示した通りですが、これがルンゲクッタではなく、違う近似解法の場合どういった挙動を示すのか気になりました。
ネットを見る限りオイラー法とルンゲクッタ法の解の誤差について、示している記事もいくつかあったのですが、1変数の一階微分方程式や単純なバネのモデルのような1変数の二階微分方程式しか見つかりませんでした。
なので、今回は3変数関数のローレンツ方程式を用いて、数値解法が異なることによる、誤差の程度を検証したいと思います(定量評価はしてません)。
Eular法とは
微分方程式の近似解法の中でも比較的計算が簡単なものですが、その分、とても粗い近似解法です。
研究等ではまず使われないですが、数値解放の概念を理解する手助けにはなります。以下のような方程式を考えます。
\begin{eqnarray}
\frac{dy}{dt}=\dot{y}=f(t,y(t))
\end{eqnarray}
次に微分の定義を用いて、
\begin{eqnarray}
\dot{y}=\lim_{\Delta{t} \to 0} \frac{y(t+\Delta t)-y(t)}{\Delta t}
\end{eqnarray}
0<t<<1の時
\begin{eqnarray}
\dot{y}=\frac{y(t+\Delta t)-y(t)}{\Delta t}
\end{eqnarray}
と表せる。これを最初の式を用いて式変形していくと、
\begin{eqnarray}
\dot{y}\Delta t = y(t+\Delta t)-y(t) \\
y(t+\Delta t) = y(t)+\Delta tf(t,y(t))
\end{eqnarray}
となり、t秒の情報(y(t)、f(t,y(t))を用いて、Δt秒後のy(t+Δt)を求めるのがEular法です。
計算ではこのΔtを小さくするほど、微分に近づくことになるので精度はよくなりますが、その分、計算時間が長くなります。
Runge-Kutta法とは
ルンゲクッタ法は、計算に少々時間がかかりますが、その精度の高さから、微分方程式の数値解析には広く用いられています。
今回計算に用いているのは、4次のルンゲクッタ法です。4次のルンゲクッタ法には公式が存在するのですが、その公式の導出をしている人はどうやらいないみたいです。
公式の導出にはオイラー法とホイン法(2次のルンゲクッタ法)を用いなければいけないみたいですが、計算量がえげつないみたいで、そんなものを書こうとするマゾな人はいないらしいです。
簡単に言うとテイラー展開の4次までの項の係数と一致するように計算するらしいです。まあいろんな人が書いているので、ここで詳しく書く必要もないでしょう。
どうしても知りたい人はこのサイトを参考にしてください笑。→http://eman-physics.net/math/differential21.html
スポンサードリンク
計算結果比較
以下に、オイラー法とルンゲクッタ法で計算したx-tグラフを示します。
最初の方(t=0から500)くらいでは、二つの数値解法にはあまり差がありませんが、徐々に差は大きくなり、やがて全く異なる挙動になります。
試しにt=0から500の範囲のグラフを載せます
解に大きなズレはなくこの程度の範囲ならば、オイラー法でも良さそうです。
しかし実際の科学計算ではオイラー法はまず用いられません。研究で使用されている方程式はより高次で、高変数のものが多く、計算範囲ももっと広い場合がほとんどです。
なので、この程度の方程式かつこの程度の時間範囲しか解の正確性を保てないオイラー法はまず使われません。
オイラー法は比較的計算も楽で、簡単に理解できるので、数値解法の玄関的な役割として、教科書等に書かれているだけです。
ついでに、1変数関数であるならば、ある程度計算時間が長くてもオイラー法とルンゲクッタ法の解に差はあまり見られません。 以下は一次元の減衰系振動モデルのx-tグラフです。

若干紫っぽくなってるのは、赤と青が重なっているためです。一次元の簡単な運動方程式程度であれば、二つの解放に誤差が出ないことが観察できます。
今回用いたプログラム
-
ローレンツ方程式
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 |
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # lorenz equation def lorenz(x, y, z, p=10, r=28, b=8/3): x_dot = -p*x +p*y y_dot = -x*z +r*x-y z_dot = x*y -b*z return np.array([x_dot, y_dot, z_dot]) t=0 ##### 初期時間###### dt = 0.01 ##### 微分間隔###### stepCnt =500 ######時間ステップ## ########ルンゲクッタ#################### # Need one more for the initial values xs = np.empty((stepCnt + 1,)) ys = np.empty((stepCnt + 1,)) zs = np.empty((stepCnt + 1,)) # Setting initial values xs[0], ys[0], zs[0] = (0.1, 0.5, 1.2) k0=[0,0,0] k1=[0,0,0] k2=[0,0,0] k3=[0,0,0] listData=[0] for i in range(stepCnt): x,y,z=xs[i],ys[i],zs[i] k0 = dt * lorenz(x,y,z) k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.) k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.) k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2]) dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0 dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0 dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0 xs[i+1] = xs[i] + dx ys[i+1] = ys[i] + dy zs[i+1] = zs[i] + dz listData.append(i) ########オイラー法#################### # Need one more for the initial values xs2 = np.empty((stepCnt + 1,)) ys2 = np.empty((stepCnt + 1,)) zs2 = np.empty((stepCnt + 1,)) # Setting initial values xs2[0], ys2[0] ,zs2[0]= (0.1, 0.5,1.2) k02=[0,0,0] k12=[0,0,0] listData1=[0] for i in range(stepCnt): x,y,z=xs2[i],ys2[i],zs2[i] k02 = dt * lorenz(x,y,z) k12 = dt * lorenz(x+k02[0], y+k02[1] ,z+k02[2]) dx = (k02[0]+k12[0])/2.0 dy = (k02[1]+k12[1])/2.0 dz = (k02[2]+k12[2])/2.0 xs2[i+1] = xs2[i] + dx ys2[i+1] = ys2[i] + dy zs2[i+1] = zs2[i] + dz listData1.append(i) ###########all-time-datasiries(graph-setting)################## plt.plot(listData, xs, c='b', label="Runge-Kutta") plt.plot(listData1, xs2, c='r', label="Eular") ax = plt.gca() ax.set_ylabel("X(t)") ax.set_xlabel("time") ax.set_title("x-t graph") plt.legend() plt.show() |
-
一次元減衰振動モデル
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 |
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import math k = 1.0 m = 1.0 w = 0.01 tr=[] xr=[] te=[] xe=[] #############ルンゲクッタ################### def g(t, x): return np.array([x[1], (-k/m*x[0])-2*w*x[1]]) x = np.array([0.01, 0.0]) # 初期条件 x(0) = 1, v(0) = 0 dt = 0.1 for i in range(10000): t = i*dt k1 = g(t, x) k2 = g(t + dt/2, x + k1*dt/2) k3 = g(t + dt/2, x + k2*dt/2) k4 = g(t + dt, x + k3*dt) x += (k1 + 2*k2 + 2*k3 + k4)*dt/6 xr.append(x[0]) tr.append(t) ############オイラー法######################### x1 = np.array([0.01, 0.0]) # 初期条件 x(0) = 1, v(0) = 0 for i in range(10000): t = i*dt k1 = g(t, x1) k2 = g(t + dt, x1 + k1*dt) x1 += (k1 + k2)*dt/2 xe.append(x1[0]) te.append(t) fig = plt.figure() ax = plt.gca() plt.plot(tr, xr, c='b', label="Runge-Kutta") plt.plot(te, xe, c='r', label="Eular") ax = plt.gca() ax.set_ylabel("X(t)") ax.set_xlabel("time") ax.set_title("x-t graph") plt.legend() plt.show() |
ルンゲクッタ法とオイラー法で別々にリストを作って、別々にfor文を回しているので、すごい効率の悪いプログラムです。
if文等用いて変数を減らし、一つのfor文ですむように書き直したいと思います。
スポンサードリンク
最後に
今回は3変数関数におけるオイラー法とルンゲクッタ法の数値計算の誤差について述べました。
数値解法は他にもホイン法、ベルレ法、かえる跳び法、修正オイラー法等あるので、それらを用いて計算し、解の比較をしていきたいと思います。
ついでにルンゲクッタ法よりも僕は、計算回数と長期正確性を加味して、速度ベルレ法が一番おすすめです笑。