FrontPage

2011/06/19からのアクセス回数 2281

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/9/

また、Sageのサーバを公開しているサイト(http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/)にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

微分方程式

Sageを使って微分方程式を解く方法を紹介します。微分方程式の性質から解析的な解が求まらない場合等に備え、 数値的に微分方程式を方法についても取り上げます。

desolveを使った微分方程式の解法

最初にSageで微分方程式を解くときに一般的に使用されるdesolve関数の使い方について説明します。

例として以下のような2次の微分方程式をSageで表してみます。 $$ \ddot z + \dot z + z = u(t) $$ ここで外部からの変位u(t)を、\(u(t) = 1, t \ge 0\)のステップ関数とすると上記の微分方程式は 以下のように定義されます。

desolveでは、変数tをvar関数で、変数tの関数zをfunction関数で定義し、 微分方程式deに\(\ddot z + \dot z + z = 1\)の関係式をセットし、 初期条件icsに\([ t_0, z(t_0), \dot z(t_0)]\)を与えて解きます。

sageへの入力:

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)

out1.png

求まった解\(z(t)\)とその微分\(\dot z(t)\)をプロットしてみます。

\(z(t)\)は、最初は振動していますが、t=10になると定常状態に収まって1となっています。

sageへの入力:

plt_sol1 = plot(sol, [t, 0, 20], linestyle ='-')
plt_sol2 = plot(diff(sol,t), [t, 0, 20], linestyle =':')
(plt_sol1 + plt_sol2).show()

out2.png

伝達関数H(s)を使った解法

次に、ラプラス変換と伝達関数を使って同じ微分方程式を解いてみましょう。 $$ \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への入力:

(s, t) = var('s t')
P = 1/(s*s + s + 1)
# 単位ステップ応答は伝達関数/sの逆ラプラス変換で得られる
ip = inverse_laplace(P/s, s, t)
view(ip)

out3.png

結果のプロット

求まった解をプロットするとdesolveの同じ結果と同じ結果を得ることができました。

sageへの入力:

plt_1 = plot(ip, [t, 0, 20], linestyle ='-')
plt_2 = plot(diff(ip,t), [t, 0, 20], linestyle =':')
(plt_1 + plt_2).show()

out4.png

微分方程式の数値解法

desolve_system_rk4を使った解法

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への入力:

# 連立一次常微分方程式を解く関数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()

out5.png

連立微分方程式の行列とベクトルを使った表現

先の連立微分方程式を行列とベクトルで表現すると以下のようになります。

eq1.png

ここでベクトル\(x(t)\)を以下のように定義します。

$$ x(t) = \left[\begin{array}{c} z(t)\\ \dot z (t) \end{array}\right] $$

\(x(t)\)で連立微分方程式を表すと以下のように簡単になります。

eq2.png

オイラー法を使った解法

オイラー法は、以下のように行列とベクトルで表現された連立一次常微分方程式を差分近似式を使って解く方法です。 $$ \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 \)で計算します。

eq3.png

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()

out6.png

コメント

選択肢 投票
おもしろかった 2  
そうでもない 0  
わかりずらい 0  

皆様のご意見、ご希望をお待ちしております。


(Input image string)


添付ファイル: fileeq3.png 392件 [詳細] fileeq2.png 376件 [詳細] fileeq1.png 388件 [詳細] fileout6.png 410件 [詳細] fileout5.png 392件 [詳細] fileout4.png 391件 [詳細] fileout3.png 424件 [詳細] fileout2.png 398件 [詳細] fileout1.png 460件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2013-01-21 (月) 16:41:18 (1735d)
SmartDoc