計算数理A演習第12回

今日の演習

  • 連結バネ系

連結バネ系

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

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

(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(), ... 等)

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

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

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

t=0

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

(2. 計算)

(時間方向への繰り返し)
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() の時間を適度に調整して, 滑らかな動きが見えるようにしてください.


課題1

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

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

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

課題2

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


課題3

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


課題4

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


課題5

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


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

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