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

今日の演習

連結バネ系

連結バネ系

前回、前々回と偏微分方程式の数値解法について演習を行いましたが、今回は今一度オイラー法とルンゲクッタ法を用いた数値解法に取り組みます。
今回は、次のように x 軸上にバネを並べ、おもりでつなげた場合どのような運動をするか、数値的に計算してみましょう。
ただし、全てのバネの自然長 R (伸びも縮みもしていないときの長さ)とバネ定数が同じで、おもりの質量も同じとします。

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

(1) $\begin{equation*} \left\{ \begin{array}{l} y''_m = -\omega^2(y_m - y_{m-1}) - \omega^2(y_m - y_{m+1}) = \omega^2(y_{m-1} - 2 y_{m} + y_{m+1}) \\ y_m(0) = 0 \\ y'_m(0) = 0 \end{array} \right. \end{equation*}$

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

(2) $\begin{equation*} \left\{ \begin{array}{l} y_0 = 0 \\ y_{M+1} = 0 \\ y'_0 = 0 \\ y'_{M+1} = 0 \end{array} \right. \end{equation*}$

とします。

実際には， $y_m = y1_m, y'_m = y2_m$ とした, 2(M+2) 個の連立微分方程式に対し, オイラー法（ルンゲ・クッタ法）を利用します.

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

オイラー法でのアルゴリズムを下に示します。ルンゲ・クッタ法ではどうすればよいかは各自考えてみてください。

アルゴリズム

解を求める時刻の上限を $T$ とし、適当な自然数 $N$ を定めて $h=t/N$ とする。
$y1_m$ を配列変数 Y1[m] で， $y2_m$ を配列変数 Y2[m] で表す.
(配列 Y1[M+2], Y2[M+2], Y1_new[M+2], Y2_new[M+2]は実数型配列。)
T, h, t 及び関数 f1(t, x), f2(t, x1, x2, x3) は実数型
M, N, j, m は整数型(int)

T, M, N, R, ω (w 等で代用)を設定 (#define文で定数として定義）
実数型配列 Y1[M+2], Y2[M+2], Y1_new[M+2], Y2_new[M+2]を用意する
h = T/N（TおよびN双方が整数型でないか注意）

GLSC の準備（g_init(), g_device(), g_def_scale(), g_sel_scale(), ... 等）

（１. 初期値を与える）
m = 0, 1, 2,...., M+1 に対し(for 文)

    Y1[m] = 0
    Y2[m] = 0

を繰り返す
Y2[1]=0.25（課題１での初期値の設定）

t=0

m = 0, 1, 2,...., M+1 に対し初期値のグラフの描画 (for 文 + g_marker(), g_circle() 等)

（２. 計算）

(時間方向への繰り返し)
j = 0,1,2,.......,N-1 の順に (for 文）

    (全てのおもりについて計算)
    m = 1,2,.......,M の順に (for 文）

        Y1_new[m] = Y1[m] + h*f1(t, Y2[m])
        Y2_new[m] = Y2[m] + h*f2(t, Y1[m-1], Y1[m], Y1[m+1])

    を繰り返す

    m = 1,2,.......,M の順に (for 文）

        Y1[m] = Y1_new[m]
        Y2[m] = Y2_new[m]

    を繰り返す

    t =(j+1)*h　（時刻を進める）

    一定間隔毎に(例えば j%1000==0 の時に, 等)

        グラフを消して(g_cls())
        m = 0, 1, 2,...., M+1 に対し新たにグラフの描画 (for 文 + g_marker(), g_circle() 等)
        M+1 まで描き終わったら少々停止（g_sleep()）

を繰り返す

f1(t, x), f2(t, x1, x2, x3) は関数として定義

描画の時間間隔や g_sleep() の時間を適度に調整して, 滑らかな動きが見えるようにしてください.

課題１

M=1、ω=1、R=1 として、式(1)と式(2)をとくプログラムを作成せよ。ただし初期条件は、y₁=0.0、y₂=0.25 とする。
オイラー法の場合, 時間刻み幅は h=0.0001 (T=300等)として、1000ステップ毎に、おもりの位置 (m番目のおもりの位置は (m*R + y_m, 0))を g_marker で表示させること。

（ルンゲクッタ法の場合, 時間刻み幅は h=0.1 として、1ステップ毎に、おもりの位置 (m番目のおもりの位置は (m*R + y_m, 0))を g_marker で表示させること。）

次図ではグラフにする横軸の範囲は[ 0, 2 ]とした。

課題２

課題１のプログラムを M = 3, 10 にしてその様子を見よ。y2₁=0.7、それ以外のy1_m、y2_mはゼロを初期条件とする。（グラフの横軸の範囲は [ 0, (M+1)*R ] とする.）

課題３

課題２のプログラムのグラフ表示部分を、点 (m*R, y1_m) を直線で結ぶように書き換えよ。また, M をさらに増やした場合の様子を見よ. （グラフの縦軸の範囲は [ -0.5, 0.5 ] 等）

課題４

課題３のプログラムをルンゲ・クッタ法を用いて書き直せ。時間刻み幅は h=0.1 として、1ステップ毎に、おもりの位置を表示させること。

課題５

M = 10、y2₁=0.7、それ以外のy1_m、y2_mはゼロを初期条件として、
次の二つの図を上下に表示させ、アニメーションを行うプログラムを作成せよ。
上の図：課題２と同様の図
下の図：課題３と同様の図
ただし、数値計算の方法はルンゲ・クッタ法であるとし、時間刻み幅はh=0.1であるとする。

連結バネ系について。。。

あらゆる物質は分子, 原子（厳密には素粒子ですが）から出来ているのはご存知かと思います。
そのうち「固体」と呼ばれているもの（氷等の結晶や金属等）は, そのような分子（原子）間に働くの引力と斥力によって、あたかもバネで繋がれたかのように, それらが規則正しく（等間隔に）配置されている状態であります.（この「バネ」という比喩は実際の分子によく合致します.)
ですので今回の課題２,３は, このような「固体」を原子レベルで見たような事になっています.（厳密には違いますが, 大まかに…)
ω によってその「固体」の固さが変わる感じです.
課題２, ３では波のようなものが見えていたかと思いますが, この波が空気に伝搬して耳に入ると「音」として我々が認識することになります.