[[FrontPage]]
#contents
2011/06/19からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/9/
また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。
* 微分方程式 [#tdd3f9bf]
Sageを使って微分方程式を解く方法を紹介します。微分方程式の性質から解析的な解が求まらない場合等に備え、
数値的に微分方程式を方法についても取り上げます。
** desolveを使った微分方程式の解法 [#od6e5d1e]
最初にSageで微分方程式を解くときに一般的に使用されるdesolve関数の使い方について説明します。
例として以下のような2次の微分方程式をSageで表してみます。
$$
\ddot z + \dot z + z = u(t)
$$
ここで外部からの変位u(t)を、\(u(t) = 1, t \ge 0\)のステップ関数とすると上記の微分方程式は
以下のように定義されます。
desolveでは、変数tをvar関数で、変数tの関数ｚをfunction関数で定義し、
微分方程式deに\(\ddot z + \dot z + z = 1\)の関係式をセットし、
初期条件icsに\([ t_0, z(t_0), \dot z(t_0)]\)を与えて解きます。
sageへの入力:
#pre{{
t = var('t')
z = function('z', t)
de = diff(z, t, 2) + diff(z, t) + z == 1
sol = desolve(de, z, ics=[0,0,0]); view(sol)
}}
&ref(out1.png);
求まった解\(z(t)\)とその微分\(\dot z(t)\)をプロットしてみます。
\(z(t)\)は、最初は振動していますが、t=10になると定常状態に収まって１となっています。
sageへの入力:
#pre{{
plt_sol1 = plot(sol, [t, 0, 20], linestyle ='-')
plt_sol2 = plot(diff(sol,t), [t, 0, 20], linestyle =':')
(plt_sol1 + plt_sol2).show()
}}
&ref(out2.png);
** 伝達関数H(s)を使った解法 [#j94e85fd]
次に、ラプラス変換と伝達関数を使って同じ微分方程式を解いてみましょう。
$$
\ddot z + \dot z + z = u(t)
$$
の両辺をラプラス変換します。zのラブラス変換をZ, uのラブラス変換をUとすると以下の様になります。
$$
s^2 Z + s Z + Z = U
$$
ここで伝達関数\(H(s) = \frac{Z}{U}\)を計算します。
$$
H(s) = \frac{1}{s^2 + s + 1}
$$
単位ステップ応答のラブラス変換は、\(\frac{1}{s}\)なので、\(Z(s)\)は以下の様に求めまります。
$$
Z(s) = H(s) U(s) = H(s) \frac{1}{s} = \frac{1}{s^3 + s^2 + s}
$$
最後に\(Z(s)\)をラプラス逆変換すると\(z(t)\)の解析解が求まります。
上記の処理をSageを使って行ってみます。
sageへの入力:
#pre{{
(s, t) = var('s t')
P = 1/(s*s + s + 1)
# 単位ステップ応答は伝達関数/sの逆ラプラス変換で得られる
ip = inverse_laplace(P/s, s, t)
view(ip)
}}
&ref(out3.png);
*** 結果のプロット [#je21b4ab]
求まった解をプロットするとdesolveの同じ結果と同じ結果を得ることができました。
sageへの入力:
#pre{{
plt_1 = plot(ip, [t, 0, 20], linestyle ='-')
plt_2 = plot(diff(ip,t), [t, 0, 20], linestyle =':')
(plt_1 + plt_2).show()
}}
&ref(out4.png);
** 微分方程式の数値解法 [#vd44088a]
*** desolve_system_rk4を使った解法 [#c17696e4]
desolve_system_rk4は、ルンゲクッタ法を使って効率よく微分方程式の数値解をもとめる関数です。
例題の微分方程式を\(\ddot z, \dot z\)について見てみます。
$$
\ddot z + \dot z + z = u(t)
$$
以下のような等価な連立微分方式を得ることができます。
$$
\left\{\begin{array}{l}
\dot z = \dot z\\
\ddot z = - z - \dot z + u(t)
\end{array} \right.
$$
これをdesolve_system_rk4関数を使って定義すると以下のようになります。
初期条件icsは\([ t_0, z(t_0), \dot z(t_0)]\)として解きます。
sageへの入力:
#pre{{
# 連立一次常微分方程式を解く関数desolve_system_rk4を使った例
(z,z_dot,t) = var('z z_dot t')
P=desolve_system_rk4([z_dot,-z - z_dot + 1],[z,z_dot],ics=[0,0,0],ivar=t,end_points=20)
plt_f = list_plot([(i,j) for i, j, k in P], plotjoined=True, linestyle ='-')
plt_g = list_plot([(i,k) for i, j, k in P], plotjoined=True, linestyle =':')
(plt_f+plt_g).show()
}}
&ref(out5.png);
*** 連立微分方程式の行列とベクトルを使った表現 [#t6e08df3]
先の連立微分方程式を行列とベクトルで表現すると以下のようになります。
&ref(eq1.png);
ここでベクトル\(x(t)\)を以下のように定義します。
$$
x(t) = \left[\begin{array}{c}
z(t)\\
\dot z (t)
\end{array}\right]
$$
\(x(t)\)で連立微分方程式を表すと以下のように簡単になります。
&ref(eq2.png);
*** オイラー法を使った解法 [#db09f249]
オイラー法は、以下のように行列とベクトルで表現された連立一次常微分方程式を差分近似式を使って解く方法です。
$$
\dot x(t) = A x(t) + B u(t), x(0) = x_0
$$
これを差分近似で表現すると以下のようになります。
$$
x(t + \Delta t) \simeq x(t) + (A x(t) + B u(t)) \Delta t
$$
それでは、オイラー法を使って先の微分方程式を解いてみましょう。
A, Bを以下のように定義し、刻み幅\(\Delta t = 0.05 \), 終了時刻\(t_f = 20 \)で計算します。
&ref(eq3.png);
sageへの入力:
#pre{{
# オイラー法での解法、インタラクティブにdtを変更できるようにした例
from pylab import arange
A = matrix([[0, 1], [-1, -1]])
B = vector([0, 1])
u = 1
tf = 20
@interact
def _(dt=selector([0.05, 0.1, 0.5, 1])):
# 微分方程式をEuler法で解く
x = vector([0, 0])
xx = []
for t in arange(0, tf, dt):
xx.append(x.list())
dx = A*x + B*u
x = x + dx*dt
# 結果をプロット
t = arange(0, tf, dt)
x1_lst = [v[0] for v in xx]
x2_lst = [v[1] for v in xx]
lst_plt1 = list_plot(zip(t, x1_lst), plotjoined=True, linestyle ='-')
lst_plt2 = list_plot(zip(t, x2_lst), plotjoined=True, linestyle =':')
(lst_plt1 + lst_plt2).show()
}}
&ref(out6.png);
** コメント [#q96a387f]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
