第8回(2009年6月17日)


今日の演習


前回までに,常微分方程式はオイラー法,ルンゲ・クッタ法の二種類の解法を用いて解く事ができました.h が十分小さければルンゲ・クッタ法を用いて得た数値解(近似解)は,真の解に十分近いといえるでしょう.今回からは,ルンゲ・クッタ法を常微分方程式を解く信頼できる方法として採用し,具体的な問題を解き,得られた解を見るだけでなく,現象に立ち返りその意味を再考してみましょう.(ルンゲ・クッタ法になじめない方はオイラー法を用いるのでも構いません. )


減衰振動再考

減衰振動についてもう一度考えてみましょう.

課題1

下の図のように、壁に一端を固定されたバネの他端におもりがつけられている状況を考えます。

バネとおもり

このおもりを少しだけ右に引っ張ってから手をはなす。おもりの中心の時刻 t におけるつり合いの位置からのずれを y(t) とします。また、時刻 0 でのおもりの位置を y(0)=1、おもりの速度を y'(0)=0 とするとつぎの単振動の問題となります。(フックの法則を用いた)

  (1)

これではおもりに働く力はバネの復元力のみで、ある意味非現実的です。今回は復元力以外におもりの速度 y' に比例する抵抗が働く場合を考えましょう(講義ではダンパを加えた場合として説明されている)。改良されたモデル方程式は次のようになります。

  (2)

ここで k,ω はそれぞれ抵抗力、復元力の強さを表す正の定数です。(2)式の問題を連立1階常微分方程式に変換し、ルンゲ・クッタ法(オイラー法)を用いて解き、GLSCを用いて解を可視化(グラフ)してください。(前回, もしくは前々回のレポートで提出されているプログラムを使い回してください. それを元に今日は現象を観察してもらいます)

以下の3通りについて,それぞれグラフを描いてください.(例えば,ωをある値に固定して,k のみを変化させ,下記3通りを満たすような k を選ぶとよい.3通りのグラフを色を変えて重ね描きするとよい.)

  1. k < ω の場合:減衰振動
  2. k > ω の場合:過減衰
  3. k = ω の場合:臨界減衰

ところで,ドアは一種のバネ振り子+ダンパ系と見る事ができる.設計する際,ダンパの強さを臨界減衰になるよう設計するのが一番良いとされる.なぜであろうか?考えてみてください.(良いドアである為には次の二つが重要である.素早くしまる事,しまるときの音が小さい事.)

自動車やバイクのサスペンションも一種のバネ振り子+ダンパ系であるが,おそらくはやはり臨界減衰に近いような設定が良いのであろうと考えられます,(実際にはより複雑な状況もあり,そこまで単純ではないかも知れませんが.)


強制振動

ここからは強制振動系を考えよう.先の減衰振動系に外から周期的な外力を加える事を考える.


扱う方程式は,以下のようになる.

  (3)

見ての通り,右辺が0であれば(2)の減衰振動系となります.つまりは,(2)の減衰振動系(放っておけば減衰して静止する)に対して,外から周期的な外力を加えているのが(3)であります.


課題2

(3)を解くプログラムを作成せよ.下の図のように y(t) と Acosωet 同時に表示するようにしてください. 例えばパラメータ,初期値としては,

とし,時間刻み幅 h は h=0.1 とする.解 y(t) を赤色,外力 Acosωet を緑色で表示することにする.

A = 0 の場合(つまり外力が無い場合)には,

となり,当然ながら減衰振動の様子がみられる.

A = 0.1 の場合には,

となり,過渡応答,定常応答が見られる.

A = 0.1, ωe = 1.0 の場合には,次のようになる.


課題3

(1) ωe = 0.2, 0.4, 0.6, ...., 2.0 に対して, y(t) と Acosωet の様子(振幅や位相の関係)を見てみよ.

(2) 各 ωe に対して, A = 0.1, 1.0, 10.0 とした場合の y(t) と Acosωet の様子(振幅や位相の関係)を見てみよ.

(3) 横軸 y1, 縦軸 y2 として表示し, 振動の様子を見よ.

(1)(2) から分かるように, バネの固有の振動数 ω (外力, 抵抗が無いときの振動の周期の逆数)と, 外力の振動数 ωe が近くなるにつれ, 大きな振幅の振動が現れます. これが「共鳴(共振)」と言われる現象で, 例えばブランコを手で押して大きく揺らす事等は(ふつうブランコの揺れに合わせて押すことで振幅を大きくすると思いますが.), 強制振動による共鳴の典型例です. (逆にこのようなバネ系は, 強く揺すってもまあそれ相応にしか振幅はふえません. 外からの揺すりの「タイミング」が肝ですね. それを実感してください.)


課題4

バネは大きく伸び縮みすると急激に硬くなるような場合もあります. そのようなバネに対する強制振動系として, 次の微分方程式

dy/dt = v

dv/dt = - 0.3*v - y3 + B*cos(t)

を考えます. この式を数値的に解いて、B = 0, 1, 2, ...., 12, ..., 60 ... 等いろいろな場合の y の時間変化を調べてください.特に y, v の初期値が少し違う場合の. y の時間変化を比較してください. 課題2,3とどのような特徴の違いがあるでしょうか. (十分大きい t までの結果を見てみてください. また横軸 y, 縦軸 v として表示した場合の様子も見てみてください.)

ハードディスクから原子炉まで、様々な機械の「振動」が、このような方程式を応用しながら解析されています.


自然界の振動

上記では, 摩擦が働く等して放っておいたら止まってしまうようなシステム (バネ系) が, 外部からの「振動(時間に依存して振動的に変化する外力)」を受ける事で定常的に振動する様子について考察しました. しかし実際には, 外からの影響が「振動的」でなくても(外力が時間変化しない. つまり方程式の右辺に「 t 」が現れない.), システム自身が振動的になる場合があります. そのような振動は「自励振動」と呼ばれます.以下, 自励振動の簡単な例題を見ていきます. 生物を含む自然界の様々な場面で見られるいろいろな「振動(リズム)」の多くが, 自励振動であります.(体内時計, 脳波, 生態系, 為替変動, 蛇口の水滴, 地震, 気象,.. 等等. *宇宙や惑星の運動等は少々違います...)


課題5

次の連立微分方程式を適当な初期条件のもと、オイラー法かルンゲ・クッタ法で解いてください。

dx/dt = 1 - (1 + b)*x + x2y

dy/dt = b*x - x2y

(1) b < 2 の場合
(2) b > 2 の場合

この方程式は、ブラッセレーター方程式といわれるもので、BZ反応と呼ばれる、時間的に化学成分の濃度が振動するような化学反応の振動をモデル化した方程式です。(BZ 反応とは何でしょう? http://www.chem.scphys.kyoto-u.ac.jp/nonnonWWW/kitahata/bz_1.html あたりが参考になります。)

最近では例えば皆さんの体の細胞が刻む体内時計の概日リズム(おおよそ25時間周期のリズム)の仕組みは, このブラッセレーター方程式の示す階と同等の数理構造から理解できると考えられています.


課題6

次の連立微分方程式を適当な初期条件のもと、オイラー法かルンゲ・クッタ法で解いてください。

dx/dt = 3x - xy

dy/dt = -y + xy

課題5の方程式同様、この方程式も振動解を見せますが、その特徴は少々違っております。例えば初期条件を変えた場合に、これらの(課題5と課題6それぞれの)方程式系の解の振る舞いが、どのように変化していくか、見比べてください。

この方程式は、ロトカ・ボルテラ方程式といわれるもので、生態系内のある被食者と補食者の数 (x, y) 増減関係を簡単にモデル化した方程式です。(各時刻、被食者はある割合(3x)で増えて、補食者に出会うと食べられ減る(-xy)。同様に捕食者は被食者に出会うとそれを食べて増えて(xy)、ある割合で(-y)で寿命等で減っていく。)

自然界の関係は常にこのように変動し続けています.


計算数学演習の期末レポートで扱った「ローレンツ方程式系」も, このような自励振動の典型例であります. (おさらいしてみてください.)
宇宙や惑星の運動は自励振動とは異なります. こちらは「摩擦がない」ので振動(回転)し続けます.でも自励振動と同じくらい複雑な振動(回転)が起こり得ます.(このことは銀河系等に様々な形があることの原因でもあります.)

課題7

 惑星の運動を考えます。万有引力の法則から2物体(質量はそれぞれ M,m)に働く力は

wpe3.jpg (1439 バイト)

となります。rは2物体間の距離、Gは万有引力定数で

wpe4.jpg (2817 バイト)

です。原点に太陽(質量 M)があり、惑星(質量 m)の位置を

wpe1.jpg (3009 バイト)

とします。惑星に働く力は

wpe6.jpg (5166 バイト)

となり、

wpe7.jpg (1519 バイト)

より

wpe8.jpg (1269 バイト)

とおくと、2階の常微分方程式の初期値問題

wpe9.jpg (4810 バイト)

が得られます。これを1階の常微分方程式系に直し、数値計算にて解をもとめて下さい。
(まずは r3 を考えずに, r1 - r2 平面上の運動方程式として考えてもよろしいです.)


課題8

水星は太陽に非常に近く、その効果が公転に影響するため、その運動は以下のような方程式に従うようになる。

課題7同様、数値的に解をもとめてください。ただし h << g とします。


次回の予告

x軸上にバネを並べ、おもりでつなげた場合どのような運動をするか、数値的に計算してみましょう。ただし、全てのバネの自然長R(伸びも縮みもしていないときの長さ)とバネ定数が同じで、おもりの質量も同じとします。


左のおもりから順に 0, 1, ..., M, M+1 と番号を振ります。全てのバネが伸びも縮みもしていない状態であれば、m番目のおもりの位置はバネの自然長Rを用いて、m R と書けます。これを用いると、m番目(m≠0,m≠M+1)のおもりの運動は、m R からのズレ ym は次のように書けます。



・・・ (1)

0 番目と M+1 番目のおもりは固定されて動かないとします。つまり、



          ・・・ (2)

とします。

補足: 先週までのバネの伸びの方程式から式(1)がでてきます。yjをおもりの場所そのものとしても式(1)と同じ方程式になることと合わせて確認してみてください。

課題9


M=3、ω=1、R=1 として、式(1)と式(2)を解くプログラムを作成せよ。ただし初期条件は、y1=0.25、y'1=0.0 とする。適度なステップ毎に、おもりの位置(m番目のおもりの位置は (m R + ym, 0))を g_marker で表示させること。次図は M = 3 の場合で, グラフにする横軸の範囲は[0,2]とした。





戻る