$\def\Ld{{\cal L}}$ $\def\pd{\partial}$ $\def\d{\delta}$ $\def\e{\epsilon}$ $\def\l{\lambda}$ $\def\mvec{\boldsymbol}$ $\def\ds{\displaystyle}$ $\def\ipro{\!\cdot\!}$
あもんノート TOP

数値計算 for JS 実行

JavaScript をオンラインで実行するアプリである JS 実行 では、 物理定数に関するオブジェクト Phys、 数値計算に関するオブジェクト Math2 を提供しています。 ここではこれらオブジェクトの使用例をいくつか紹介します。
JS 実行

黒体放射

温度 $T$ の黒体の放射強度は、$c=\mu_0=\hbar=k_B=1$ の自然単位系で、 $$ I = {\pi^2 \over 60}\, T^4 $$ で与えられます( 統計力学 の章参照)。 よって $500 {}^\circ{\rm C}$ の黒体の放射強度の値を ${\rm W}/{\rm m}^2$ (ワット毎平方メートル)の単位で求めるコードは以下のようになります。

Phys.set()
T = ( 273.15 + 500 ) * K
I = PI^2 / 60 * T^4
I / ( W / m^2 )  ///
結果は $ I \sim 2.0261 \times 10^4 \,\,{\rm W}/{\rm m}^2 $ であることがわかります。 Phys.set() は Phys のメンバ変数全てをグローバル領域( window オブジェクト配下)にコピーします。 これにより単位込みの計算を簡単に行えます。 また、式の後に /// を付けるとその式の値が改行出力されます。

単振り子の周期

単振り子の周期 $T$ は、振幅 (最大振れ角) を $\Theta$ としたとき、 $$ T = {2T_0 \over \pi}\, I(\Theta), \ \quad T_0 = 2\pi \sqrt{l/g}, \ \quad I(\Theta) = \int_0^{\pi/2} { d\psi \over \sqrt{ 1 - \sin^2 (\Theta/2) \cos^2 \psi } } $$ で与えられます( ニュートン力学 の章参照)。 $l$ は棒の長さ、$g$ は重力加速度です。 よって $T/T_0$ のグラフを出力するコードは以下のようになります。

Math2.set()
T = function( Theta ){
	var f = function( psi ){
		return ( 1 - sin( Theta/2 )^2 * cos( psi )^2 )^(-0.5)
	}
	return 2/PI * integration( f, 0, PI/2, 200 )
}

c = Coord
c.range( 0, PI, 0, 3 )
c.axis( PI/18, 0.1 )
c.line( 0, 1, PI, 1 )
c.graph( T )
目盛は、横方向を $\pi/18=10^\circ$ きざみ、縦方向を $0.1$ きざみとしています。 Math2.set() は Math のメンバの一部( PI や sin )、 および Math2 のメンバをグローバル領域にコピーします。 Math2.integration は台形法による数値積分を与えます。

シュヴァルツシルト時空における粒子の軌道

シュヴァルツシルト時空 $(t,r,\theta,\phi)$ における粒子の軌道は、$\theta = \pi/2$ の平面内において、 $$ \frac{d^2u}{d\phi^2} = \frac{D}{2} -u + \frac{3a}{2}\, u^2, \quad u = \frac{1}{r} $$ という非線形微分方程式で与えられます。 ここで $a$ はシュヴァルツシルト半径、 $D$ は運動の定数です。 $D$ は、近日点(または遠日点)の動径座標を $R$、 このときの固有速度を $U$ とすれば、 $ D = a/(R^2 U^2) $ で与えられます( 一般相対論 の章参照)。

上の微分方程式は、 $x_0 = \phi$, $x_1 = u$, $x_2 = du/d\phi$ という変数を導入すれば、 $a=1$ の単位系において、 $$ \dot{x}_0 = 1, \quad \dot{x}_1 = x_2, \quad \dot{x}_2 = \frac{D}{2} - x_1 + \frac{3}{2}\, x_1^2 $$ と書けます。ここでドットは $x_0$ による微分を意味します。 よって軌道を描くコードは以下のようになります。

Math2.set()
R = 5 // 初期半径 / シュヴァルツシルト半径 
U = 0.42 // 初期固有速度 ( 光の場合は Infinity )
n = 3 // 周回数 
h = 0.01 // ステップ幅 

der = function( x, i ){ 
 if ( i == 0 ) return 1 
 if ( i == 1 ) return x[2] 
 if ( i == 2 ) return D/2 - x[1] + 3/2 * x[1]^2 
} 
D = 1 / ( R^2 * U^2 )
N = 2 * PI * n / h 
x = evolution( [ 0, 1/R, 0 ], der, h, N ) 

c = Coord
c.range( -R*3, R*3 ) 
c.axis( 1 )
c.circle( 0, 0, 1 ) 
for ( i = 0; i <= N; i++ ){ 
 phi = x[0][i]; r = 1 / x[1][i] 
 if ( r < 1 || r > 1e4 * R ) break 
 c.pset( r * cos( phi ), r * sin( phi ) ) 
}  
是非 JS 実行で実行し、試してみてください。 パラメータを色々と変えると面白いです。 Math2.evolution はルンゲクッタ法による微分方程式の数値解を与えます。

二重振り子のアニメーション

二重振り子のラグランジアンは、 $$ L = \frac{ml^2}{2} \left( 2\dot\theta^2 + \dot\phi^2 + 2\cos(\theta-\phi) \dot\theta \dot\phi \right) + mgl \left( 2\cos\theta + \cos\phi \right) $$ でした( 解析力学 の章参照)。 $g=l=1$ とする単位系では、ラグランジュ方程式は、 $$ \begin{pmatrix} \ddot\theta \\ \ddot\phi \end{pmatrix} = \frac{1}{2-c^2}\, \begin{pmatrix} 1 & -c \\ -c & 2 \end{pmatrix} \begin{pmatrix} -s\, \dot\phi^2 - 2\sin\theta \\ s\, \dot\theta^2 - \sin\phi \end{pmatrix}, $$$$ s = \sin( \theta - \phi ), \quad c = \cos( \theta - \phi ) $$ と整理されるでしょう。 $z_0 = t$, $z_1 = \theta$, $z_2 = \phi$, $z_3 = \dot\theta$, $z_4 = \dot\phi$ という変数を導入すれば、 $$ \dot z_0 = 1, \quad \dot z_1 = z_3, \quad \dot z_2 = z_4, $$$$ \begin{pmatrix} \dot z_3 \\ \dot z_4 \end{pmatrix} = \frac{1}{2-c^2}\, \begin{pmatrix} 1 & -c \\ -c & 2 \end{pmatrix} \begin{pmatrix} -s\, z_4^2 - 2\sin z_1 \\ s\, z_3^2 - \sin z_2 \end{pmatrix}, $$$$ s = \sin( z_1 - z_2 ), \quad c = \cos( z_1 - z_2 ) $$ と書けることに注意して、 二重振り子のアニメーションのコードは以下のようになります。

Math2.set()
T = 40   // 時間最大値
h = 0.001   // 時間ステップ幅
th0 = (3/4) * PI   // θ 初期値
phi0 = (3/4) * PI   // φ 初期値

der = function( z, i ){
 var s = sin( z[1] - z[2] ), c = cos( z[1] - z[2] )
 var A = -s * z[4]^2 - 2 * sin( z[1] )
 var B = s * z[3]^2 - sin( z[2] )
 if ( i == 0 ) return 1
 if ( i == 1 ) return z[3]
 if ( i == 2 ) return z[4]
 if ( i == 3 ) return ( A - c * B ) / ( 2 - c^2 )
 if ( i == 4 ) return ( -c * A + 2 * B ) / ( 2 - c^2 )
}
z = evolution( [ 0, th0, phi0, 0, 0 ], der, h, T/h )

putAnimation = function(){
 var i = round( tCount / h )
 drawPendulum( z[1][i], z[2][i] )
 cls(); print( "t = " + floor( tCount ) )
 tCount = Number( ( tCount + 0.1 ).toFixed( 1 ) )
 if( tCount > T ) clearInterval( iid )
}

drawPendulum = function( th, phi ){
 var x1 = sin( th ), y1 = - cos( th )
 var x2 = sin( th ) + sin( phi )
 var y2 = - cos( th ) - cos( phi )
 var c = Coord
 c.clear(); c.axis( 1 )
 c.line( 0, 0, x1, y1 ); c.line( x1, y1, x2, y2 )
 c.circle( x1, y1, 0.1 ); c.circle( x2, y2, 0.1 )
}

Coord.range( -2.2, 2.2 )
if( typeof iid != "undefined" ) clearInterval( iid )
iid = setInterval( "putAnimation()", 100 )
tCount = 0
解析力学の章では、振れ幅が小さいとし、線形近似を行いましたが、 振れ幅が大きいときの運動は非常に複雑になることが確かめられます。
inserted by FC2 system