束縛条件下の運動 – 自由度がうまく落とせない運動

このページは
束縛条件下の運動 – ホロノミックな束縛と非保存力
の続きです。

拘束がある場合の2次元の運動


ラグランジュの方程式を解く際に適切な座標変換が見付からず自由度が落とせない場合を考えます。その時の運動方程式は

の形で止まります。式(75)において、未知の関数は\(x(t), y(t), \lambda(t)\)の3つであり、条件式は3つなので解くことは出来るはずです。
方針は、式(75)から\(\lambda(t)\)を消去すれば良いのです。

束縛条件\(f(x,y,t)=0\)が成立する点\(x=x_0, y=y_0\)周りで、微小時間\(\Delta t\)の変化を考えます。
点\((x_0, y_0)\)まわりでテーラー展開すれば、

を得ます。よって

が満たされる変化でなければなりません。すなわち、\(x,y\)の時間変化は式(77)をいつも満たしていなければならないのです。両辺を\(\Delta t\)で割って極限をとれば、

が成立することが分かります。式(75)を書きかえれば、

を得ます。式(79c)もまた時間変化しても成立し続けなければなりません。よって、左辺の時間微分を取れば、

を得ます。整理して

を得ます。ここで関係式

を用いました。式(79a), (79b)から未定乗数\(\lambda(t)\)を消去すれば

を得ます。よって、式(81)と式(84)から連立方程式

を立てることが出来ます。これをあらわに解けば、運動方程式

を得ます。ここで、式(86), (87)の右辺の\(-\frac{\partial V}{\partial x}, -\frac{\partial V}{\partial y}\)は\(x, y\)方向の保存力であり、右辺の残りの項は束縛条件によって決まる束縛力であることが分かります。

拘束がある場合の3次元の運動


スカラー関数\(f\)で表される拘束がある条件のもとで, 運動方程式

またはベクトル表記で

について考えます。ここで表記を略すために

の略記を用いました(\(y, z\)方向についても同様)。
また、\(\nabla\)はナブラ演算子で、

です。2次元の場合と同様に、式(88d)の時間微分から

が導けます。ここで、\(\hat{H}f\)は関数\(f\)のヘッセ行列であり、

を表します。式(92)は、式(81)を3次元に拡張した形になっています。

式(88a),(88b),(88c)から未定乗数を消去すると従属な3つの関係式

を得ることが出来ます。この3式で独立なのは2つです。
独立な式として式(94a)と式(94b)を選ぶことにし、式(92)と共に書けば、連立方程式

を得ます。この式をあらわに解けば、運動方程式

を得ることが出来ます。ここで、\(n_x, n_y, n_z\)は壁の法線ベクトル\(\mathbf{n}=\frac{\nabla f}{|\nabla f|}\)の\(x,y,z\)成分です。偏微分ではないことに注意してください。

ベクトル表記であれば

と書くことが出来ます。

単振り子の例


ここでは単振り子を例に挙げ、式(86)に代入した時本当に振り子を表す方程式になっているのか、また束縛条件の選び方に依らず、同じ運動を記述しているのかを例に挙げます。
保存力は重力で

を考えます。

考える束縛条件は
\(
f(x,y,t)=x^2+y^2-l^2=0
\)


\(
f(x,y,t)=\sqrt{x^2+y^2}-l=0
\)

です。

束縛条件1


束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

束縛条件2


束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

考察


同じ束縛条件なので、式(105)と(109)の示す運動は同じになるはずです。
それを示すには、2式の差異である

を示せればよい、ということです。そのために拘束力に対して座標変換

を考えます。右辺と左辺をそれぞれ計算すれば、

となり、同じ運動を記述していることが分かります。

Fortran90によるプログラム


Fortran90で2次元平面を動く質点の束縛運動のプログラムを作ります。
そのプログラムは
https://slpr.sakura.ne.jp/qp/supplement_data/constraint.tar.gz
に置いてあります。

一応補足しておきますが、初期状態の位置は\(f(x, y, t)=0\)を満たすx,yでなければなりませんし、初期速度も\(f(x, y, t)=0\)の傾きと一致していなければなりません。
上の条件を満たさずとも何も警告は出さず、計算自体は行われますが、その計算は束縛運動ではなくなります。

上のプログラムの中で、変更すると思われるのは以下の2つです。
1. fp2d
2. fw2d

1. fp2dは、ポテンシャル\(V(x,y)\)を記述しており、デフォルトのプログラムでは、重力による位置エネルギー
\(
V(x,y)=mgy
\)
に設定しています。

2. fw2dは、束縛条件\(f(x,y, t)=0\)を記述しており、デフォルトのプログラムでは、
\(
f(x,y,t)=y – [2(cos(1.5*x)-1) + 0.1x^2]=0
\)
に設定しています。

もしも非保存力を入れたい場合は、grkの最後の方を変更してください。

数値計算では、
時間発展は刻み幅制御の陽的ルンゲクッタ、
gradient, Hesse行列、束縛条件に関する時間の二階微分は、Hyper-dual numbersを利用してほぼ厳密な値を得ています。

program main
  ! Author : sikino
  ! Date   : 2019/10/12 (yyyy/mm/dd)
  ! URL    : http://slpr.sakura.ne.jp/qp/constraint-condition3/
  use Parameters
  use Hyperdualmod
  implicit none
  integer::i,N,info,Nt
  double precision::t,h,tol,ta,tb
  double precision,allocatable::x(:),tgrid(:)
  external::grk
  double precision::E,Ek,Ep

  N=4
  allocate(x(1:N))

  ! Time range [ta, tb]
  ta = 0d0
  tb = 30d0
  ! Divide time range equally among Nt
  Nt = 201
  ! Tolerance of the RK method
  tol=1d-8

  ! Initial condition
  x(1)=0d0 ! x (ta)
  x(2)=4d0 ! x'(ta)
  x(3)=0d0 ! y (ta)
  x(4)=0d0 ! y'(ta)

  ! Time grid
  allocate(tgrid(1:Nt))
  do i=1,Nt
     tgrid(i) = (i-1)*(tb-ta)/dble(Nt-1)+ta
  enddo
 
  do i=2,Nt
     info = 0
     t = tgrid(i-1)
     h = (tgrid(i) - t)*0.1d0
     do while(info.le.0)
        call drkf45(grk,t,h,N,x,tgrid(i),info,tol)
     enddo
     
     ! Kinetic energy
     Ek = 0.5d0*mass*(x(2)**2+x(4)**2)
     ! Potential energy
     Ep = mass*g*x(3)
     ! Total energy
     E  = Ek + Ep
     write(10,'(8e25.10e3)')t,x(1),x(2),x(3),x(4),E,Ek,Ep
  enddo
 
  stop
end program main

subroutine fp2d(N,xH,fH)
  use Parameters, only:mass, g
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  type(Hyperdual),intent(in)::xH(1:N)
  type(Hyperdual),intent(out)::fH

  ! Potential

  type(Hyperdual)::x,y,t
  x = xH(1)
  y = xH(2)
  t = xH(3)
  ! If conservative force,
  !  fH doesn't depend on the time.
  !-----------------------------
 
  fH = mass * g * y
 
  return
end subroutine fp2d

subroutine fw2d(N,xH,fH)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  type(Hyperdual),intent(in)::xH(1:N)
  type(Hyperdual),intent(out)::fH

  ! Wall of the function

  type(Hyperdual)::x,y,t
  x = xH(1)
  y = xH(2)
  t = xH(3)
  !-----------------

  fH = y - (2d0*(cos(1.5d0*x)-1d0) + 0.1d0*x*x)
!  fH = (y-5d0)**2 + x**2 - 25d0
 
  return
end subroutine fw2d

デフォルトのまま計算すると、fort.10に以下の出力が現れます。アニメーションは、同じファイルに入っているanime.pltで実行できます。

エネルギーは計算の範囲ではほとんどありません。全力学的エネルギーは、計算時間内でおおよそ5桁目が変化していました。

参考文献


6. 振り子 -物理学の見つけ方
9. 自由な座標 -物理学の見つけ方

束縛条件下の運動 – ホロノミックな束縛と非保存力

このページは
束縛条件下の運動 – ホロノミックな束縛 1
の続きです。

単振り子



ポテンシャルが存在し、ホロノミックな束縛が課されている代表的な問題である、単振り子を考えましょう。
単振り子は、束縛条件下の運動 – ホロノミックな束縛 1で考えた円に束縛されている問題に対してポテンシャルが加わった問題です。

座標系は以下のような極座標系を選びます。

通常の極座標系とは異なります。\(\theta=0\)で鉛直下方向にとる座標系にしたかったので、通常の極座標系とは変えています。

ホロノミックな束縛条件\(f(x,y,t)=0\)を含めたラグランジアンは、デカルト座標系から極座標系に変えて

と書くことが出来ます。注記しますが、ラグランジュの未定乗数は時間に依存していて構いません。
ラグランジアンの偏微分をそれぞれ計算すれば、

ですので、運動方程式と束縛条件

が得られます。連立微分方程式(31)の3本の方程式に対して未知の関数は\(r(t),\theta(t),\lambda(t)\)の3つなので、方程式を立てて\(r(t),\theta(t)\)を求めることが出来ます。

では、伸び縮みしない振り子を考えましょう。束縛条件は

をして表現することが出来ます。束縛条件の一般化座標による偏微分は

ですので、運動方程式

を得ます。式(34)の第2式を変形すると

ですので、これはまさに振り子を表す運動方程式です。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
 \(m=3, l=1, g=9.8\)
\(\theta(0)=\frac{4\pi}{5}, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

支点が動く振り子



次に、束縛条件が2つある場合を考えましょう。その例題として、
振り子の支点が摩擦の無いレールの上に載っていることを考えます。
ラグランジアンは

です。2つのホロノミックな束縛条件

を考えます。すると、ホロノミックな束縛条件を含めたラグランジアンは、

と書くことが出来ます。
これから伸び縮みしない振り子を考えることを見越して、座標系

を考えます。この新たな座標系でのラグランジアンは、

と書けます。それぞれの偏微分である



を用いると、運動方程式

を得ます。
連立微分方程式(46)の4本の方程式と2本の束縛条件に対して、未知の関数は\(x(t), y(t), r(t), \theta(t), \lambda_0(t), \lambda_1(t)\)の6つなので、方程式を立ててそれぞれの関数を求めることが出来ます。

では、支点が直線状のレールの上に置かれている状況を考えましょう。
支点に関する束縛条件\(f_0\)と, 振り子の紐の長さが変わらないという束縛条件\(f_1\)を考えると

と書くことが出来ます。なので、

とそれぞれ計算できるので、運動方程式

を得ます。求めたい未知の関数は\(x(t), \theta(t)\)です。
式(51a)と(51b)から、\(x(t), \theta(t)\)に関する運動方程式

を得ます。変形すると

ですので、あらわに求めれば、

を得ます。
特に支点の質量が振り子の質量に比べて非常に重い場合、すなわち\(m_0 \gg m_1\)の場合、

になります。第2式はまさに単振り子の方程式となっています。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
 \(m_0=3, m_1=1, l=1, g=9.8\)
\(x(0)=0, x'(0)=0, \theta(0)=\frac{4\pi}{5}, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

支点が放物線上に束縛されている振り子



座標系は1つ前の問題と同じく、

にとることにします。すると運動方程式は式(46)で表されるので、運動方程式の導出過程を省略することが出来ます。
問題設定から、束縛条件は

と書くことが出来ます。なので、それぞれの偏微分は

と書けますので、運動方程式は

です。式(60a)と(60b)から\(\lambda_0(t)\)を消去すると、

を得ます。更に、束縛条件から座標\(y(t)\)の時間変化は\(x(t)\)を用いて

と書くことが出来ます。式(61)と式(60d)から、\(x(t), \theta(t)\)の運動方程式

を得ます。あらわに解くと式が長くなってしまうので、ここで止めておきます。

実際に解きますと、

というような振る舞いになります。

・パラメータ
 \(m_0=3, m_1=1, l=1, g=9.8\)
 \(x(0)=0, x'(0)=5, \theta(0)=0, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

単振り子と空気抵抗


束縛条件と、非保存力である空気抵抗が存在する場合を考えます
座標系は

です。ホロノミックな束縛条件を含めたラグランジアンは、

と書けます。

さて、非保存力である空気抵抗が存在することを考えましょう。
非保存力なのでラグランジアンで書くことは出来ませんので、個別に求める必要があります。
空気抵抗の力\(F_x, F_y\)はデカルト座標系において

で掛ける事を既知とします。今知りたいことは、式(64)で表される一般化座標\(r, \theta\)で式(67)がどのように表されるのか?ということです。
\(r, \theta\)において、一般化座標における空気抵抗\(Q_x, Q_y\)は、

で変換することが出来ます。具体的にそれぞれの力を示せば、

ですので\(Q_x, Q_y\)は、

と表せられます。非保存力がある場合、ラグランジュの方程式の右辺にそれを入れれば良いので、

と表されます。よって運動方程式

を得ます。
長さ\(l\)の伸び縮みしない紐に繋がれた振り子をを考えれば、束縛条件は

ですので、運動方程式

を得ます。

実際に解きますと

・パラメータ
 \(m=3, l=1, g=9.8, c_1=0.3, c_2=0.3\)
 \(\theta(0)=0.8\pi, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

このページはここで終わりです。
続いて、自由度がうまく落とせない、うまい座標系が見付からない場合を考えます。

次のリンク
束縛条件下の運動 – 自由度がうまく落とせない運動

参考文献


田辺行人、品田正樹著『理・工基礎 解析力学』 裳華房

束縛条件下の運動 – ホロノミックな束縛

ホロノミックな束縛を受けている質点の運動を考えます。

ホロノミックな束縛とは、
”束縛条件が一般化座標と時間のみで決まる方程式によって解析的に表現できるもの”
です。空気抵抗のない単振り子などは、半径(一般化座標)が一定という条件のみで決まるので、ホロノミックな束縛に相当します。

ラグランジュの方程式を基本として考えていきます。

スタートはラグランジアン、ゴールは運動方程式を立てる事です。

2次元平面上に束縛されている質点


1つの質量を持った質点が、ある軌道上しか動けない運動に対する運動方程式を導きましょう。
ある軌道上しか動けないという条件は、時刻\(t\)における質点の位置\((x,y)\)が、関数\(f(x,y,t)=0\)を満たす動きしか許されないということで表現することが出来ます。

束縛条件を含めた実効的なラグランジアンは, 未定乗数\(\lambda(t)\)を用いて

と書けます。ラグランジュの方程式

に代入すれば運動方程式

を得ます。

円に束縛されている運動


運動が半径\(l\)の円に束縛されている場合、束縛条件は

として表現されます。もちろん、\(f(x,y,t)=\sqrt{x^2+y^2}-l=0\)としても構いません。運動は束縛条件の形には依らないからです。

式(3)に代入すると、運動方程式

を得ます。方程式(3),(4),(5)の3本の方程式に対して未知の関数は\(x(t),y(t),\lambda(t)\)の3つなので、方程式を立てて\(x(t),y(t)\)を求めることが出来ます。

式(5)の第1式に\(x\)、第2式に\(y\)を掛けて両式を引くと

です。変形すれば

を得ます。どういう解き方でも構わないのですが、半径\(l\)が変わらないことから極座標

を用いることを考えます。すると、式(7)は

と書けます。よって

ここで、\(\omega,\theta_0\)は初期条件から来まる定数です。よって、解

を得ます。ちなみに、\(\lambda(t)\)は

です。

角周波数\(\omega\)で等速回転する筒に入れられた質点が描く軌道



続いて、ホロノミックな束縛条件が時間\(t\)に依存する問題を考えてみましょう。
そこで、角周波数\(\omega\)で等速回転する筒に入れられた質点の軌道を考えます。

極座標表示の方が扱いやすいので、極座標で考えていきます。

ホロノミックな束縛条件を含めた実効的なラグランジアンは、デカルト座標系のラグランジアンから出発して

を得ます。ここで、\(\lambda(t)\)は未定乗数を表します。
計算していけば、運動方程式

を得ます。方程式(18a),(18b),(18c)の3本の方程式に対して未知の関数は\(r(t),\theta(t),\lambda(t)\)の3つなので、方程式を立てて\(x(t),y(t)\)を求めることが出来ます。

角周波数\(\omega\)で等速回転する筒を表現する関数を考えてみましょう。
筒の上の任意の点\((x,y)\)は、

と書くことで表現することが出来ます。よって、束縛条件から

という関係式を得ます。式(21)を運動方程式(18)に代入すれば、

を得ます。式(22a)を解けば、

であるので、解

を得ます。ここで、\(A,B\)は初期条件から決まる定数です。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
\(A=0.01, B=0, w=2\)

このページはここで終わりです。
続いて、ポテンシャルが加わる問題である単振り子と束縛条件を考えていきます。
次のリンク
束縛条件下の運動 – ホロノミックな束縛と非保存力

参考文献


田辺行人、品田正樹著『理・工基礎 解析力学』 裳華房