2015年度 数式処理実習 [第 12回] 連立微分方程式とベクトル場 1 連立

July 9, 2015
[p.1/4]
2015 年度 数式処理実習 [第 12 回] 連立微分方程式とベクトル場
1
連立微分方程式
前回は,1つの未知関数に対する微分方程式を解きました.今回は,2つの未知関数に対する連立微分方程式 (“ 微分
方程式系 ”とも言います) を解きましょう.
2つの未知関数 x(t), y(t) に対する連立微分方程式


 dx(t) = −y(t)
dt
(1)

 dy(t) = −x(t)
dt
について考えます.未知関数が2つで,微分方程式も2つです.
Maxima で (1) の一般解を求めましょう.いつも通り,su-shiki ディレクトリに移動してから,maxima を立ち上げ
ます.まず,変数 eq1 に第一の微分方程式を代入し,変数 eq2 に第二の微分方程式を代入します.
maxima (%i1) eq1:’diff(x(t),t)=-y(t);
変数 eq1 に (1) の第一の微分方程式を代入
(%i2) eq2:’diff(y(t),t)=-x(t);
変数 eq2 に (1) の第二の微分方程式を代入
単独微分方程式 (通常の 1 つの微分方程式) を解くには ode2 コマンドを使いましたが,連立微分方程式を解くには,
desolve というコマンドを用います.書式は desolve(連立微分方程式, 未知関数達); です.ここで,連立微分方程式は
リストで,つまり四角い括弧 [ と ] で囲んで与えます.未知関数達も同じくリストで与えます.
desolve で,連立微分方程式 (1) の一般解を求めると次のようになります.
maxima (%i3) desolve([eq1,eq2],[x(t),y(t)]);
- t
t
(y(0) + x(0)) %e
(y(0) - x(0)) %e
(%o3) [x(t) = ------------------- - -----------------,
2
2
従って,(1) の一般解は
連立微分方程式 (1) を解く
t
- t
(y(0) - x(0)) %e
(y(0) + x(0)) %e
y(t) = ----------------- + -------------------]
2
2
(y (0) + x (0)) e−t
(y (0) − x (0)) et
−
,
2
2
(y (0) − x (0)) et
(y (0) + x (0)) e−t
y (t) =
+
2
2
x (t) =
となります.ここでは,2 つの任意定数が %k1, %k2 などではなく,x(0) と y(0) で表されています.(ode2 の時とは違う
ことに注意しましょう.) この x(0) と y(0) に具体的な値を入れると特殊解が得られる,というわけです.
初期条件を与えて,それをみたすような特殊解を求めましょう.前回の単独微分方程式の時には,まず ode2 を用いて
一般解を求めておき,一般解と初期条件を ic1 や ic2 に引き渡すことによって特殊解を求めました.
desolve コマンドを用いる場合は,少し考え方が違い,次のようにします.まず atvalue というコマンドで初期条件
を設定しておいてから,desolve コマンドを実行させます.これで初期条件をみたす特殊解を得ることができます.
具体的にやってみましょう.初期条件
x(0) = 0, y(0) = 4
をみたす連立微分方程式 (1) の特殊解を求めてみます.
「x(t) という式は t = 0 において 0 という値をとる」,
「y(t) という
式は t = 0 において 4 という値をとる」というわけですから,次のように設定します.
maxima (%i4) atvalue(x(t),t=0,0);
(%i5) atvalue(y(t),t=0,4);
このように設定した上で,再び連立微分方程式 (1) の解を求めてみると,次のようになります.
2015 年度 数式処理実習 [第 12 回] 連立微分方程式とベクトル場
July 9, 2015
[p.2/4]
maxima
(%i6) desolve([eq1,eq2],[x(t),y(t)]);
- t
t
t
- t
(%o6)
[x(t) = 2 %e
- 2 %e , y(t) = 2 %e + 2 %e
]
従って,初期条件 x(0) = 0, y(0) = 4 をみたす (1) の特殊解は
x(t) = 2e−t − 2et
y(t) = 2et + 2e−t
と求まりました.
なお,一度 atvalue コマンドで設定した初期値は,ずっと記憶されます.新しい問題を解く前に,
(%i6) kill(all);
と打ち込み,これらの設定を消去しておきましょう.
2
maxima
ベクトル場と解曲線
連立微分方程式 (1) のことを理解するには,解を図に表してみることが大切です.
t を時刻のパラメータとし,関数 x(t), y(t) を平面内の点 P の x 座標,y 座標と考えましょう.時刻
t が動くにつれて,


dx(t)
(
)
x(t)


点P
は平面内の曲線を描きます.時刻 t における,この点 P の速度ベクトルは  dt  となります.速
dy(t)
y(t)
dt
度ベクトルはもちろん曲線に接しますから,速度ベクトルのことを接ベクトルとも言います.
(
)
(
)
x(t)
x(t)
さて,t をパラメータとする曲線
が連立微分方程式 (1) をみたすと仮定しましょう.この時,曲線
y(t)
y(t)
は,連立微分方程式 (1) の解曲線である,と言います.
(
)
(
) (
)
a
x(t0 )
a
時刻 t0 において,この解曲線が点
を通るとします.つまり,
=
です.すると,連立微分方
b
y(t0 )
b
程式 (1) が成り立つことから,この点における曲線の接ベクトルは


dx(t) (
) (
)
 dt t=t 
−y(t0 )
−b

0 =
=
 dy(t) 
−x(t0 )
−a
dt t=t0
(
)
(
)
a
−b
となります.一言でいうと「解曲線が点
を通るなら,その点における接ベクトルは
となる」というこ
b
−a
とです.
(
)
(
)
x
−y
そこで,xy 平面上の各点
に対して,この点を始点とするベクトル
を,あらかじめ描いておきます.
y
−x
すると下の図のようになります.(ただし図を見易くするため,ベクトルの長さを縮小して描いています.) そして,こ
れらのベクトルに常に接しているような曲線を描けば,それが解曲線になる,ということが分かります.
y
8
4
0
-4
-8
-8
-4
0
4
8
2015 年度 数式処理実習 [第 12 回] 連立微分方程式とベクトル場
July 9, 2015
[p.3/4]
なぞっていけば手書きでもできそうですが,これを Maxima に描かせてみましょう.まずは,上の図を描きます.こ
のように,平面の各点に対して,その点を始点とするベクトルが一つ定まっている時,ベクトル場が与えられていると
言います.
y
8
ベクトル場は
plotdf
というコマンドで描きます.今,描きたいベクトル場の
)
(
−y
成分は
ですから,次のように打ち込みます.
−x
maxima (%i8) plotdf([-y,-x]);
右のような図が,新しいウィンドウで立ち上がります.
4
0
-4
-8
-8
今作ったベクトル場の図の上でマウスをクリックすると,そ
の点を通る解曲線 (このベクトル場に常に接するような曲線)
を赤い線で描いてくれます.t が正の方向に動くとき,どち
らの方向に動くかを示す矢印も書いてくれます.
図の上でカーソルを動かすと,右下に示されている座標の
(
)
5.5
数値も動きます.例えば,点
の辺り (大体でいい
6.5
です) でクリックしてみましょう.すると,右図のような赤
い線
( (2 つの内の上の方
) (
) ) が描けますね.これが,初期条件
x(0)
5.5
=
をみたす,連立微分方程式 (1) の特殊
y(0)
6.5
(
)
x(t)
解の解曲線です.t が 0 から正の方向へ動く時に,解
y(t)
がどちらの方向へ動くかが赤い矢印の向きで示されています
ね.時間がどんどん経つと (つまり t → ∞ の時),解は左上
の方向に行くこともわかります.
他にも色々な点をクリックしてみてください.クリックする
点 (初期条件) が違うと,時間がどんどん経った時の解の振
る舞い (左上に行くとか,右下に行くとか) が全く変わって
しまうということも分かります.
-4
0
4
8
y
8
4
0
-4
-8
-8
-4
0
4
8
なお,上では描画する範囲が自動的に決まっていましたが,これをオプションで指定することも可能です.例えば,同
じベクトル場を 0 ≦ x ≦ 1, 0 ≦ y ≦ 3 の範囲で描きたければ
maxima (%i9) plotdf([-y,-x],[x,0,1],[y,0,3]);
とします.
2015 年度 数式処理実習 [第 12 回] 連立微分方程式とベクトル場
July 9, 2015
[p.4/4]
描いたベクトル場の図を保存するには,ディスクマークのボタンをクリックし
ます.
Postscript filename の所に保存先とファイル名を指定します.
ファイル名の拡張子は.ps にして (例えば 12_ex1.ps など),su-shiki ディレ
クトリ内に保存します.(なお,デフォルトではホームディレクトリに保存さ
れてしまうので注意しましょう.)
解曲線が描かれていれば,それも一緒に保存されます.
(1) 以外の連立微分方程式でも,上と同じ考察ができます.一般的な形で述べましょう.f (x, y) 及び g(x, y) を既知の
2 変数関数とし,x(t), y(t) を未知関数とする次の連立微分方程式を考えます.(例えば (1) の場合なら,f (x, y) = −y ,
g(x, y) = −x です.)


 dx(t) = f (x(t), y(t))
dt
(2)

 dy(t) = g(x(t), y(t))
dt
(
)
(
)
(
) (
)
x(t)
a
x(t0 )
a
時刻 t0 において,(2) の解曲線
が点
を通るとします.つまり,
=
です.すると,
y(t)
b
y(t0 )
b
(2) が成り立つことから,この点における曲線の接ベクトルは


dx(t) (
) (
)
 dt t=t 
f (x(t0 ), y(t0 ))
f (a, b)

0 =
=
 dy(t) 
g(x(t0 ), y(t0 ))
g(a, b)
dt t=t0
(
となります.一言でいうと「解曲線が点
ことです.
a
b
)
(
を通るなら,その点における接ベクトルは
f (a, b)
g(a, b)
)
となる」という
)
(
)
x
f (x, y)
そこで,xy 平面上の各点
に対して,この点を始点とするベクトル
を対応させるようなベクトル
y
g(x, y)
場を描いておきます.このベクトル場は,plotdf([f(x,y),g(x,y)]) とすれば描くことができます.上と同様に,でき
た図の上の点をマウスでクリックすれば,その点を初期点 (初期条件) とする解曲線を赤い線で描いてくれます.
3
(
レポート課題
• http://www.cis.fukuoka-u.ac.jp/˜kawakubo/ の今日の所の sm140000 12.pdf の設問に答えよ.
• sm140000 12.tex (sm140000 12.pdf の原稿ファイル) をダウンロードし,˜/su-shiki ディレクトリに保存せよ.
• ファイル名「sm140000 12.tex」の「sm140000」の部分は各自の学籍番号に変更すること.また,Emacs でファイ
ルの中身を編集し,上段の氏名欄に自分の学籍番号と氏名を書くこと.
• Maxima で (手計算と書いてある部分は手計算で) 問題を解け.解答は原稿ファイルを編集して LATEX で作成せよ.
• PDF ファイルを作成し,印刷したものを提出すること.締切は 7 月 16 日.
• 疑問点等があれば,書いて下さい.