計算数理A演習第7・8回

第７・８回（第４週）

ルンゲ・クッタ法

前回まではオイラー法を用いた数値解法の演習を行いました。
今回は計算数学演習でも学んだルンゲ・クッタ法を用いて、次のような常微分方程式の初期値問題を数値的に解き、結果をグラフにしてもらいます。

(1)   \begin{align*}\dot{y}&= f(t, y) \\ y(0) &= a\end{align*}

プログラミング言語は違えど、アルゴリズムは同じですので、ルンゲ・クッタ法についての説明は省略します。計算数学演習を履修していない、あるいは、履修したはずだけど忘れてしまった、という人は、計算数学演習のページを確認するか、このページの最後の方に付録として載せておきますので、勉強しておいてください。

以下の課題では、4次のルンゲ・クッタ法で実装して下さい。

課題１

次の初期値問題を0\le t \le 5の範囲でh=0.5とし、実際にオイラー法とルンゲ･クッタ法を用いて解き、結果を比較せよ。

(2)   \begin{align*}\dot{y} &= \lambda y \\y(0) &= 1\\\lambda &= 1\end{align*}

課題２

ルンゲ・クッタ法を用いて、

(3)   \begin{eqnarray*}\ddot{y}+2\gamma\dot{y}+\omega^2y=0\\y(0)=1\\\dot{y}(0)=0\end{eqnarray*}


について、前回の課題8と同等のアニメーションを表示するプログラムを作成せよ。
ただし、\gamma=0, \omega=1, h=0.1とすること。また、オイラー法による解も同時にアニメーションで描画し、違いを確認すること。

振動現象について

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

減衰振動再考

前回行った減衰振動についてもう一度考えてみましょう。下の図のように、壁に一端を固定されたバネの他端におもりがつけられている状況を考えます。

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

(4)   \begin{equation*}\left\{ \begin{array}{rcl}\ddot{y} &=& -y(t) \\y(0) &=& 1 \\\dot{y}(0) &=& 0\end{array}\right.\end{equation*}

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

(5)   \begin{equation*}\left\{ \begin{array}{rcl}\ddot{y}(t) + 2k\dot{y}(t)+\omega^2y(t) &=& 0 \\y(0) &=& 1 \\\dot{y}(0) &=& 0\end{array}\right.\end{equation*}

ここでk, \omegaはそれぞれ抵抗力、復元力の強さを表す正の定数です。
式(5)の問題を連立１階常微分方程式に変換し、ルンゲ・クッタ法（できなければオイラー法）を用いて解き可視化（グラフ）してください。

課題3

上記の減衰振動に関して、以下の３通りについて、それぞれグラフを描け。
（例えば、\omegaをある値に固定して、kのみを変化させ、下記３通りを満たすようなkを選ぶとよい。３通りのグラフを色を変えて重ね描きするとよい。）

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

小ネタ：ドアのダンパ

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

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

強制振動

ここからは強制振動系を考えます（計算数理Aのmoodleを参照）。
先の減衰振動系に外から周期的な外力を加える事を考えます。扱う方程式は、以下のようになります。

(6)   \begin{equation*}\ddot{y} + 2k\dot{y}+\omega^2y=A\cos \omega_e t\end{equation*}

見ての通り、右辺が０であれば式(5)の減衰振動系となります。つまりは、式(4)の減衰振動系（放っておけば減衰して静止する）に対して、外から周期的な外力を加えているのが式(6)です。

課題4

式(6)を解くプログラムを作成せよ。下の図のようにy(t)とA\cos \omega_e t同時に表示せよ。例えばパラメータ、初期値としては、y(0)=1, \dot{y}(0)=0, k=0.1, \omega=1, \omega_e=1.3とし、時間刻み幅hはh=0.1とする。解y(t)とA\cos \omega_e tを同じグラフに表示することにする。

A=0の場合（つまり外力が無い場合）には、

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

A=0.1の場合には、

となり、過渡応答、定常応答が見られる。

A=0.1, \omega_e=1.0の場合には、次のようになる。  

課題5

課題4について、横軸をyとして、質点の動きをアニメーションで示せ。また、外力の項（A\cos \omega_e t）の力の向きと大きさも可視化せよ。

参考用動画(k=0.1, \omega=1.0, \omega_e=1.0, A=0.1）

動画中では、青いバーは\cos \omega_e tを表しており、右側にあるときは右向きに、左側にあるときは左向きに外力がかかっている。
（Aの大きさを青いバーの長さに反映させるかどうかは、各自の判断に任せます。）

課題6

(1)\omega_e =0.2, 0.4, 0.6, \cdots 2.0に対して、y(t)とA \cos \omega_e tの様子（振幅や位相の関係）を見てみよ。

(2)各\omega_eに対して、A=0.1, 1.0, 10.0とした場合のy(t)とA \cos \omega_e tの様子（振幅や位相の関係）を見てみよ。

(3)横軸y、縦軸\dot{y}として表示し、振動の様子を見よ。

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

今回の強制振動については、実は解析的に解（解析解）を求めることができます。課題6が合っているかは、解析解を求め、それと比較してみてください。
（実は解析解がわかってしまえば、下のレポート問題のパラメータはすぐにわかってしまう…）

レポート課題

強制振動の式(6)を用いて、課題5のようなアニメーション
（t=0から100までのアニメーション、時刻tもwindow内に表示）を表示するpythonプログラム作成せよ。ただし、以下の仕様を満たすこと。

【仕様1】
t\ge70のときは、
・t=70から現在の時刻tまでの範囲におけるyの最小値、最大値の位置をそれぞれ縦線で図示する。
（つまりt=80のときは70\le t\le 80の範囲におけるyの最小値、最大値の位置を縦線で図示する。）

【仕様2】
アニメーション終了後、70\le t\le 100におけるy(t)の最大値と最小値の差をDとして、そのDの値を画面に表示すること。
ただし、D>1となるように\omega_eを設定し、その値はアニメーションのウィンドウ内に表示すること。他のパラメータ（k, \omega）と初期値は課題4と同じものとし、A=0.1とする。

注意点：

  • プログラムは関数を用いてアルゴリズムが分かりやすくすること。
  • 質点の運動が画面内に収まるように座標系を設定すること。
  • \cos \omega_e tを表すバーは、見やすさを優先して大きさを調整しても構わない）
  • オイラー法の場合は１割程度減点する。
  • プログラムを〇〇.py形式（名前は自由）で作成し、moodleから提出すること。
  • インデント（時下げ）、コメントは必ず行うこと。
  • プログラムの最初の方にコメントとして学生番号と氏名は（最低限）必ず記述すること。
    （# bxxxxxx 氏名）

提出期限2024年5月9日（木）5月11日（土）0時00分

課題7

バネは大きく伸び縮みすると急激に硬くなるような場合もある。通常のバネが線形バネと呼ばれるのに対して、このようなバネは非線形バネと呼ばれるが、非線形バネに対する強制振動系として、次の微分方程式

(7)   \begin{align*}\dot{y} &= v \\\dot{v} &= -0.3 v - y^3+B \cos(t)\end{align*}

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

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

自然界の振動

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

課題8

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

(8)   \begin{align*}\dot{x} &= 1 - (1+b)x + x^2 y \\\dot{y} &= bx - x^2y\end{align*}

※ bの値が2より大きいか小さいかで振る舞いが変わります。

この方程式は、ブラッセレーター方程式といわれるもので、BZ反応と呼ばれる、時間的に化学成分の濃度が振動するような化学反応の振動をモデル化した方程式です。(BZ 反応とは何でしょう？その存在は古くから知られているものの、最近でも当時の高校生たちが新しい発見をしています。調べてみましょう。)

課題9

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

(9)   \begin{align*}\dot{x} &= 3x - xy \\\dot{y} &= -y + xy\end{align*}

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

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

課題10

惑星の運動を考える。万有引力の法則から２物体（質量はそれぞれ M、m）に働く力は

    \[\frac{GMm}{r^2}\]

となる。rは２物体間の距離、Gは万有引力定数で

    \[G \approx 6.674 \times 10^{-11} [m^3 / (kg\cdot s^2)]\]

である。原点に太陽（質量 M\approx 1.989 \times 10^{30} [kg] ）があり、惑星（質量 m）の位置を

    \[\vec{r}(t) = \left( \begin{array}{c} r_1(t) \\ r_2(t) \\ r_3(t) \end{array} \right)\]

とする。惑星に働く力は

    \[\vec{F} = - \frac{GMm}{|\vec{r}(t)|^2}\cdot \frac{\vec{r}(t)}{|\vec{r}(t)|} = -GMm\frac{\vec{r}(t)}{|\vec{r}(t)|^3}\]

となり、

    \[\vec{F} = m\vec{a} = m\ddot{\vec{r}}\]

よりg = GMとおくと、２階の常微分方程式の初期値問題

(10)   \begin{align*}\ddot{\vec{r}}(t) &= -g \frac{\vec{r}(t)}{|\vec{r}(t)|^3}\\\vec{r}(0) &= \vec{r}_0\\\dot{\vec{r}}(0) &= \vec{v}_0\end{align*}


が得られる。これを１階の常微分方程式系に直し、ルンゲ・クッタ法を用いて数値的に解き、惑星が太陽に衝突せずに周回軌道を取れるような初期条件を求めよ。

メモ：

実際にはGM\approx 1.3275\times 10^{20} [m^3 / s^2] ではありますが、簡単のためにGM=1、初期時刻での|\vec{r}(0)|は1から10程度の値を用いて良いです。実際の値を使うのであれば、地球の太陽からの距離が1.4960\times 10^{8} [km]、公転周期が約365日であることを考えて、

(11)   \begin{equation*}\left. \begin{array}{rcl}GM&\approx& 9.9097\times 10^{-4}[(10^8km)^3/Day^2]\\&\approx& 1.3202\times 10^{2}[(10^8km)^3/Year^2]\end{array}\right.\end{equation*}


あたりを使って計算すると良いでしょう。

難しいようであれば、まずは r_3 を考えずに、r_1–r_2 平面上の運動方程式として考えていくのも良いでしょう。

また、3次元のプロットは教えていません。r_1–r_2平面・r_2–r_3平面・r_1–r_3平面を並べるなどで断面図を描きましょう。どうしても3次元のプロットをしたい場合は、各自で調べてみてください。

課題11

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

    \[\ddot{\vec{r}}(t) = -g\frac{\vec{r}(t)}{|\vec{r}(t)|^3} + h\frac{\vec{r}(t)}{|\vec{r}(t)|^4}\]

課題5と同様に数値的に解をもとめよ。ただし、h \ll gとする。

付録：ルンゲ・クッタ法

1階常微分方程式のルンゲ・クッタ法

オイラー法は単純でそれはそれでよいのですが、実際の解に近い数値解を出すためには、時間刻み幅hをかなり小さくとらねばならないことが課題の結果からわかりました。一方、時間刻み幅を小さくしすぎると、計算コストがかかる（計算を大量にする必要になる）ことが問題です。そこで、より精度の高い解法が必要となります。そのような問題点の解法として、ルンゲ・クッタ法を扱い、それを用いた数値解法の演習を行います。

この演習ではせいぜい3、4変数の連立微分方程式を解くので、計算コストについてはそこまで問題になりませんが、例えば、生体分子のシミュレーションでは1,000や10,000のオーダーの変数を扱います。そうすると計算コストをきちんと考えてシミュレーションしないと、シミュレーション終了まで月・年単位でかかってしまうこともあります。そのようなことを防ぐため、様々な解法やクラスターマシンやスパコンといった大規模な計算システムを有効に活用して計算を行います。

ルンゲ・クッタ法の導出については、このページの最後を参照してください（詳細は参考書を参照してください）。
ルンゲ・クッタ法は次に示すアルゴリズムで、オイラー法に比べれば多少複雑なものとなっています。
比較の為に両方のアルゴリズムを載せておきます（赤字の部分が異なります）。
よく比較してください。

オイラー法 ルンゲ・クッタ法
T ：解を求める時間の上限
h ：時間刻み幅
N ：時間刻み数
y：変数
a：初期値
y_new：実数型の変数		 T ：解を求める時間の上限
h ：時間刻み幅
N ：時間刻み数
y：変数
a：初期値
y_new, r_1, r_2, r_3, r_4 ：実数型の変数

T, N, h を定める
a の入力を受ける(input関数とfloat)
y = a

i = 0, 1, 2, … , Nの順に
t = i*h

y_new = y + h*f(t, y)

 

y = y_new

一定間隔毎にグラフを描画

(ax.set_data関数とplt.pause関数)
を繰り返す

f(t, y) は関数として定義するとよい.

T, N, h を定める
a の入力を受ける (input関数とfloat)
y = a

i = 0, 1, 2, … , Nの順に
t = i*h

r_1 = f(t, y)
r_2 = f(t + h/2, y + (h/2)*r_1)
r_3 = f(t + h/2, y + (h/2)*r_2)
r_4 = f(t + h, y + h*r_3)
y_new = y + (h/6)*(r_1 + 2*r_2 + 2*r_3 + r_4)

y = y_new

一定間隔毎にグラフを描画

(ax.set_data関数とplt.pause関数)


を繰り返す

f(t, y) は関数として定義するとよい.

＊ オイラー法では, （とびとびの）時刻 t = jh での値 y から 次の（とびとびの）時刻 t = (j+1)h での値 y_new を直接計算しています。
本来、微分方程式では時刻は連続ですが、この場合は時刻 t = jh ～ (j+1)h の間の事は一切考慮されておらず、その分誤差が現れてきます。

＊ ルンゲ･クッタ法では, この t = jh ～ (j+1)h の間の情報を (r_1, ～ r_4 を使って ) いくらか取り入れる事で、よりよい近似計算が可能になっています。
勿論誤差は現れるのですが、オイラー法と比べて雲泥の差です。
r_1〜r_4を求めることで処理は増えますが、それと比較しても同じ計算精度を実現するときの計算回数は非常に少なくなります（上の課題1、2あたりで同じ計算精度を実比較するときの時間刻み数Nを比較してみてください）。

2階（n階）常微分方程式のルンゲ・クッタ法

前回の課題で扱った

(12)   \begin{align*}\ddot{y}(t) &= -\omega^2 y(t) \\y(0) &= 1 \\\dot{y}(0) &= 0\end{align*}


をy_1(t)=y(t), y_2(t)=\dot{y}(t)とおいて、次のように連立の１階の常微分方程式に帰着したもの

(13)   \begin{align*}\dot{y}_1(t) &= y_2(t) \\\dot{y}_2(t) &= -\omega^2 y_1(t) \\y_1(0) &= 1 \\y_2(0) &= 0\end{align*}

を再び考えます。
このような連立の常微分方程式をルンゲ･クッタ法を用いて解く場合のアルゴリズムを次に示します。

2階常微分方程式のルンゲ・クッタ法
T ：解を求める時間の上限
h ：時間刻み幅
N ：時間刻み数
Y1, Y2：変数
a1, a2：初期値
Y1_new, Y2_new, r1_1, r1_2, r1_3, r1_4, r2_1, r2_2, r2_3, r2_4 ：実数型の変数

T, N, h を定める
a1, a2 の入力を受ける(input関数とfloat)
y1 = a1
y2 = a2


i = 0, 1, 2, …. ,Nの順に

t = i*h


r1_1 = f1(t, y1, y2)
r2_1 = f2(t, y1, y2)

r1_2 = f1(t + h/2, y1 + (h/2)*r1_1, y2 + (h/2)*r2_1)
r2_2 = f2(t + h/2, y1 + (h/2)*r1_1, y2 + (h/2)*r2_1)

r1_3 = f1(t + h/2, y1 + (h/2)*r1_2, y2 + (h/2)*r2_2)
r2_3 = f2(t + h/2, y1 + (h/2)*r1_2, y2 + (h/2)*r2_2)

r1_4 = f1(t + h, y1 + h*r1_3, y2 + h*r2_3)
r2_4 = f2(t + h, y1 + h*r1_3, y2 + h*r2_3)

y1_new = y1 + (h/6)*(r1_1 + 2*r1_2 + 2*r1_3 + r1_4)
y2_new = y2 + (h/6)*(r2_1 + 2*r2_2 + 2*r2_3 + r2_4)

y1 = y1_new
y2 = y2_new

一定間隔毎にグラフを描画 (lines.set_data関数とplt.pause関数)

を繰り返す

f1(t, y1, y2), f2(t, y1, y2) は関数として定義するとよい.

また、f1(t, y1, y2) および f2(t, y1, y2) はそれぞれ実数型の関数です。
（12）式の場合には f1 が第一式の右辺, f2 が第二式の右辺となります。

(今回, ω = 1 とします.)

小ネタ：変数名の付け方

言語によっては、意外な変数名がすでに使用されていることがあります。例えばC言語では、y0, y1, ynはグローバル変数の名前として好ましくありません（調べてみましょう）。y_0, y_1, y_nや、Y0, Y1, Ynで代用しましょう。また、すでに気づいている人もいるかもしれませんが、pythonでは「lambda」は無名関数で使うため、使えません。「Lambda」などで代用しましょう。

ルンゲ・クッタ法の仕組み（2段公式を例に）

簡単にルンゲ･クッタ法について説明します。演習で取り扱ったルンゲ･クッタ法は４段公式とよばれるものですが、ここでは簡単の為、２段公式を取り扱います。

次の常微分方程式の初期値問題の数値解を求めたいとしましょう。

(R1)   \begin{align*}\dot{x} &= f(t, x) \\x(t_0) &= x_0\end{align*}


２段公式とは、t=t_nで（R1）の近似解x_nが求められたとして、時間ステップ幅hだけ進んだt_{n+1}=t_n + hにおける値x_{n+1}を次の形で計算する公式です。

(R2)   \begin{align*}x_{n+1} &= x_n + (\gamma_1 k_1 + \gamma_2 k_2) \\k_1 &= hf(t_n, x_n)\\k_2 &= hf(t_n+\alpha h, x_n+\beta k_1)\end{align*}


つまり、x_nからx_{n+1}へ進む１ステップにおいて方程式（R1）の第一式の右辺のf(t, x)を２回計算するので、２段公式と呼びます。
この２段公式には、パラメータ\alpha, \beta, \gamma_1, \gamma_2が含まれますが、これらは方程式（R1）の真の解x(t)のテイラー展開

(R3)   \begin{align*}x(t_n+h)=&x(t_n)+hf(t_n,x(t_n)) \\&+ \frac{1}{2!}h^2\{f_t(t_n, x(t_n)) + f_x(t_n, x(t_n))f(t_n, x(t_n))\} + O(h^3)\end{align*}


と（R2）式の第一式の同様の展開とがhのなるべく高い次数の項まで一致するように定めます。
つまり、x_n=x(t_n)と仮定して（R2）式のk_2の式の右辺をf(t, x)周りでテイラー展開すると、

(R4)   \begin{align*}x_{n+1}=&x_n+h(\gamma_1 + \gamma_2)f(t_n,x_n) \\&+ h^2\gamma_2\{\alpha f_t(t_n, x_n) + \beta f_x(t_n, x_n)f(t_n, x_n)\} + O(h^3)\end{align*}


となります。
これと（R3）がh^2の項まで一致するように

(R5)   \begin{align*}\gamma_1 + \gamma_2 = 1\\\gamma_2 \alpha =\gamma_2 \beta = \frac{1}{2}\end{align*}


とおく。公式（R2）で定められる近似解のテイラー展開がhの２次の項まで真の解のテイラー展開に一致するという意味で、この公式を２次の公式と呼びます。
この公式は２段公式でもあるので、２段２次の公式とも呼ば れます。パラメータを定める方程式（R5）には自由度が２残っているのでその定め方によって色々な公式が導かれます。（改良オイラー法、ホイン法、修正オイラー法など）
演習で用いた４段４次のルンゲ･クッタ法の公式も同様にして求めることが出来ます。

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University