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

July 9, 2015
[p.1/4]
2015 年度数式処理実習 [第 12 回] 連立微分方程式とベクトル場
1
連立微分方程式
前回は，１つの未知関数に対する微分方程式を解きました．今回は，２つの未知関数に対する連立微分方程式 (“ 微分
方程式系 ”とも言います) を解きましょう．
２つの未知関数 x(t), y(t) に対する連立微分方程式


 dx(t) = −y(t)
dt
(1)

 dy(t) = −x(t)
dt
について考えます．未知関数が２つで，微分方程式も２つです．
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 日．
• 疑問点等があれば，書いて下さい．

Download Report