第9回(2009年6月24日)


今日の演習



連結バネ系


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


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



・・・ (1)

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



          ・・・ (2)

とします。

実際には, ym = y1m, y'm = y2m とした, 2(M+2) 個の連立微分方程式に対し, オイラー法(ルンゲ・クッタ法)を利用します.

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

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


アルゴリズム

解を求める時刻の上限を T とし、適当な自然数 N を定めて h = T/N とする。
y1mを配列変数 Y1[m] で, y2mを配列変数 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.25、それ以外のy1m、y2mはゼロを初期条件とする。(グラフの横軸の範囲は [ 0, (M+1)*R ] とする.)

課題3

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

課題4

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

レポート課題

課題2〜4のうちの一つの課題についてのプログラムを rep0624.c というファイル名で report フォルダにおいてください. (プログラムの先頭にコメント文で名前と学籍番号と課題の番号を書く事. またプログラムの説明もコメント文として適宜記入する事.)

自動回収日時:2009年7月7日(火)午後4時



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

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

前回の課題4の補足。。。

前回の課題4では以下のような, バネは大きく伸び縮みすると急激に硬くなるような場合の強制振動系を扱いました.

dy/dt = v

dv/dt = - 0.3*v - y3 + B*cos(t)

この方程式は「ダフィン方程式」と呼ばれる, カオスを見るのに適した典型的な微分方程式(を少し変形したもの)です. このような方程式は、一般的に「解が解析的な関数では記述できない(非可積分)」という特徴を持っています. むしろ世の中の事を微分方程式で記述(モデル化)した場合(物理や工学で行われる事です. ), そこに出てくる微分方程式の多くは「非可積分」です. (微分方程式はある種自然界の「ルール」を記述しますが, そのルールが分かっただけでは実際に何が起こるのかまでは分からない, という事です.)しかし厳密な解が得られくても、その方程式の数理構造は理解できて、その理解に数値計算は重要な役割を持ってます. (これは数値計算の有用性の一部です.)


戻る