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 日. • 疑問点等があれば,書いて下さい.
© Copyright 2024 ExpyDoc