[[FrontPage]]
#contents
2012/04/03からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www.pwv.co.jp:8000/home/pub/23/
私のSageサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートをアップロードし、実行したり、変更していろいろ動きを試すことができます。
* Sageでグラフを再現してみよう：微分方程式（オイラー法）編 [#kd0f2935]
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。
今回は、島根大学電子制御システム工学科の吉田和信氏の
[[Octaveによる動的システムシミュレーション入門>http://www.ecs.shimane-u.ac.jp/~kyoshida/octave(2002).pdf]]
の例題2.1の図2.1を題材にします。
&ref(fig21.png);
** 例題2.1について [#y5692924]
以下のような連立一次常微分方程式をオイラー法を使って解く問題です。
$$
\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, uは以下の様な行列とベクトルで与えられています。刻み幅\( \Delta t \) = 0.05, 終了時刻\( t_f \) = 20とします。
&ref(eq1.png);
*** octave連携の準備 [#tec3c795]
sageからoctaveを使う場合のグラフ処理を助ける関数を定義します。
sageへの入力:
#pre{{
import time
# octaveのグラフ出力ターミナルをpngにセット
#junk = octave.eval('putenv("GNUTERM", "png")')
# octave3.6.1の場合コメントアウトしてください。
# sageとoctaveがインストールされているマシンで、octaveのグラフを確認しながら
# 実行する場合には、上記の設定をコメントアウトしてください。
# octaveのグラフをファイルに保存し、表示する関数を定義（画像のキャッシュをしない版）
def saveAndShowPlot(filename, fac=0.75):
output = DATA + filename
octave.eval("print -dpng '%s'"%output)
width = int(640*fac)
html('<img src="%s?%s" width="%spx">'%(filename, time.time(), width))
}}
** octaveを使ったオイラー法の解析例 [#da86c0f2]
最初は、octaveを使った解析例（解ex2_1.m）をそのままsageで実行します。ex2_1.mはワークシート上部のDataタグで内容を確認することができます。
ベクトル\(x\)は、\( z, \dot z \)を要素に持ちます。
$$
x(t) = \left[\begin{array}{c}
z(t)\\
\dot z (t)
\end{array}\right]
$$
図の青の線が、\(z(t)\)、緑の線が\(\dot z(t)\)に相当し、\(z(t)\)のピークで\(\dot z(t)\)が0になっているのが分かります。
sageへの入力:
#pre{{
file = DATA + "ex2_1.m"
junk = octave.eval('source("%s")'%file);
saveAndShowPlot('octave_euler.png');
}}
&ref(fig1.png);
** Sageでのオイラー法の解析例 [#id7aaad8]
次にSageを使って同じEuler法で解いてみます。ここではSageの特徴である、インタラクティブ機能を使って、
刻み幅dtを0.05, 01, 0.5, 1から選択できるようにしました。刻み幅によってグラフがどのように変化するかをみてください。
オイラー法は、連立一次常微分方程式を数値微分を使って解く方法であり、刻み幅によって解が正しく求まらないことがあることに注意が必要です。
sageへの入力:
#pre{{
# ex2_1.mをsageで記述し、インタラクティブに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(fig2-1.png);
&ref(fig2-2.png);
** Sageの関数を使った例 [#wc6f2842]
Sageにも連立一次常微分方程式を解く関数がいくつか提供されています。ここでは、
desolve_system_rk4を使った例を示します。
例題を\( z, \dot z \)で表すと、
&ref(eq2.png);
これを展開すると、以下の様になります。
$$
\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{{
# Sageの連立一次常微分方程式を解く関数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(fig3.png);
** 伝達関数H(s)を使った解法 [#sec1544e]
\(\ddot z, \dot z, z, u\)をまとめると以下の２次の微分方程式になります。
$$
\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(fig4.png);
*** 結果のプロット [#r34495e0]
求まった解をプロットすると図2.1と同じ結果となり、解析解と数値計算の結果が一致することが確認できました。
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(fig5.png);
** 微分方程式の解 [#dc295304]
最後に数式処理システムSageの微分方程式を解く関数desolveを使った例を示します。
desolveでは、変数tをvar関数で、変数tの関数ｚをfunction関数で定義し、
微分方程式deに\(\ddot z + \dot z + a = 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(fig6.png);
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(fig7.png);
** まとめ [#p1f81fe3]
行列とベクトルで表される連立一次常微分方程式をオイラー法で数値微分することによって、
n次（今回は2次）の常微分方程式の解を求めることができる。
さらにSageを使った解析解の求め方を紹介し、微分方程式をSageを使ってどのように解くかを図2.1
を例に説明しました。
是非、この機会にSageを使っていろいろな微分方程式を解いてみてください。
** コメント [#na90a136]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[4])
#vote(おもしろかった[2],そうでもない[0],わかりずらい[4])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
