コンテンツ

オイラー法２

今日の演習

2階常微分方程式の数値解法（オイラー法）
n階常微分方程式についても連立1階常微分方程式に帰着できること知る

オイラー法による２階常微分方程式解法

２階常微分方程式をオイラー法を用いて解いてもらいます。次の問題を考えます。

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

これは単振動を記述する２階の常微分方程式です。
２階の常微分方程式は次のように連立の１階の常微分方程式に帰着できます。
$y_1(t) = y(t)$ , $y_2(t) = \dot{y}(t)$ とおくと、式(1)は

(2) $\begin{equation*} \left\{ \begin{array} $\dot{y}_1 = y_2$ \\ \dot{y}_2 = -y_1 \\ y_1(0) = 1, y_2(0) = 0 \end{array} \right. \end{equation*}$

となります。
こうなれば先程と同様、オイラー法を用いて解を求めることができます。

もう少し一般的な書き方をすれば

$\left\{ \begin{array} $\dot{y}_1 = f_1(t, y_1, y_2)$ \\ \dot{y}_2 = f_2(t, y_1, y_2) \end{array} \right.$

$f_1(t, y_1, y_2)$ および $f_2(t, y_1, y_2)$ は、それぞれ関数で式(2)の場合には $f_1$ が第一式の右辺、 $f_2$ が第二式の右辺となります。
今回は $f_1(t, y_1, y_2) = y_2$ 、 $f_2(t, y_1, y_2) = -y_1$ という形になります。
アルゴリズムは次のようになります。
ただし、T を解を求める時刻の上限とし、N を時間刻み数、h = T/(N-1) を時間刻み幅とします。
変数は Y1, Y2 として、a1, a2 を初期値とします。

オイラー法による2階常微分方程式のアルゴリズム

define文を使った定数の定義：T, h, a1, a2
関数の宣言：関数 f1, f2 は実数型(double)

変数の宣言  j, Y1, Y2, Y1_new, Y2_new, t（時刻）, N（時間刻み数）
N を定める
Y1 = a1, Y2 = a2

j = 0,1,....,N-1の順に
    t = j*h // 次の時刻
    一定間隔毎に t, Y1, Y2 を画面に表示（printf関数）
    Y1_new = Y1 + h*f1(t, Y1, Y2)
    Y2_new = Y2 + h*f2(t, Y1, Y2)
    Y1 = Y1_new
    Y2 = Y2_new    
を繰り返す

各ステップで, Y1_new, Y2_new が両方求まってから Y1 = Y1_new, Y2 = Y2_new と値を更新するようにします。
ここの順序が変わってくると, Y1 と Y2 で対応する時間がずれてしまうので, 正確な計算になりません。

課題１

(2) をオイラー法で解き、グラフを描け。グラフは横軸が t であり、縦軸を y(t) とせよ。
また、時間刻み幅 h を変えたときに結果がどのようにかわるか観察せよ。
時刻 0≦t≦20の範囲において、h = 0.0001, 0.001, 0.01, 0.1それぞれの場合(時間刻み数 N はTとhを使って直接調整)について、微分解 $y(t)=\cos t$ とどの程度あっているか調べよ。
色々な h に対する計算結果は次のようになる。

この課題の結果から分かるように、オイラー法では h を十分小さくとらないと本当の解をうまく近似できません。
h を小さく取るということは、計算時間が多くかかることを意味します。
そこでより高精度で大きな h でもうまく解を近似できるような解法が求められます。
そのような要求に応えるルンゲ・クッタ法があります（次回に説明します）。

課題２

次の微分方程式

(3) $\begin{equation*} \ddot{y} + 0.1\dot{y} + y = A\cos(Bt) \end{equation*}$

を、初期値 $y(0)=1, \dot{y}(0)=0, A=2.0$ のもと数値的に解き $B = 0,0.5,1,2$ の場合の $y$ の時間変化を調べて下さい(時刻の上限を少し大きめにして下さい)。
また、 $y$ と $\dot{y}$ の関係（gnuplot で横軸を $y$ , 縦軸を $\dot{y}$ として）も plot してみて下さい。
これは摩擦抵抗のあるバネ振動子に、周期的な外力を加えたもので、強制振動系と呼ばれます。
様々な機械の「振動」が、このような方程式を応用しながら解析されています。

レポート課題

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

このおもりを少しだけ右に引っ張ってから手をはなす。
おもりの中心の時刻 t におけるつり合いの位置からのずれを $y(t)$ とする。
また、時刻 0 でのおもりの位置を $y(0)=1$ 、おもりの速度を $\dot{y}(0)=0$ とすると、
課題１と全く同じ、つぎの単振動の問題となる。

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

これではおもりに働く力はバネの復元力のみで、ある意味非現実的である。
今回は復元力以外におもりの速度 $\dot{y}$ に比例する抵抗が働く場合を考える。
（おもりと床との間の摩擦力のようなものをイメージ）。
改良されたモデルは次のようになる（ $2k\dot{y}$ の項が新たに加わったもので、摩擦のようなものをあらわしている）。

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

ここで $k, \omega$ はそれぞれ抵抗力、復元力の強さを表す正の定数である。
課題は、（5）式の問題を連立１階常微分方程式に変換し、オイラー法を用いて解くプログラムを作成せよ。

ただし、 $\omega=1.0$ , $h=0.0001$ とし、出力は $t, y(t), \dot{y}(t)$ の３成分で、

0.00 $y$ (0.0) $\dot{y}$ (0.0)
0.01 $y$ (0.01) $\dot{y}$ (0.01)
0.02 $y$ (0.02) $\dot{y}$ (0.02)
…
49.99 $y$ (49.99) $\dot{y}$ (49.99)
50.0 $y$ (50.0) $\dot{y}$ (50.0)

というふうに、 $t$ は0.01刻みで、各列の間は「スペースを入れて」出力するようにすること。
（ヒント:100回に１回出力するので、%を使うと良いでしょう）

$k$ の値について： $y(t)$ のグラフが下図のようになるkの値を見つけ、その $k$ の値を使ってプログラムを作成せよ。
（ヒント：正解データと一致しているか確認）

提出するファイルは、見つけた $k$ の値が変数に代入されているCプログラム（そのままgnuplotで描画できるような形式に結果を出力する）および下記のようなグラフ（pdf形式、eps形式、png形式のいずれか）とする。
ただし、グラフのタイトルに学生番号を入れること。
Cプログラムのファイルの先頭にコメント文で学生番号と氏名を記載すること。
Cプログラムのファイル内のコメント文やインデントなどの読みやすさ、グラフの縦軸・横軸なども採点対象に含める。

提出締め切り：2025年1月29日(水)午後6時

n階常微分方程式のオイラー法

ここまでに、2階常微分方程式は2変数1階常微分方程式へと変換し、1変数１階常微分方程式のアルゴリズムと同様にオイラー法で数値計算できることを学びました。
この変換方法は簡単に一般化できて、n階微分されている変数に対してn-1階微分を新たな変数と置き換えることで、その変数の1階微分としてn階微分されている変数を置き換えていけば全ての高階微分を１階微分で表現できます。
結果として、n変数の1階微分方程式へと変換することができます。
例えば、次の様な一般的な場合を考えます。

(6) $\begin{equation*} y^{(n)} = f( x, y, \dot{y}, \cdots, y^{(n-1)}) \end{equation*}$

このn階常微分方程式に対して、 $z_i(x)=y^{(i)}$ とすると、

(7) $\begin{eqnarray*} \dot{z}_{1}(x) &=& z_2(x) \\ \dot{z}_{2}(x) &=& z_3(x) \\ &\vdots& \\ \dot{z}_{n-1} &=& y_{n} \\ \dot{z}_n &=& f(z, z_1, z_2, \cdots, z_n) \end{eqnarray*}$

となりn変数1階常微分方程式とできます。

Masashi FUJII

計算数学演習第１２回