今回は微分方程式をプログラムで解いていきます。
答えを直接求める解析解、そして離散化させることによって求める数値解(オイラー法)、さらに精度の高いルンゲクッタ法でシミュレーションをしていきます。
解析解と数値解(オイラー法)を比べる
\begin{align}
\frac{dx}{dt}=-5x \tag{1}
\end{align}
上記の微分方程式を解析解(高校で習う紙の上で解くやつ)と数値解(離散化させてパソコンで解くやつ。)で比べたいと思います。数値解が解析解と近くなれば、シミュレーションが上手くいくと判断でき、誤差が大きいくなるとこの数値解は使えないと判断されます。
解析解
上記の微分方程式の解を、
\begin{align}
x(t)=Ce^{at}
\end{align}と仮定します。
すると、その微分は
\begin{align}
\frac{dx}{dt}=Cae^{at}=ax(t)
\end{align}
より$a$=-5だとわかります。
以上より一般解は
\begin{align}
x(t)=Ce^{-5t}
\end{align}
と分かります。
仮に初期値t=0の時に$x=10$とすると、
\begin{align}
10=Ce^{-5*0}=C
\end{align}
したがって解析解は
\begin{align}
x(t)=10e^{-5t}
\end{align}と求まります。
数値解
微分方程式
\begin{align}
\frac{dx}{dt}=-5x
\end{align}
をプログラムで解いていきます。まずは単純な近似を用いたオイラー法で解きます。
オイラーを用いると、次のように近似できます。
\begin{align}
\frac{dx}{dt} \approx\frac{x(t+\Delta T)-x(t)}{\Delta T}=-5x(t)\\
tをk\Delta Tとしてk=(0,1,2…)と離散化すると\\
\frac{x(k+1)-x(k)}{\Delta T}=-5x(t)\\
上記の式は次のように展開できます。\\
x(k+1)=x(k)+\Delta T(-5x(k))
\end{align}
これをプログラムで表すとしたのようになります。
import numpy as np
import matplotlib.pyplot as plt
Tsim = 3 # Simulation Time [s]
Ts = 0.1
a = -5
x0 = 10
# Analytical Solution
t_a = np.arange(0, Tsim + Ts, Ts)
x_a = x0 * np.exp(a * t_a)
# Numerical Solution (Euler Method)
Tn = round(Tsim / Ts)
t_s = np.zeros(Tn + 1)
x_s = np.zeros(Tn + 1)
x_s[0] = x0
t_s=np.linspace(0,Tsim,Tn+1)
for k in range(Tn):
x_s[k+1]=x_s[k]+Ts*(a*x_s[k])
# Plotting
plt.figure()
plt.plot(t_a, x_a, linewidth=2, label='Analytical')
plt.plot(t_s, x_s, linewidth=2, label='Numerical')
plt.grid()
plt.xlabel('Time [s]', fontsize=20)
plt.ylabel('x', fontsize=20)
plt.legend(fontsize=20)
plt.show()
今は初期値を10として設定しています。またサンプリング時間が0.1秒と荒くなっているため解析解と離れてしまっています。
そこでサンプリング時間を0.01秒と小さくしてみます。その結果がしたのように差が小さくなっています。
研究を行っていくとサンプリング時間が非常に重要なことがわかります。
リアルタイム制御の場合はサンプリング時間が0.05秒(50ms)でも遅く感じます。
できれば0.001~0.005秒(1ms~5ms)欲しいところです。
このサンプリング時間はセンサによって変わりますがセンサの仕様書にサンプリング時間または周波数として掲載されています。
サンプリング時間の逆数がサンプリング周波数というのは当たり前ですが、
身の回りものの周波数を知ると驚くことが多いです。
優秀なセンサですら数メガHzなのに、それに匹敵する周波数を出す機器がたくさん存在しており、
「お前、こんなに高い周波数だしていたのか」と周波数を聞いたときにわかる情報が増える感じです。
また並列化処理を施している場合はセンサからサンプリング時間を送っていたとしても
常に同じ周期かわからないため計算で使うときはPython上でのサンプリング時間(for分内でサンプリング時間を作る)を信用した方がいいケースもあります。
オイラー法とルンゲクッタ法の精度の違い
今度はオイラー法とより精度の高いルンゲクッタ法で(1)を解いていきます。
初期値は変わらずに10からスタートします。
下のグラフを見るとオイラー法はかくついていますがルンゲクッタ法はまだましであることがわかります。そして解析解とも一致しているように見えます
import numpy as np
import matplotlib.pyplot as plt
# Parameters of Simulation
Tsim = 10 # Calculation Time [s]
Ts = 0.2 # Sampling Time [s]
Tn = round(Tsim / Ts)
x1_01 = 10 # Initial Value of x
a = -5 # Parameter of Differential Equation
# Analytical Solution
t = np.arange(0, Tsim + Ts, Ts)
x_a = x1_01 * np.exp(a * t)
# Runge-Kutta method
def runge_kutta( x1_01, a):
Tn = round(Tsim / Ts)
x = np.zeros(Tn + 1)
x[0] = x1_01
for k in range(Tn):
p1 = a * x[k]
p2 = a * (x[k] + (Ts / 2) * p1)
p3 = a * (x[k] + (Ts / 2) * p2)
p4 = a * (x[k] + Ts * p3)
x[k + 1] = x[k] + Ts / 6 * (p1 + 2 * p2 + 2 * p3 + p4)
return x
# Runge-Kutta method simulation
runge_result = runge_kutta( x1_01, a)
# Euler Method simulation
x_s = np.zeros(Tn + 1)
x_s[0] = x1_01
for k in range(Tn):
x_s[k + 1] = x_s[k] + Ts * a * x_s[k]
# Plotting
plt.figure()
plt.plot(t, x_a, label='Analytical Solution')
plt.plot(t, x_s, label='Euler Method')
plt.plot(t, runge_result, label='Runge-Kutta Method')
plt.grid(True)
plt.xlabel('Time [s]', fontsize=12)
plt.ylabel('x', fontsize=12)
plt.legend(fontsize=12)
plt.show()
私はオイラー法やルンゲクッタ法を、何かを積分したいときによく使います。
例えばセンサ情報を積分したいときです。
しかしただ積分しただけですと積分誤差が蓄積してしまうため何か工夫を施す必要があります。
対処法を以下に３つ挙げます。
1.微小時間をかけ続けたものを積分したものから引く
2.ハイパスフィルタやローパスフィルタを移動平均フィルタとして使う(FIRフィルタ)
3. 求めたい値を別の手法でも計測し、その差にゲインをかけてフィードバックする。
(カルマンフィルタなど)
これらの手法についてはまた後々記事を書いていこうと思います。
関連記事
上のMEMO「1.微小時間をかけ続けたものを積分したものから引く」の詳細記事
コメント