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

数値計算 for JS 実行

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

黒体放射

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

Phys.set()
T = degC( 500 )
I = pi**2 / 60 * T**4
Text.putlnN( I / ( W / m**2 ) )
Phys.set() は Phys のメンバをグローバル領域( window オブジェクト配下)にコピーします。 これにより単位込みの計算を簡単に行えます。

誤差関数

誤差関数 ${\rm erf}$ は、 $$ {\rm erf}(x) = {2 \over \sqrt{\pi}} \int_0^x dt \, \exp( -t^2 ) $$ で与えられます( 関数論と応用数学 の章参照)。 これを数値計算するコードは以下のようになります。

pi = Math.PI

erf = function( x ){
	var f = function( t ){
		return Math.exp( -1 * t**2 )
	}
	return Numerical.integral( f, 0, x, 5000 ) * 2 / pi**0.5
}

Text.clear()
Text.putlnN( erf( 1.96 / 2**0.5 ) )
Numerical.integral はシンプソン法による数値積分値を与えます。 第4引数の 5000 は積分区間の分割数です。

単振り子の周期

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

pi = Math.PI

J = function( Th ){
  var f = function( psi ){
    var v = 1 - Math.sin( Th/2 )**2 * Math.cos( psi )**2
    return 1 / v**0.5
  }
  return 2 / pi * Numerical.integral( f, 0, pi / 2, 200 )
}

c = Coordinates
c.range( 0, pi, 0, 3 )
c.axis( pi / 18, 0.1 )
c.line( 0, 1, pi, 1 )
c.graph( J )
Coordinates は論理座標のオブジェクト。 目盛は、横方向を $\pi/18=10^\circ$ きざみ、縦方向を $0.1$ きざみとしています。

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

シュヴァルツシルト時空 $(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 = Numerical.evolution( [ 0, 1/R, 0 ], der, h, N )

c = Coordinates
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.point( r * Math.cos( phi ), r * Math.sin( phi ) )
}
是非 JS 実行で実行し、試してみてください。 パラメータを色々と変えると面白いです。 Numerical.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 ) $$ と書けることに注意して、 二重振り子のアニメーションのコードは以下のようになります。

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

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

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

drawPendulum = function( th, phi ){
  var m = Math, c = Coordinates
  var x1 = m.sin( th ), y1 = - m.cos( th )
  var x2 = m.sin( th ) + m.sin( phi )
  var y2 = - m.cos( th ) - m.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 )
}

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

SU(3)リー代数の構造定数

SU(3)のリー代数の構造定数は、 ゲルマン行列を $\l_a$ ($a=1 \sim 8$) として、 $$ f_{abc} = \frac{-i}{4}\, {\rm tr} \left( [\l_a, \l_b] \l_c \right) $$ で与えられます(連続群論入門 の章参照)。 よって、複素数・複素行列のクラスを用いると、 構造定数の数値計算は以下のようになります。

iu = new Complex( 0, 1 )
miu = iu.times( -1 )

gmn = []
gmn[ 1 ] = CMatrix.byArray( [[0,1,0],[1,0,0],[0,0,0]] )
gmn[ 2 ] = CMatrix.byArray( [[0,miu,0],[iu,0,0],[0,0,0]] )
gmn[ 3 ] = CMatrix.byArray( [[1,0,0],[0,-1,0],[0,0,0]] )
gmn[ 4 ] = CMatrix.byArray( [[0,0,1],[0,0,0],[1,0,0]] )
gmn[ 5 ] = CMatrix.byArray( [[0,0,miu],[0,0,0],[iu,0,0]] )
gmn[ 6 ] = CMatrix.byArray( [[0,0,0],[0,0,1],[0,1,0]] )
gmn[ 7 ] = CMatrix.byArray( [[0,0,0],[0,0,miu],[0,iu,0]] )
gmn[ 8 ] = CMatrix.byArray( [[1,0,0],[0,1,0],[0,0,-2]] )
	.times( 3 ** ( -1/2 ) )

Text.clear()
for( a = 1; a <= 8; a++ )
for( b = a + 1; b <= 8; b++ )
for( c = b + 1; c <= 8; c++ ){
	f = gmn[ a ].commute( gmn[ b ] ).times( gmn[ c ] )
		.trace().times( miu ).times( 1/4 )
	if( f.abs() > 0 ) Text.putln( "" + a + b + c + " -> " + f )
}
構造定数 $f_{abc}$ のうち、$a \lt b \lt c$ を満たす代表元で、 特に $0$ でないものだけ表示しています。 iu は虚数単位、gmn[ a ] はゲルマン行列 $\l_a$ を意味します。

ガンマ関数とゼータ関数

解析接続されたガンマ関数は、 $$ \Gamma(z) = \begin{cases} \ \ds \sum_{k=0}^\infty \, { (-1)^k \over k!\, (z+k) } + \int_1^\infty dt \ t^{z-1} e^{-t} \quad \left( {\rm Re}\, z \ge {1 \over 2} \right) \\ \ \ \ds {\pi \over \sin(\pi z) \Gamma(1-z)} \ \ \, \quad \left( {\rm Re}\, z < {1 \over 2} \right) \end{cases} $$ ゼータ関数は、 $$ \zeta(z) = \Gamma(1-z) \left( \int_0^{2\pi} {d\th \over 2\pi}\, {e^{iz(\th-\pi)} \over e^{\ds e^{i\th}}-1} +{\sin(\pi z) \over \pi}\, \int_1^\infty \! dt \, {t^{z-1} \over e^t -1} \right) $$ で与えられるのでした (ゼータ関数 の章参照)。 よって複素数のクラスを用いると、 これらを数値計算するコードは以下のようになります。

one = new Complex( 1 )
pi = one.times( Math.PI )
iu = new Complex( 0, 1 )
miu = iu.times( -1 )
Nm = Numerical
M = Math

cpow = function( a, z ){
	if( typeof z == "number" ) return one.times( a ** z )
	return z.times( M.log( a ) ).exp()
}

gamma = function( z ){
	if( typeof z == "number" ) z = one.times( z )
	if( z.re < 1/2 ){
		var g = gamma( one.sub( z ) )
		return pi.div( pi.times( z ).sin() ).div( g )
	}
	var L = 30, d = 200
	var f1 = function( t ){
		return cpow( t, z ).times( M.exp( -t ) / t )
	}
	var f2 = function( k ){
		var a = z.add( k ).times( Nm.factorial( k ) )
		return one.times( (-1) ** k ).div( a )
	}
	return Nm.integral( f1, 1, L, d )
		.add( Nm.summation( f2, 0, 10 ) ) 
}

zeta = function( z ){
	if( typeof z == "number" ){
		if( z > 0 && M.abs( z - M.round( z ) ) < 1e-10 )
			z += 1e-9
		z = one.times( z )
	}
	var L = 30, d = M.floor( 100 + z.abs() * 300 )
	var f1 = function( t ){
		return cpow( t, z ).div( t * ( M.exp( t ) - 1 ) )
	}
	var f2 = function( th ){
		var a = miu.times( z ).times( pi.sub( th ) ).exp()
		var b = iu.times( th ).exp().exp().sub( 1 )
		return a.div( b )
	}
	var i1 = Nm.integral( f1, 1, L, d )
		.times( pi.times( z ).sin() ).div( M.PI )
	var i2 = Nm.integral( f2, 0, 2 * M.PI, d ).div( 2 * M.PI )
	return i1.add( i2 ).times( gamma( one.sub( z ) ) )
}

c = Coordinates
c.range( -5, 5 )
c.axis( 1 )
c.graph( function( x ){ return zeta( x ).re } )
cpow は正の実数の複素数乗を与える関数です。 確認のため、 実軸上のゼータ関数のグラフを出力しています。

電子・陽電子の対消滅

最後に、 電子・陽電子の対消滅の散乱振幅や断面積を、 ヘリシティ固定で数値計算するコードを載せておきます (素粒子論の計算 の章参照)。

M = Math
CM = CMatrix
iu = new Complex( 0, 1 )
miu = iu.times( -1 )

// ガンマ行列(下付き添字)
gm = [] 
gm[0] = CM.byArray( [[0,0,1,0],[0,0,0,1],[1,0,0,0],[0,1,0,0]] )
gm[1] = CM.byArray( [[0,0,0,1],[0,0,1,0],[0,-1,0,0],[-1,0,0,0]] )
gm[2] = CM.byArray( [[0,0,0,-1],[0,0,1,0],[0,1,0,0],[-1,0,0,0]] )
	.times( iu )
gm[3] = CM.byArray( [[0,0,1,0],[0,0,0,-1],[-1,0,0,0],[0,1,0,0]] )
gm[5] = CM.byArray( [[1,0,0,0],[0,1,0,0],[0,0,-1,0],[0,0,0,-1]] )

// ディラック粒子の波動関数
// 引数 s はヘリシティ符号、k は運動量の大きさ、th は天頂角
Dirac_u = function( s, k, th ){
	var k_0 = M.sqrt( k ** 2 + 1 )
	var p = M.sqrt( k_0 + k ), m = M.sqrt( k_0 - k )
	var cs = M.cos( th / 2 ), sn = M.sin( th / 2 )
	if( s > 0 )
		return CM.byArray( [[cs*p],[sn*p],[cs*m],[sn*m]] )
	return CM.byArray( [[-sn*m],[cs*m],[-sn*p],[cs*p]] )
		.times( iu )
}
Dirac_v = function( s, k, th ){
	return gm[2].times( Dirac_u( s, k, th ).conj() ).times( miu )
}

// 光子の波動関数
// 引数 mu は時空添字、s はヘリシティ符号、th は天頂角
photon = function( mu, s, th ){
	var s2 = M.sqrt( 2 )
	if( mu == 0 ) return 0
	if( mu == 1 ) return M.cos( th ) / s2
	if( mu == 2 ) return iu.times( s ).div( s2 )
	if( mu == 3 ) return -M.sin( th ) / s2
}
photon_sla = function( s, th ){
	var a = new CM( 4, 4 )
	for( var mu = 1; mu <= 3; mu++ )
		a = a.add( gm[ mu ].times( photon( mu, s, th ) ) )
	return a
}

// 以下、対消滅の直接グラフ、散乱振幅、微分断面積、断面積
// 引数 p は電子・陽電子の運動量の大きさ、th は散乱角
//  s, sd は電子・陽電子のヘリシティ符号、l, ld は光子のそれ

graph_A = function( p, th, s, sd, l, ld ){
	var E = M.sqrt( p ** 2 + 1 )
	var p_sla = gm[ 0 ].times( E ).add( gm[ 3 ].times( p ) )
	var q_sla = gm[ 0 ] .add( gm[ 1 ].times( M.sin( th ) ) )
		.add( gm[ 3 ].times( M.cos( th ) ) ).times( E )
	var a = Dirac_v( sd, p, M.PI ).hermit().times( gm[ 0 ] )
	a = a.times( photon_sla( -ld, th + M.PI ) )
	a = a.times( p_sla.sub( q_sla ).add( CM.identity( 4 ) ) )
	a = a.times( photon_sla( -l, th ) )
	a = a.times( Dirac_u( s, p, 0 ) )
	return a.trace().times( iu ).times( 2 * M.PI )
		.div( E ** 2 - E * p * M.cos( th ) )
}
s_amp = function( p, th, s, sd, l, ld ){
	var a = graph_A( p, th, s, sd, l, ld )
	var b = graph_A( p, th + M.PI, s, sd, ld, l )
	return a.add( b )
}
dif_cs = function( p, th, s, sd, l, ld ){
	var a = s_amp( p, th, s, sd, l, ld ).abs() ** 2
	return a / ( 256 * M.PI**2 * p * M.sqrt( p ** 2 + 1 ) ) 
}
c_section = function( p, s, sd, l, ld ){
	var th_max = M.PI / 2 ; if( l != ld ) th_max *= 2
	var f = function( th ){
		return M.sin( th ) * dif_cs( p, th, s, sd, l, ld )
	}
	return Numerical.integral( f, 0, th_max, 20 ) * 2 * M.PI
}

// 確認グラフ出力
c = Coordinates
c.range( 0, 5 )
c.axis( 1 )
c.graph( function( p ){ return c_section( p, 1, 1, 1, 1 ) } )

inserted by FC2 system