計算数理A演習第6回 – Masashi Shiraishi Homepage

今日の演習

今日は以下の内容で演習を行います.

2階常微分方程式の数値解法（オイラー法）
$n$ 階常微分方程式の数値解法（オイラー法）

１階常微分方程式の数値解法

前回，簡単な１階常微分方程式をオイラー法にて解いてもらいました．
いくつか講義で出てきた問題をあげておきます．
オイラー法で解いてみて，解をグラフとして表示してみてください．
初期値やパラメータを適当にきめて，いくつかシミュレーションしてみてください．

(1) $\begin{equation*} \left\{ \begin{array}{rcl} \frac{dx}{dt} &=& \gamma x \\ x(0) &=& x_0 \end{array} \right. \end{equation*}$

(2) $\begin{equation*} \left\{ \begin{array}{rcl} \frac{du}{dt} &=& u(1-u) \\ u(0) &=& u_0 > 0 \end{array} \right. \end{equation*}$

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

２階常微分方程式もオイラー法を用いて解くことができます。
次の問題を考えて行きましょう。

(3) $\begin{equation*} \left\{ \begin{array}{rcl} y'' &=& -\omega^2 y \\ y(0) &=& 1 \\ y'(0) &=& 0 \end{array} \right. \end{equation*}$

講義の説明の通り(下の枠中)、これは単振動を記述する２階の常微分方程式です。
ただし簡単の為、時間 $t$ での微分を ‘ を用いてあらわしています。

バネ定数 $k$ のバネの一端を固定し、他端に質量 $m$ のおもりをつけた状況を考える。

このおもりを少しだけ右に引っ張ってから手をはなす。
おもりの中心の時刻 $t$ におけるつり合いの位置からのずれを $y(t)$ とすると。
$y(t)$ が満たす微分方程式は

(4) $\begin{equation*} \frac{d^2 y}{dt^2} = - \omega^2 y \end{equation*}$

ただし、

(5) $\begin{equation*} \omega = \sqrt{k/m} \end{equation*}$

とした。

２階の常微分方程式は次のように連立の１階の常微分方程式に帰着できます。
$y_1(t)=y(t), y_2(t)=y'(t)$ とおくと、式(3)は、

(6) $\begin{equation*} \left\{ \begin{array}{rcl} y_1' &=& y_2 \\ y_2' &=& -\omega^2 y_1 \\ y_1(0) &=& 1 \\ y_2(0) &=& 0 \end{array} \right. \end{equation*}$

となります。
この形になれば、前回紹介したオイラー法を用いて解を求めることができます。
アルゴリズムは次のようになります。
但し、 $T$ を解を求める時間の上限、 $h$ を時間刻み幅、 $N$ を時間刻み数とします。
また、変数は $y_1, y_2$ として、 $a_1, a_2$ を初期値とします。
また、 $f_1(t, y_1, y_2)$ および $f_2(t, y_1, y_2)$ はそれぞれ関数で式(6)の場合には $f_1$ が第一式の右辺、 $f_2$ が第二式の右辺となります。

オイラー法のアルゴリズム

解を求める時刻の上限を T とし、適当な自然数 N を定めて h=T/N とする。

T, h, a1, a2, y1, y2, y1_new, y2_new, t 及び関数 f1, f2 は実数型(double)

N, j は整数型(int)

T,N,h を定める

y1 = a1, y2 = a2

t = 0 でのグラフの描画 (g_move 関数)

j = 0,1,2,….,N-1の順に

t = j*h

y1_new = y1 + h*f1(t, y1, y2)

y2_new = y2 + h*f2(t, y1, y2)

y1 = y1_new

y2 = y2_new

一定間隔毎にグラフを描画 (g_plot　関数)

を繰り返す

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

課題１

式(3)を解析的に解け。

課題２

$\omega=1$ として式(6) をオイラー法で解き、グラフを描け。グラフは横軸が $t$ であり、縦軸を $y(t)$ とせよ。
また、時間刻み幅 $h$ を変えたときに結果がどのようにかわるか観察せよ。
0≦t≦20の範囲において、h=0.0001,0.001, 0.01, 0.1それぞれの場合について、解析解とどの程度あっているか各自調べよ。

色々な h に対する計算結果は次のようになる（緑はh=0.1、青はh=0.01、ピンクはh=0.001、黄色はh=0.0001のときの結果。赤は解析解のグラフ）。
この結果から、h は十分小さい値でないとまずいことがよく分かる。

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

課題３

課題２のプログラムを用いて、

(7) $\begin{equation*} E = \frac{1}{2}\left( y' \right)^2 + \frac{1}{2} y^2 \end{equation*}$

を計算せよ。結果は、横軸を $t$ 、縦軸を $E$ とするグラフで確認すること。
この $E$ は $m=k=1$ のときの、バネの力学エネルギーに対応する。

課題４

先程単振動を考察しましたが, 実際のバネではバネの復元力以外におもりの速度 $y'$ に比例する抵抗が働く。
その抵抗を考慮すると、式(3)ではなく次の式になる。( $2\gamma y'$ の項が新たに加わった。)

(8) $\begin{equation*} \left\{ \begin{array}{rcl} y'' + 2\gamma y' + \omega^2 y &=& 0 \\ y(0) &=& 1 \\ y'(0) &=& 0 \end{array} \right. \end{equation*}$

ここで $\gamma$ は抵抗力の強さを表す正の定数です。
式(8)の問題を連立１階常微分方程式に変換し、オイラー法を用いて解き、GLSCをつかって、下のように「左端から線の繋がった赤い質点が振動する」動画を作成してください。

ただし、 $\omega = 2, \gamma = 0.4$ とし、時間の刻み幅 $dt$ は、0.01とします。
ます。また、左上の”t= ..”は現在の $t$ を表示している。
デモプログラムでは、時間の刻み幅は $dt=0.01$ として、 $t=50$ まで計算しています。
また、提出するプログラムでは、アニメーションの描画を $t$ が $0.1 ( =10dt)$ 動く毎におこなうようにしてください。（ヒント１： “%”をつかったif文）
これは、無駄にプログラムの実行時間が長くなるのを防ぐためです。
以下、アドバイス等：
“t=”の部分のアニメーションのさせ方：カウントダウンを行うサンプルプログラムを参考にしてください。
質点の赤い領域：g_area_color(G_RED)とg_circleで塗りつぶしにG_YES。

課題５

下図のような重力下で棒の長さ $L$ が一定の振り子運動について、振り子の角度(ラジアン）を $X$ とすると $X$ は

(9) $\begin{equation*} LX'' = -g \sing(X) \end{equation*}$

という運動方程式を満たす（ $g$ は重力加速度、摩擦は無視しています）。
$L=1, g=9.8$ として、この方程式を適当な初期条件、 $X(0) = 0, X'(0) = a$ でシミュレーションし、GLSC を用いた振り子の動画を作成せよ。

今回の課題は、オイラー法では十分に小さい $h$ を用いないと計算結果がすぐにずれてしまいますが, 練習としてやってみてください.

次回, オイラー法より計算精度の高いルンゲ・クッタ法を紹介しますが, そのときに比較してみてください.

課題６

課題４の振り子が２つ相互作用する場合を考える.
例えば, ２つの振り子の角度を $2\pi X_1, 2\pi X_2$ とすると、 $X_1, X_2$ が

(10) $\begin{equation*} \left\{ \begin{array}{rcl} X_1'' &=& - A \sin(2\pi X_1) - C \sin(2\pi(X_1 - X_2)) \\ X_2'' &=& - B \sin(2\pi X_2) - C \sin(2\pi(X_2 - X_1)) \end{array} \right. \end{equation*}$

という運動方程式を満たすとする。
この方程式を適当な $A, B, C$ と初期条件、 $X_1(0), X_1'(0), X_2(0), X_2'(0)$ でシミュレーションし、GLSC を用いた振り子のアニメーションを作成せよ.
また振り子の数をもっと多くした場合も考察せよ.