第6回(2009年6月3日)


今日の演習

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


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

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


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


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



          ・・・ (7)


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

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



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

          ・・・ (8)


ただし、
          ・・・ (9)


とした。

2階の常微分方程式は次のように連立の1階の常微分方程式に帰着できます。y1(t)=y(t), y2(t)=y'(t) とおくと、式(7)は、




          ・・・ (10)


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


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

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

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

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) は関数として定義するとよい.


課題1

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


課題2

ω=1として式(10) をオイラー法で解き、グラフを描け。グラフは横軸が 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 でもうまく解を近似できるような解法が求められます。そのような要求に応える解法としては、ルンゲ・クッタ法があります。



課題3

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


          ・・・ (11)


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


課題4

下図のような棒の長さが一定の振り子の、重力下での運動について、振り子の角度を 2πX とすると、X は

X'' = -sin(2πX)

 という運動方程式を満たす(摩擦は無視しています)。この方程式を適当な初期条件、X(0) = 0, X'(0) = a でシミュレーションし、GLSC を用いた振り子の動画を作成せよ.(デモをお見せします.)


レポート課題

先程単振動を考察しましたが, 実際のバネではバネの復元力以外におもりの速度 y' に比例する抵抗が働く。(おもりと床との間の摩擦力のようなものをイメージ)抵抗を考慮すると、式(8)ではなく次の式になる。(2γy' の項が新たに加わった。)



          ・・・ (12)


ここで γは抵抗力の強さを表す正の定数です。式(12)の問題を連立1階常微分方程式に変換し、オイラー法を用いて解き、GLSCをつかって「課題2のようなグラフ」か「課題4のような動画」を作成せよ. 各自 γ < ω での様子と γ > ω での様子を見比べること。作成したプログラムを rep0603.c というファイル名で report フォルダにおいてください. (プログラムの先頭にコメント文で名前と学籍番号を書く事. またプログラムの説明もコメント文として適宜記入する事.)

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


今回の課題は、オイラー法では十分に小さい h を用いないと計算結果がすぐにずれてしまいますが, 練習としてやってみてください. 次回, オイラー法より計算精度の高いルンゲ・クッタ法を紹介しますが, そのときに比較してみてください.

課題5

課題4の振り子が2つ相互作用する場合を考える. 例えば, 2つの振り子の角度を 2πX1, 2πX2 とすると、X1, X2, が

X1'' = -Asin(2πX1) - Csin(2π(X1-X2))

X2'' = -Bsin(2πX2) - Csin(2π(X2-X1))

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


今日学んだ事


戻る