力学演習 II

シミュレーション (1) 参考資料
福井・大信田 (応用数理工学科)

1

2006-10-18

共同レポートについて

2 人で共同レポートを出してよいのは、コンピュータを使わない場合に限ります。コンピュータを使う場合は 1 人で 1 通のレポートを書いてください。言い換えれば「人間 (指令&記録係)+人間 (計算係)」
または「人間+コンピュータ」というチームは OK ですが、「人間+人間+コンピュータ」は反則です。

2無次元化

力学などの物理的な問題の数値計算では、無次元化をおこなう場合が多い。その理由を説明し、一般的な手順を示そう。
いくつかの例外を除き、すべての物理量は次元をもっている。次元が異なる量どうしを足したり引いたり等号で結んだりすることはできない。特に、物理量 (次元のある量) と普通の数を等号で結ぶことはできない。単位というものが必要になるのは、そのギャップを埋めるためであって、

L = 30 cm,

S = 2πR2(#1a)

S = 2π(#1b)

は OK だが

L = 30,はダメだ (30 も 2π も普通の数であることに注意)。一般的に言えば、物理量と数は[物理量] = [数] × [単位](#2)という形で結びつけられ、式 (#1a) は確かにこの形になっている。式 (#1b) のように「数」と「物理量」を直接結びつけるのはルール違反であって、次元のない世界 (数) と次元のある世界 (物理量) を正しく結びつけるには、単位というものを仲介することが必要なのである。
さて、数値計算は読んで字のごとく「数」の計算なので、次元がない世界を扱う。したがって、物理的な問題の数値計算をおこなうには、物理量 (有次元) を数 (無次元) に変換しなければならない。これを無次元化という。


ポイントその 1
無次元化とは、 物理量 (次元あり) → ふつうの数 (次元なし) という変換のことだ。



最も安易な無次元化の方法は、単位を黙って消してしまうことだが、これはいくらなんでも乱暴すぎる (15cm と 15mm の区別がつかなくなる)。まともな方法は、式 (#2) を変形して[数 (無次元量)] =[物理量][単位](#3)とすること — たとえば位置座標 x という物理量に対して、無次元化された位置座標をx̃ = x/cm(#4)

力学演習 II (2006-10-18) シミュレーション (1) 参考資料: 無次元化2のようにあらわすことだ。この考え方はさらに一般化できる。式 (#3) の右辺の分母は、必ずしも cm とか mm とかいった普通の単位である必要はなく、適切な次元をもつ定数であれば何でもよいのだ。つまり[数 (無次元量)] =[物理量][基準となる定数](#5)というようにして、物理量を無次元量に変換できる。基準となる定数をうまく選ぶことで、方程式を簡単化したり、計算すべき問題の量を減らしたりすることができるため、無次元化といえば (#5) のような変換を意味するのが普通である。
ポイントその 2
式 (#5) を用いて、物理量 (有次元の量) を無次元量に変換する。

無次元化の手順
マニュアル的にまとめてみよう:(1) 無次元化を必要とする変数はどれなのか明らかにする。
(2) 問題の物理的な意味を考え、基準量の候補を検討する。
(3) 式 (#5) を用いて運動方程式を変換する。
(4) 上記 (2)(3) の結果をふまえて基準量を設定する。
(5) 方程式全体を無次元変数と無次元パラメータで書かれた形に整理する。
これだけでは何だか分からないので、例をまじえてもう少し詳しく説明しよう。
ステップ 1: 無次元化を必要とする変数はどれなのか明らかにする。
たとえば人工衛星の運動では、位置ベクトルも時間も有次元の変数なので、無次元化が必要である。
他方、振り子の運動方程式は、角度 θ についての微分方程式になる。角度は無次元だから、θ については無次元化は不要だ。時間だけ無次元化すればよい。
ステップ 2: 問題の物理的な意味を考え、基準量の候補を検討する。
このステップは、問題に対する見通しがどれだけあるかに依存する。全く見当がつかない場合はとばして次のステップに進んでもよい (が、そのぶん、あとで苦労することになるかもしれない)。
他方、問題を大胆に簡単化することにより手計算で解けるようになるなら、その結果を利用することが考えられる。
いくつか具体例を見てみよう。

• 月に向かうロケットの運動の場合、長さの基準量の候補としては、地球の半径か、または月の軌道半径が考えられる。時間の基準量としては、これらの軌道半径に対応する公転周期 (またはそれから係数を除いたもの) をとるのが自然で、これは地球の引力による円運動を仮定すれば手計算で求められる。

力学演習 II (2006-10-18) シミュレーション (1) 参考資料: 無次元化

3

• バネでつるされた質点の 2 次元運動は、振り子運動とバネの伸縮が複合した複雑な運動になpるが、単純な振り子におきかえて考えると、周期は 2π `/g で、これ (または係数を除いたp`/g) が時間の基準量のひとつの候補になる。他方、1 次元のバネ質点系におきかえて考えpると、周期は 2π m/k で、これも時間の基準量の候補になる。
長さの基準量としては、バネの自然長をとるのがよいだろう。

なお、これらの例にあるように、基準量には複数の候補がある場合が珍しくない。このことについてはステップ 4 およびステップ 5 で再び触れる。
ステップ 3: 式 (#5) を用いて運動方程式を変換する。
たとえばバネ振り子の方程式を例に説明しよう:pd2 xx2 + y2 − ` p=−kx2dtx2 + y 2
pyd2 y− mgx2 + y 2 − ` pm 2 = −k2dtx + y2m(#6a)(#6b)

今の場合、長さの基準量は ` 以外に候補がない。時間に関しては基準量の候補が複数あって一つに絞れないため、とりあえず基準量を τ とおく。式 (#5) にしたがい、(x̃, ỹ) =(x, y),`t̃ =tτすなわち(x, y) = (x̃`, ỹ`),ddt̃ dd== τ −1dt
dt dt̃
dt̃(#7)のように変数変換して、運動方程式に代入する。結果を整理すると

px̃
m` d2 x̃
2 + ỹ 2 − 1 p=−k`x̃
222
τ dt̃
x̃ + ỹ 2pỹm` d2 ỹ2 + ỹ 2 − 1 p=−k`x̃
− mg22τ dt̃
x̃2 + ỹ 2(#8a)(#8b)となり、さらに両辺を左辺の係数 m`/τ 2 で割ると

d2 x̃となる。

px̃
2 k2 + ỹ 2 − 1 p=−τx̃
2mx̃2 + ỹ 2dt̃

pỹ
d2 ỹ
g2 kx̃2 + ỹ 2 − 1 p− τ22 = −τ m22`
x̃ + ỹdt̃

(#9a)
(#9b)

ステップ 4: ステップ 2 および 3 の結果をふまえて基準量を設定する。
上記の例の場合、τ を決める必要があるが、どう決めるかには選択の余地がある。そこで、ステップ 3 の最後の式 (#9) の右辺を見て、どこかの係数が 1 になるようにする。たとえば τ 2 k/m = 1となるようにするにはp(#10)
τ = m/kとすればよい。これはステップ 2 でのふたつの候補のうちの片方にちょうど一致している (ただし2π をつけないほうが式がきれいになる)。

力学演習 II (2006-10-18) シミュレーション (1) 参考資料: 無次元化
もちろん、式 (#9) の右辺の最後の項に着目して
pτ = `/g4(#11)とすることも可能だ。どちらを選ぶかは数値計算をする本人が決断すべきことであり、何にピントをあわせて解析したいかによる。
ステップ 2 を飛ばした人は、ここでいったん振り返って、基準量の物理的意味を考えてみること。
今の例で言えば、式 (#10) や式 (#11) の意味を考えてみろ、ということだ。式 (#10) は 1 次元のバネ質点系の周期のオーダーの量で、式 (#11) は振り子の周期のオーダーの量であることが分かれば OK だ。
ステップ 5: 方程式全体を無次元変数と無次元パラメータで書かれた形に整理する。
たとえば (#10) のほうを選択したとしよう。これにより方程式 (#9) を書き直し、無次元の量 (無次元化された変数および無次元パラメータ) だけであらわす。結果は

d2 x̃となり、ここでpx̃
2 + ỹ 2 − 1 px̃
=−22x̃ + ỹ 2dt̃p
ỹd2 ỹ2 + ỹ 2 − 1 p−G=−x̃
2x̃2 + ỹ 2dt̃(#12a)(#12b)

mg(#12c)k`と置いた。この G は、問題ごとに異なる値をとる可能性のある定数 (m, ` など) を含む無次元の定数であり、このようなものを無次元パラメータという。
G=

ところで、ここで考えたバネ振り子の例では、時間の基準量の候補 (時間の代表スケール) が 2 つあったことを思い出そう。式 (#10) の右辺を τS とおき、式 (#11) の右辺を τG とおくと、p 2m/k √τSτS
p(#13)
== G
つまり
G=
τG
τ
`/gGである。 このように、無次元パラメータはふたつの代表スケールの比であらわすことができる。このことと、τS および τG の物理的意味をあわせて考えると、たとえば G が大きい場合と小さい場合の違いをある程度まで予想することが可能になる。

• もし G  1 なら、τS  τG ということだから、振り子が何度も往復するあいだにバネがゆっくり伸縮するような運動が予想できる。
• もし G  1 なら、τS  τG だから、振り子が 1 往復するあいだにバネは何度も伸縮を繰り返す。
• もし G が「オーダー 1」なら、バネの伸縮と振り子運動がからみあった複雑な運動になるだろう。



ポイントその 3

無次元化するだけならマニュアルどおりやればできる。でも、本当に大事なのは、基準量や無次元パラメータの意味を考えることだ。



力学演習 II (2006-10-18) シミュレーション (1) 参考資料: 常微分方程式の数値解法

5

常微分方程式の数値解法

3

例題として、スカイダイバーの運動方程式∗ を数値的に解く方法を説明しよう。
パラシュートをつけたスカイダイバーの落下速度 v は、

mdv
= mg − βv 2dt(#14)という式に従う (ただし v > 0)。ここで β は定数で、空気の密度やパラシュートの大きさなどで決まる。
まずは式 (#14) を無次元化しよう。

(1) 無次元化を必要とする変数は v と t である。
(2) 速度の基準量の候補としては、終端速度 V =pmg/β が考えられる。

時間の基準量としては、等加速度運動 (自由落下) から等速運動に変化するのにかかる時間をとることが考えられるが、すぐにはその値は分からない。

(3) 時間の基準量をとりあえず τ と置き、速度の基準量として終端速度 V をとることにして、
tvt̃ =ṽ = ,Vτのように変数変換し、式 (#14) を書き直す。
(4) 上記 (3) の結果を見て、係数が 1 になるように τ を定める。
(5) 結果を整理すると、ṽ についてのdṽ
= 1 − ṽ 2dt̃(♦)という常微分方程式になる。式 (♦) には無次元パラメータは含まれない。
次に式 (♦) の数値解法を考える。以下では、面倒なので ˜ を省略し、式 (♦) の右辺を f と置いてdv= f,dt

f = 1 − v2という形で考える。もちろん v = v(t) である。
式 (♦) の解を数値的に求めるには、高校数学でおなじみの 「将棋倒し」の考え方 を用いる† 。つまり、ある時刻における v の値が分かっているとして、これから ∆t だけ経過したあとの v の値を求める方法を示せばよい。時間を ∆t ごとに区切り、j 番目の時刻を tj としよう。初期条件により v(t0 ) は分かっているから、v(tj ) から v(tj+1 ) を求める方法さえあれば、すべての j に対して v(tj ) を求めることができる。


ポイントその 4

一般的に v(tj ) 7→ v(tj+1 ) の計算方法が分かれば、v(t0 ) 7→ v(t1 ) 7→ v(t2 ) 7→ · · · のような「将棋倒し方式」により、すべての値が計算できる。

∗†

今年度の常微分方程式演習を履修した人は、4/25 のプリントも参考にするとよい。
数列の漸化式とか数学的帰納法とかの考え方を思い出そう。

力学演習 II (2006-10-18) シミュレーション (1) 参考資料: 常微分方程式の数値解法

6

このような「将棋倒し」の方法を得るために、式 (♦) をv(tj+1 ) − v(tj )=f∆t(#15)のように書き直し、この式を v(tj+1 ) について解く。結果はv(tj+1 ) = v(tj ) + f ∆t(#16)となり、ここで f = f (v(tj )) と考えれば、右辺は v(tj ) から計算できるので、v(tj ) 7→ v(tj+1 ) が得られたことになる。
式 (#16) を用いた計算の手順は、たとえば次のようになるだろう:∆t を適当に決める (たとえば ∆t = 0.05)
t, v, f という記入欄をもつ表を用意する
表の何段目を見ているかということをあらわす変数 j を用意
とりあえず j = 0 と設定 (いちばん最初の段を 0 段目とする)
表の最初の段に t および v の初期値を記入※反復開始
現在の v の値から f = 1 − v 2 を計算
vnew = v + f ∆t を計算
tnew = t + ∆t を計算
vnew および tnew を新しい段に書き込む
j をひとつ増やす (新しい段に移動)
記入欄がなくなったら終了※に戻る

この手順書では f の値を f (v(tj )) として計算している。この方法は確かに簡単だが、あまり精度がよくないのが欠点である。
精度を高めるには、次のような 2 次の Taylor 展開を利用する:

1 d2 v 2dv∆t +∆t + · · ·dt
2 dt2
1 df 2= v(tj ) + f ∆t +∆t + · · ·2 dt 1 df∆t ∆t≒ v(tj ) + f +2 dtv(tj+1 ) = v (tj + ∆t) = v(tj ) +(Taylor 展開)(f = dv/dt)(#17)

ここで f = f |t=tj = f (v(tj )) であり、また tj + 21 ∆t = tj+1/2 とすると


1 df
1 df
f v(tj+1/2 ) = f |t=tj +∆t/2 = f +∆t + · · · ≒ f +
∆t
2 dt
2 dt

(#18)

だから、これを式 (#17) に代入して

v(tj+1 ) = v(tj ) + f v(tj+1/2 ) ∆t(#19)

力学演習 II (2006-10-18) シミュレーション (1) 参考資料: 運動方程式の 1 階化7という式を得る。これは中点公式などと呼ばれる方法で、式 (#16) の f を f (v(tj+1/2 )) でおきかえた形になっている。式 (#19) は 2 次の Taylor 展開に基づいているので、もとの式 (#16) よりも精度が高くなっているはずである。
中点公式を用いるには、手順書の反復部分を次のように書き直す:※反復開始
現在の v の値から f = 1 − v 2 を計算
vhalf = v + f ∆t/2 を計算
2fhalf = 1 − vhalfを計算
vnew = v + fhalf ∆t を計算
tnew = t + ∆t を計算
...
※に戻る

2 次の Taylor 展開を利用する方法としては、ここで示した中点公式のほかに、台形公式というものも存在する。さらに、3 次や 4 次の Taylor 展開を利用する方法も知られている。詳しくは各自で調べてみること。

運動方程式の 1 階化

4解法 (#16) や (#19) などを用いれば 1 階の常微分方程式の数値解を求めることができるが、この方法は、2 階以上の常微分方程式には適用できない。ところが運動方程式は 2 階の常微分方程式なので、このままでは数値解が求められないことになる。
そこで、変数をうまく置き直して、運動方程式を 1 階の常微分方程式系に直す。



ポイントその 5運動方程式 (2 階) → 1 階の常微分方程式に変換すれば、式 (#16) などが利用できる。



たとえば 1 次元のバネ質点系の運動方程式は

md2 x= −kxdt2

(#20)と書けるが、これを 1 階の連立常微分方程式系に直すために

p=mdxdt(#21)と置く‡ 。これにより、運動方程式は" # "#d xp/m=dt p−kx(#22a)

これは質点の運動量である。位置座標に対しては運動量、角度に対しては角運動量というように、正準変数を用いて 1 階化するのがよい。
‡

力学演習 II (2006-10-18) シミュレーション (1) 参考資料: 運動方程式の 1 階化という形に書き直せる。無次元化した形で書くと" # " #d x̃p̃=dt̃ p̃−x̃

8(#22b)のようになる。
式 (#22) は 連立常微分方程式系 である。このような場合の数値計算では、(x, p) を常にひ̇と̇ま̇と̇ま̇り̇のものとして扱うようにする。例を以下に示す:※反復開始
現在の (x, p) の値から (v, f ) = (p, −x) を求める
(xhalf , phalf ) = (x, p) + (v, f )∆t/2 を計算(vhalf , fhalf ) = (phalf , −xhalf ) を求める
(xnew , pnew ) = (x, p) + (vhalf , fhalf )∆t を計算
tnew = t + ∆t を計算
...
※に戻る

初心者の場合、つい「x だけ先に計算しよう」などと考えたくなるようだが、これでは順番が混乱してしまう。この方程式の従属変数は、(x, p) という 1 つのもの (位相空間上の 1 点をあらわす 1 つのベクトル) なのだから、むりに分けようとせず、いつでも並列に扱うようにすれば、混乱を防止できるはずだ。
ポイントその 6
連立常微分方程式系の数値計算では、従属変数をひ̇と̇ま̇と̇ま̇り̇のものとして扱う。

