九州大学理学部数学科2年
「計算数学概論」
12. 微分方程式

2009 年7月9日

1. ボールの運動方程式

重力が働く空間内で投げたボールの軌跡を微分方程式を解くことでシミュレーションする.
ここでは, 原点 (0, 0) から, 速度 (10, 10) で投げたボールの軌跡を求める.

1.1 アニメーションの準備

最初に空間内のボールを表示するためのアニメーション画面を準備する.
Animate関数の詳細は[ヘルプ] 等で確認せよ.
下の例では, ボールを赤丸で表示し, そのx座標に変数x1を使い, y座標に変数y1を使っている.
それぞれの値は, 0 から100までの間でスライダーで可変可能である.

In[3]:=

"20090709_1.gif"

Out[3]=

"20090709_2.gif"

1.2 ラグランジアン (L0)

"20090709_3.gif"

質点の質量 (m) と重力加速度 (g) の値を与えたときのラグランジアンを返す関数 L0[m, g] を定義する.

In[4]:=

"20090709_4.gif"

1.3 ラグランジュの運動方程式

"20090709_5.gif"

質量 10, 重力加速度 9.8 のラグランジアンを lex0に代入しておく.

In[5]:=

"20090709_6.gif"

Out[5]=

"20090709_7.gif"

"20090709_8.gif"

In[6]:=

"20090709_9.gif"

Out[6]=

"20090709_10.gif"

"20090709_11.gif"

In[7]:=

"20090709_12.gif"

Out[7]=

"20090709_13.gif"

"20090709_14.gif"

In[8]:=

"20090709_15.gif"

Out[8]=

"20090709_16.gif"

1番目の運動方程式

In[9]:=

"20090709_17.gif"

Out[9]=

"20090709_18.gif"

すなわち, x''[t] = 0 となる.

2番目の運動方程式

In[10]:=

"20090709_19.gif"

Out[10]=

"20090709_20.gif"

すなわち, y''[t] = -9.8 となる.

1.4 厳密解

問題 1. 微分方程式 x''[t] =  0, y''[t] = -9.8 の一般解を求めよ. また, 初期位置 (x[0], y[0]) = (0, 0), 初期速度 (x'[0], y'[0]) = (10, 10)
とするとき, 解は, y = - 0.049 "20090709_21.gif"+ x となることを確認せよ.

解のグラフは以下のようになる.

In[11]:=

"20090709_22.gif"

Out[11]=

"20090709_23.gif"

1.5 オイラー法による近似解法

"20090709_24.gif"

U[t] が与えられたとき, U[t + d] の値は, U[t+d] = U[t] + "20090709_25.gif" であるが, これを U[t+d] = U[t] + d * F[t, U[t]] で近似するのがオイラー法である.

ここでの F[t, u1, u2, u3, u4] をMathematicaの関数 EX1F[t, u1, u2, u3, u4] で定義する.

In[12]:=

"20090709_26.gif"

u1, u2, u3, u4 の値を代入する配列を準備する.

d を 0 .0001 とし,
u1[0] の値を ex1u1[[1]] に, u1[n*d] の値を ex1u1[[n + 1]] に入れる.
u2[0] の値を ex1u2[[1]] に, u2[n*d] の値を ex1u2[[n + 1]] に入れる.
u3[0] の値を ex1u3[[1]] に, u3[n*d] の値を ex1u3[[n + 1]] に入れる.
u4[0] の値を ex1u4[[1]] に, u4[n*d] の値を ex1u4[[n + 1]] に入れる.
配列の大きさ, は 20000 とする.

In[13]:=

"20090709_27.gif"

初期値を代入する.

ここでは, 初期値を(u1[0], u3[0])= (0.0, 0.0) 初期速度を (u2[0], u4[0])=(10.0, 10.0) とする.

In[16]:=

"20090709_28.gif"

オイラー法の計算

刻み幅 d , 現在値 t, u1[t], u2[t], u3[t], u4[t], 微分方程式の右辺の関数 F[t, u1, u2, u3, u4] が与えられたとき, オイラー法での U[t + d] の近似値  U[t] + d * F[t, U[t]] を計算する関数 EX1Euler[d, t, u1, u2, u3, u4, F] を定義する.

In[17]:=

"20090709_29.gif"

刻み幅をd = 0.0001 にして配列の値を計算する.

初期時間を t = 0 とし, 刻み幅を d = 0.0001 として, 準備した配列が一杯になるまで計算する.
配列の添字が n の値が, 時間 t = n*d のときの値である.
従って, 0.0001*20000 = 2 秒後までの軌跡が計算される.

In[18]:=

"20090709_30.gif"

値の座標 (x[t], y[t])=(u1[t], u3[t]) をグラフ表示する.

In[19]:=

"20090709_31.gif"

Out[19]=

"20090709_32.gif"

誤差(厳密解との差)を調べる.

In[20]:=

"20090709_33.gif"

Out[20]=

"20090709_34.gif"

値の座標をアニメーション表示する.

1.1 節で準備したアニメーションの中で, x1に ex1u1[[Floor[n]]] を,
y1に ex1u3[[Floor[n]]] を代入する文をAnimate関数内の最初に入れる.
これに伴い, スライダーで調整するのは, x1, y1 ではなく, n の値になる.
  最小のnは1で, 最大のnは20000である.
スライダーで指定された n の値は整数とは限らないので,
配列要素を参照するときは Floor[n] として小数点以下を切り捨てて参照する.

In[21]:=

"20090709_35.gif"

Out[21]=

"20090709_36.gif"

1.6 定理 (解の存在と一意性)

"20090709_37.gif"

(注. ) リプシッツ条件は, 微分方程式の解の存在や一意性だけでなく, 近似計算法において局所打ち切り誤差の評価にも利用される.

2. 単振り子の運動方程式

"20090709_38.gif"

2.1 アニメーションの準備

振り子の状態を描画するアニメーションを準備する.
ここでは, パラメータ t を鉛直方向, すなわち, y軸の負の方向からの偏角として,
0 度から360度までスライダーで動かせるようにしている.
また, 棒の長さを1としている. 偏角が s1 ラジアンのとき,
質点の座標は (x1, y1) = (Sin[s1], -Cos[s1]) となる.

In[22]:=

"20090709_39.gif"

Out[22]=

"20090709_40.gif"

2.2 ラグランジアン (L1)

In[23]:=

"20090709_41.gif"

Out[23]=

"20090709_42.gif"

"20090709_43.gif"

棒の長さ(l), 質点の質量 (m) と重力加速度 (g) の値を与えたときのラグランジアンを返す関数 L1[l, m,  g] を定義する.

In[24]:=

"20090709_44.gif"

2.3 ラグランジュの運動方程式

"20090709_45.gif"

弦の長さ10, 質量 10, 弦重力加速度 9.8 のラグランジアンを lex1に代入しておく.

In[25]:=

"20090709_46.gif"

Out[25]=

"20090709_47.gif"

ラグランジュの運動方程式の計算

In[26]:=

"20090709_48.gif"

Out[26]=

"20090709_49.gif"

すなわち, s''[t] = - 0.98 Sin[s[t]] となります.

2.4 厳密解

問題 2.   微分方程式 s''[t] = - 0.98 Sin[ s[t] ] について,
(1) s[t]が小さいと仮定し, Sin[s[t]]≈s[t] と仮定した場合の一般解を求めよ.
(2) s[t]が小さいとは限らない場合の一般解を求めよ. (3) 初期位置 s[0] = "20090709_50.gif", 初期速度 s'[0] = 0 とするとき, 周期を求めよ.

2.5 オイラー法による近似解法 (1)

"20090709_51.gif"

ここでの F[t, u1, u2] をMathematicaの関数 EX2F[t, u1, u2] で定義する.

In[27]:=

"20090709_52.gif"

u1, u2, u3, u4 の値を代入する配列を準備する.

d を 0 .01 とし,
u1[0] の値を ex2u1[[1]] に, u1[n*d] の値を ex2u1[[n + 1]] に入れる.
u2[0] の値を ex2u2[[1]] に, u2[n*d] の値を ex2u2[[n + 1]] に入れる.
配列の大きさ, は 20000 とする.

In[28]:=

"20090709_53.gif"

初期値を代入する.

ここでは, 初期位置を u1[0] = "20090709_54.gif" 初期速度を u2[0] = 0 とする.

In[30]:=

"20090709_55.gif"

オイラー法の計算

刻み幅 d , 現在値 t, u1[t], u2[t], 微分方程式の右辺の関数 F[t, u1, u2] が
与えられたとき, オイラー法での U[t + d] の近似値  U[t] + d * F[t, U[t]] を
計算する関数Ex2Euler[d, t, u1, u2,  F] を定義する.

In[31]:=

"20090709_56.gif"

刻み幅をd = 0.003 にして配列の値を計算する.

初期時間を t = 0 とし, 刻み幅を d = 0.003 として,
準備した配列が一杯になるまで計算する.
配列の添字が n の値が, 時間 t = n*d のときの値である.
従って, 0.003*20000 = 60 秒後までの軌跡が計算される.

In[32]:=

"20090709_57.gif"

角度θ = s[t] = u1[t]  (d=0.003, t<60) をグラフ表示する.

誤差が積もって途中(n=10000 (30秒後) くらい)から解が発散していることがわかる.

In[33]:=

"20090709_58.gif"

Out[33]=

"20090709_59.gif"

2.6 オイラー法による近似解法 (2)

刻み幅をd = 0.002 にして配列の値を計算する.

初期時間を t = 0 とし, 刻み幅を d = 0.002 として,
準備した配列が一杯になるまで計算する.
配列の添字が n の値が, 時間 t = n*d のときの値である.
従って, 0.002*20000 = 40 秒後までの軌跡が計算される.

In[34]:=

"20090709_60.gif"

角度θ = s[t] = u1[t]  (d=0.002, t<40) をグラフ表示する.

最後(40秒後)でも発散はしていないが, 触れ幅がだんだん大きくなっていることがわかる.

In[35]:=

"20090709_61.gif"

Out[35]=

"20090709_62.gif"

アニメーション表示する.

2.1 節で準備したアニメーションの中で, s1に ex2u1[[Floor[n]]] を代入する文を
Animate関数内の最初に入れる. これに伴い, スライダーで調整するのは, x1, y1 ではなく,
n の値になる. 最小のnは1で, 最大のnは20000である.
スライダーで指定された n の値は整数とは限らないので,
配列要素を参照するときは Floor[n] として小数点以下を切り捨てて参照する.
再生ボタンで再生すると早く動きすぎるので, 右から2番目の[遅く] ボタンで再生速度を遅くする.

In[36]:=

"20090709_63.gif"

Out[36]=

"20090709_64.gif"

2.7 二次のルンゲクッタ法による近似解法

微分方程式 "20090709_65.gif" に対して,
U[t] が与えられたとき, U[t + d] の値は, U[t+d] = U[t] + "20090709_66.gif" であるが, これを U[t+d] = U[t] + d * F[t, U[t]] で近似するのがオイラー法であった.
ここでは, k1 = d * F[t, U[t]], k2 = d*F[t+d, U[t]+d*F[t,U[t]]] = d * F[t+d, U[t] + k1] として,
U[t+d] = U[t] + "20090709_67.gif" (= "20090709_68.gif")
で近似する方法を考える. この方法を二次のルンゲクッタ法と呼ぶ.

二次のルンゲクッタ法の計算

刻み幅 d , 現在値 t, u1[t], u2[t], 微分方程式の右辺の関数 F[t, u1, u2] が与えられたとき,
二次のルンゲクッタ法での U[t + d] の近似値を計算する
関数 EX2RungeKutta2[d, t, u1, u2,  F] を定義する.

In[37]:=

"20090709_69.gif"

刻み幅をd = 0.003 にして配列の値を計算する.

初期時間を t = 0 とし, 刻み幅を d = 0.003 として,
準備した配列が一杯になるまで計算する.
配列の添字が n の値が, 時間 t = n*d のときの値である.
従って, 0.003*20000 = 60 秒後までの軌跡が計算される.

In[38]:=

"20090709_70.gif"

角度θ = s[t] = u1[t]  (d=0.003, t<60) をグラフ表示する.

2.5 オイラー法による近似解法 (2) の刻み幅 d=0.002 よりも大きいにも関わらず最後(60秒後)まで安定している.

In[39]:=

"20090709_71.gif"

Out[39]=

"20090709_72.gif"

刻み幅をd = 0.0037 にして配列の値を計算する.

初期時間を t = 0 とし, 刻み幅を d = 0.0037 として,
準備した配列が一杯になるまで計算する.
配列の添字が n の値が, 時間 t = n*d のときの値である.
従って, 0.0037*20000 = 74 秒後までの軌跡が計算される.

In[40]:=

"20090709_73.gif"

角度θ = s[t] = u1[t]  (d=0.0037, t<74) をグラフ表示する.

最後の2000ステップだけを表示している.
誤差が積もって途中(n=19200 (70秒後) くらい)から解が発散していることがわかる.

In[41]:=

"20090709_74.gif"

Out[41]=

"20090709_75.gif"

2.8 二次のルンゲクッタ法の精度について

"20090709_76.gif"

U[t + d] のテーラー展開

"20090709_77.gif"

F[t + d1, u1[t] +d2, u2[t]+d3] のテーラー展開

"20090709_78.gif"

二次のルンゲクッタ法の近似式

"20090709_79.gif"

2.8 (発展課題) 四次のルンゲクッタ法による近似解法

次の関数 Ex2RungeKutta4[d, t, u1, u2, F] を用いて近似計算を行う方法を
四次のルンゲクッタ法を呼ぶ.

In[42]:=

"20090709_80.gif"

(注.) 四次のルンゲクッタ法の局所打ち切り誤差は "20090709_81.gif"以下である.

刻み幅をd = 0.1 にして配列の値を計算する.

初期時間を t = 0 とし, 刻み幅を d = 0.1 として,
準備した配列が一杯になるまで計算する.
配列の添字が n の値が, 時間 t = n*d のときの値である.
従って, 0.1*20000 = 2000 秒後までの軌跡が計算される.

In[43]:=

"20090709_82.gif"

角度θ = s[t] = u1[t]  (d=0.1, t<2000) をグラフ表示する.

最後の200ステップだけを表示している.
大きな刻み幅にも関わらず最後(2000秒後)まで安定している.

In[44]:=

"20090709_83.gif"

Out[44]=

"20090709_84.gif"

3. バネの運動方程式

摩擦のない平面上に横に置いたバネの片方の端に重り (質点) をつなぎ
もう片方の端を壁に固定する. 重りを横に動かして放すと左右の振動を始める.
この運動を微分方程式を解くことでシミュレーションする.
ここでは, 長さ10 cm のバネの固定端を原点にし, 重りの初期位置を15する.
また, バネ定数の値は 10 N/cm とする.

3.1 アニメーションの準備

バネの延び縮みに伴う質点の運動を表示するためのアニメーション画面を準備する.
下の例では, ボールを赤丸で表示し, そのx座標に変数x1を使っている.
x1の値は, 0 から30までの間でスライダーで可変可能である.
初期位置を x1 = 10 としている.

In[45]:=

"20090709_85.gif"

Out[45]=

"20090709_86.gif"

3.2 ラグランジアン (L2)

In[46]:=

"20090709_87.gif"

Out[46]=

"20090709_88.gif"

"20090709_89.gif"

バネの長さ(l), 質点の質量 (m) とバネ定数(k) の値を与えたときのラグランジアンを返す関数 L2[l, m,  k] を定義する.

In[47]:=

"20090709_90.gif"

3.3 演習問題

問題 3. 3.2のラグランジアンから, ラグランジュの運動方程式(x[t] の満たす微分方程式)を求めよ.

問題4. 初期位置 x[0]=15, 初期速度 x'[0]=0 として, 200秒のシミュレーションを行いたい.
(1) オイラー法により近似解を求めよ. (2) 二次のルンゲクッタ法により近似解を求めよ.

問題5. 問題4. で得られた近似解を 3.1のアニメーションを使って表示せよ.

4. 課題提出方法

練習問題1,2,3,4,5についてMathematicaのノートブックで解答を作成する.
完成したノートブックファイルをWebCTの課題提出箱に添付ファイルとして提出する.
答案欄には感想などをひとこと入力する.

5. 発展課題

発展問題1. 二重振り子のシミュレーションを行え

アニメーションの準備

In[48]:=

"20090709_91.gif"

Out[48]=

"20090709_92.gif"

ラグランジアン (L3)

In[49]:=

"20090709_93.gif"

発展問題2. 連接バネのシミュレーションを行え

アニメーションの準備

In[50]:=

"20090709_94.gif"

Out[50]=

"20090709_95.gif"

ラグランジアン (L4)

In[51]:=

"20090709_96.gif"

6. (おまけ) シミュレーション結果を動画ファイルに保存する.

Mathematica には画像のリストを動画ファイルとして保存する関数Exportがある.
ここでは, 2.6 節の単振り子のシミュレーションの結果を動画 (AVIファイル) に
保存する例を紹介する.
2.6 節の Animate関数を実行している部分において, DynamicModule を Module へ,
Animate を Table へ変更すると出力画像のリストが得られる.
動画は1秒間に15枚の画像を使うので, 10 秒間の動画作成には150枚で十分である. 従ってパラメータの範囲を {n, 1, 20000} でなく, {n, 1, 150} と小さくする.
出来た画像のリストを Export 関数に渡す. ファイル名として, "Documents/sample.avi"
を指定すると, [書類] フォルダ中に "sample.avi" が保存される.

(注意) 学生は100MBのファイル容量制限があるので, 大きな動画は保存出来ない.
動画ファイルの作成にも時間がかかるので, パラメータの範囲には注意が必要である.
( 10 秒の動画の場合, AVI形式の動画ファイルは作成に15秒程度かかり大きさが約60MBになる.
  フラッシュ動画 (FLV形式) の作成には1分程度かかり, 大きさは約300KBになる.)

In[52]:=

"20090709_97.gif"

Spikey Created with Wolfram Mathematica 7.0