$\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 実行 では、 数値計算に関するオブジェクト Math2、 物理定数に関するオブジェクト Phys を提供しています。 ここではこれらオブジェクトの使用例をいくつか紹介します。
JS 実行

黒体放射

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

pi = Math.PI   // 円周率
K = Phys.K   // ケルビン
W = Phys.W   // ワット
m = Phys.m   // メートル

T = ( 273.15 + 500 ) * K
I = pi**2 / 60 * T**4
Text.putlnE( I / ( W / m**2 ) )   // 2.0261269e+4
結果は $ I \sim 2.0261269 \times 10^4 \,\,{\rm W}/{\rm m}^2 $ であることがわかります。 Phys オブジェクトの利用により、このように単位込みの計算が可能です。 Text は JS 実行が用意しているテキスト入出力に関するオブジェクトで、 Text.putlnE は引数の値を小数点以下 7 桁の指数表記で行出力します。

単振り子の周期

単振り子の周期 $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$ のグラフを出力するコードは以下のようになります。

pi = Math.PI; sin = Math.sin; cos = Math.cos
T = function( Theta ){
	var f = function( psi ){
		return ( 1 - sin( Theta / 2 )**2 * cos( psi )**2 )**-0.5
	}
	return 2 / pi * Math2.integ( 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.integ により台形法による数値積分を行っていることに注意してください。

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

シュヴァルツシルト時空 $(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$ による微分を意味します。 よって軌道を描くコードは以下のようになります。

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 * Math.PI * n / h 
x = Math2.evol( [ 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 * Math.cos( phi ), r * Math.sin( phi ) ) 
} 
是非 JS 実行で実行し、試してみてください。 パラメータを色々と変えると面白いです。 Math2.evol によりルンゲクッタ法による微分方程式の数値解を取得しています。

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

二重振り子のラグランジアンは、 $$ 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 ) $$ と書けることに注意して、 二重振り子のアニメーションのコードは以下のようになります。

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

der = function( z, i ){
 var s = sin( z[1] - z[2] ), c = cos( z[1] - z[2] ),
 A = - s*z[4]**2 - 2*sin( z[1] ), 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 )
}
sin = Math.sin; cos = Math.cos
z = Math2.evol( [ 0, th0, phi0, 0, 0 ], der, h, T / h )

putAnimation = function(){
 var i = Math.round( tCount / h )
 drawPendulum( z[1][i], z[2][i] )
 cls(); print( "t = " + Math.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),
 x2 = sin(th) + sin(phi), y2 = - cos(th) - cos(phi)
 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 )
}

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