7.動く振り子

拘束条件が時間とともに変化する場合、拘束力 fc は式(5)となり、初期値 x0,˙x0 に対する拘束条件は式(1)となる。
この章のシミュレーション(クリックで計算開始)

前章の議論を、時間変化する拘束条件の場合に拡張したい

動く振り子の運動 x(t) を計算したい。これは、前章で扱った「振り子」において、中心点 a が時間とともに動くようになったものである。拘束条件 GG(t,x)|xa(t)|l=0 となり、時刻 t に依存することになる(第4章で扱った動く壁との衝突の場合と同様)
おもりの運動 x(t) を計算するには、前章同様、「初期値 x0,˙x0 に対する拘束条件」と「拘束力 fc 」の2つが分かればよい。そこでこの章では、 x(t) を実際に計算するために、以下の2つの節に分けて議論を行う: 7.1:7.2:x(t) 加えて、7.3節で、計算方法の補足を行う。

7.1 初期値に対する拘束条件と拘束力

この節では、初期値 x0,˙x0 に対する拘束条件(1)、および、拘束力の公式(5)を導く。議論の流れは、前章と全く同じである。

初期値 x0,˙x0 に対する拘束条件:式(1)

初期値 x0,˙x0 に対する拘束条件は、前章と同様、 t=0 において G=˙G=0 を満たすことである。ただし、 ˙G は、おもりの位置における時間微分なので ˙GddtG(t,x(t)) となる。よって、以下の条件を得る:(第2式は以下の【7.1-注1】の式(2)そのもの) G=0˙GtG+(G)T˙x=0}(1) ただし、 tG=GtG=(Gx)T である。

【7.1-注1】合成関数の微分公式

任意の関数 G(t,x) を考え、 xt に依存するとする。この時 G(t,x(t)) は、 t に関する1変数関数とみなせる。これを t で微分すると、以下のようになる: ddtG(t,x(t))=Gt+Gxdxdt(2)

導出

時間を δt だけ変化させたときの G の変化量 δG は、以下のように書ける:(第4章の【4.1-注4】参照)
δGGtδt+GxδxGtδt+Gxdxdtδt=(Gt+Gxdxdt)δt これを以下の dGdt の定義と見比べればよい: δGdGdtδt

拘束力 fc :式(5)

拘束力 fc の導出も、前章の議論をなぞればよい。まず、初期値 x0,˙x0 が拘束条件(1)を満たしている時、後はそれ以降の時刻で、拘束条件の2階微分 ¨G=0(3) が成り立っていれば、拘束条件は常に成立する。この式は加速度 ¨x を含むので、拘束力 fc に対する条件を与えている。しかし、これ1つだけでは fc を完全に決めることはできないので、物理的な条件を加える必要があり、それは、ダランベールの原理なのであった:λ は未定乗数) fc=λG(4) この式と式(3)を連立して λ を求めれば、 fc が確定する。実際に解くと、以下の【7.1-注2】の式(5)のようになる。

【7.1-注2】時間に依存する場合の拘束力:式(5)

時間に依存する拘束条件 G(t,x)=0 が課されている時、拘束力 fc は以下のようになる: fc=λGλm[ddttG+(ddtG)T˙x]+(G)Tf|G|2(5) ただし、 f は重力などの外力である。

導出

拘束条件(3)の微分を実際に行うと、以下のようになる:赤字部分は式(1)の第2式) 0=ddt˙G=ddt(tG+(G)T˙x)=ddttG+(ddtG)T˙x+(G)T¨x
この式に、運動方程式と式(4)を合わせたもの: m¨x=f+λG を代入して ¨x を消去すれば、 λ が決まり、式(5)のものに一致する。

7.2 物体の運動 x(t) の計算

物体の運動 x(t) を計算するには、初期値 x0,˙x0 を拘束条件(1) G=˙G=0 を満たすのようにとったうえで、運動方程式 m¨x=f+fc を解けばよい。ただし、拘束力 fc は式(5)で与えられる。
この節では例題として、冒頭で述べた動く振り子と、うねる床を考える。

動く振り子:運動方程式(7)、初期条件(6)

冒頭で述べた通り、動く振り子の拘束条件は G|xa(t)|l=0 となり、 G の微分は以下のようになる:ta のみに作用し、 d/dtx,a 両方に作用することに注意) tG=ˆxaT˙a=(xa)T˙alddttG=(˙x˙a)T˙a+(xa)T¨alG=ˆxa=xalddtG=˙x˙al
初期値 x0,˙x0 に対する拘束条件は、上式を式(1)に代入すればよい:a0=a(0)˙a0=˙a(0) |x0a0|=l(x0a0)T(˙x0˙a0)=0}(6) 拘束力 fc は、式(5)に G の微分を代入すればよい:
fc=m|˙x˙a|2+(xa)T(fm¨a)lˆxa 運動方程式 m¨x=f+fc に、この式と重力 f=mg を代入すれば、最終的に以下が得られる: m¨x=mgm|˙x˙a|2+(xa)T(g¨a)lˆxa(7)
振り子の中心を a(t)=[0asinbt]T のようにとりa,b は定数)、数値計算を行うと右図のようになる。グリッドの間隔は 1m である。

波打つ床

もう1つの例として、2次元平面上のグラフ y=f(t,x) で表される曲線上に、おもりを拘束する場合を考える。拘束条件 GGf(t,x)y=0 のようにとることができ、その微分は以下のようになる: tG=tf,G=[xf1]ddttG=(2t+˙xxt)f,ddtG=[txf+˙x2xf0]
初期値 x0,˙x0 に対する拘束条件は式(1)から決まり、拘束力 fc は式(5)から決まる。しかし、きれいな形になるわけでもないので結果については割愛する。
関数 f として例えば f(t,x)=asinbtsincx+dx2 としてa,b,c,d は定数)、数値計算を行うと右図のようになる。グリッドの間隔は 1m であり、下向きに重力がかかっている。

7.3 (計算手法の補足)

拘束力は、式(5)で与えられるわけだが、拘束条件 G の2階微分を含むので、手計算で求めるのは大変である(前章で扱った時間に依存しない場合も同様である)。そこでこの節では、微分を数値的に計算する方法と、差分化することで運動方程式から2階微分を消す方法について述べる。

微分を数値的に計算する方法

数値計算を行う際、拘束条件 G(t,x) の関数形だけを与えるようにして、拘束力(5)の計算に必要な G の微分を、全てプログラム中で数値的に計算してしまうことができる。
微分を数値的に計算するには、第1章の【1.2-注1】で述べた、速度や加速度を測定する時の方法が使える。例えば、任意の関数 G(t,x)t で微分したものは、 δt を可能な限り小さくとって、以下のように計算できる:(第1式は偏微分、第2式は全微分) GtG(t+δt,x)G(t,x)δtdGdtG(t+δt,x+˙xδt)G(t,x)δt 第2式の全微分は、 x に含まれる t も考慮するので、赤字部分が追加されている。2階微分、例えば 2tG を求めるには、この式の第1式において GtG で置き換えたものを使えばよい。ただし、丸め誤差(第1章の【1.3-注1】)のため、現実的には δt をいくらでも小さくできるわけではなく、多少誤差が増える。

差分化した運動方程式:式(9)

第3章から第5章までで扱った衝突の場合、拘束条件 G は出てきたが、2階微分は必要なかった。今考えている拘束された運動も、衝突の特別な場合(=拘束面に沿ってなめらかに衝突が起きている)とみなせるので、2階微分が現れないようにできるのではないだろうか。
衝突の場合、加速度は考えずに、衝突前後の速度の変化を直接考えた(一瞬で速度が変化するので加速度を考えることがそもそもできなかった)。同じ形にするために、速度 ˙x の1次近似を考える:(右辺の各量は時刻 t でのもの) ˙x(t+δt)˙x+¨xδt 右辺に運動方程式 m¨x=f+fc (を m で割ったもの)を代入して ¨x を消去すると、以下を得る: ˙x(t+δt)˙x+f+fcmδt(8) この式から拘束力 fc を消去すると、以下の【7.3-注1】の式(9)のように、 G の2階微分が現れない形になる。

【7.3-注1】差分化した拘束系の運動方程式:式(9)

拘束条件 G(t,x) が課されている時、運動方程式 m¨x=f+fc は、以下のような差分方程式として書くことができる:{}t+δt での値でああり、それ以外は t での値) {˙x}{P}(˙x+m1fδt)+{vwall}P1ˆG(ˆG)TvwalltG|G|ˆG(9) 右辺の {} の量は、 x(t+δt)x(t)+˙x(t)δt を用いて計算できる。
上式の P は拘束面への正射影行列である。 vwall は拘束面の速度であり、動く壁との衝突(第4章)で与えたものと同じである。式(9)は、幾何学的に解釈することができる。実際、右辺の P を含む項は拘束面に平行な方向の速度成分、同第2項は垂直な方向の成分である。特に、第1項は、拘束力を無視して δt だけ進めた時の速度を、拘束面上に正射影したものとなっている。

導出

式(8)の両辺に {P} を左乗すると、拘束力 fc が消える: {P}()={1ˆG(ˆG)T}{˙x}={˙x}+{tG|G|ˆG}(G)T˙x=tG{P}()={P}(˙x+f+fcmδt){P}(˙x+m1fδt+0)⎪ ⎪⎪ ⎪{P}fcδtPfcδt=PλGδt=0 これを合わせると式(9)そのものとなる。