簡易数学公式集（その1） 12/3/10更新
ベッセル関数や楕円積分を実際に自分のパソコンで計算してみたい人は次のページへどうぞ
珍しい電卓
珍しい電卓(JavaScript版)
ベッセル関数を連分数を用いた簡単アルゴリズムで、第1種と第3種の楕円積分をガウス変換を用いた簡単アルゴリズムで普通のパソコンで計算する
普通の電卓(関数電卓)
普通の電卓(関数電卓)(JavaScript版) 珍しい電卓のお供にどうぞ
おすすめ項目索引
ベッセル関数計算プログラム
第1種楕円積分計算プログラム（ガウス変換による）
第2種楕円積分計算プログラム（21a展開による）
第3種楕円積分計算プログラム（ガウス変換による）
第1種完全楕円積分計算プログラム（ガウス変換による）
第2種完全楕円積分計算プログラム（21a展開による）
第3種完全楕円積分計算プログラム（ガウス変換による）
ヤコビの振幅関数計算プログラム
第3種完全楕円積分の母数とパラメタによる微分
二元連立方程式
a x+b y=f ;c x+d y=g のとき x=(d f-b g)/A ;y=(a g-c f)/A: (A=a d-b c)
二次元のニュートンラフソン法
f(x,y)=0 ;g(x,y)=0 のとき根の補正項は
Δx=-((∂g/∂y) f-(∂f/∂y) g)/A; Δy=-((∂f/∂x) g-(∂g/∂x) f)/A ( A=(∂f/∂x)(∂g/∂y)-(∂f/∂y)(∂g/∂x) )
級数
1/(1-x)=1+x+x^2+x^3+ …
1/√(1-x)=1+(1/2)x+(1/2)(3/4) x^2+(1/2)(3/4)(5/6) x^3+ …
√(1-x)=1-(1/2)x-(1/2)(3/4)x^2/3-(1/2)(3/4)(5/6)x^3/5- …
ベッセル関数 (アンダーバー _ は以降が下付添え字であることを示しています)
J_0(x)=1-x^2 /2/2+x^4 /2/2/4/4- …
J_1(x)=x/2(1-x^2 /2/4+x^4 /2/4/4/6- …)
J_2(x)=x^2/2/2/2(1-x^2 /2/6+x^4 /2/6/4/8- …)
J_3(x)=x^3 /2/2/2/2/3(1-x^2 /2/8+x^4/2/8/4/10- …)
…
変形されたベッセル関数
I_0(x)=1+x^2 /2/2+x^4 /2/2/4/4+ …
I_1(x)=x/2(1+x^2 /2/4+x^4 /2/4/4/6+ …)
I_2(x)=x^2/2/2/2(1+x^2 /2/6+x^4 /2/6/4/8+ …)
I_3(x)=x^3/2/2/2/2/3(1+x^2 /2/8+x^4/2/8/4/10+ …)
…
ベッセル関数の漸化式 J_n-1(x)/J_n(x)=2n/x-1/(J_n(x)/J_n+1(x)) (1)
ベッセル関数の級数 J_0(x)-2J_2(x)+2J_4(x)- …= 1 (2)
(1),(2)の性質からベッセル関数の値を計算する事が出来る。
VBによるプログラム例
Function jn(n, x) 'ベッセル関数
If n < 0 Or n > 100 Then jn = "nが不適切": Exit Function
ni = Int(n): If n - ni <> 0 Then jn ="nが不適切": Exit Function
If x < 0 Or x > 200 Then jn = "xが不適切": Exit Function
If x = 0 Then
If ni = 0 Then jn = 1 Else jn = 0
Exit Function
Else
End If
i = 2 * Int(0.6 * x + ni+ 5): a = 2 * i / x: b= 1E-50: w = (1 + 1 / a ^ 2) * b: k = 0: il = i
While i > 0
a = 2 * i / x - 1 / a: b = b * a: w = w + k * b: k = 1 - k: i = i - 1
If i = ni Then c = b
Wend
jn = c / (2 * w - b)
End Function
Function ijn(n, x) '変形されたベッセル関数
If n < 0 Or n > 100 Then ijn = "nが不適切": Exit Function
ni = Int(n): If n - ni <> 0 Then ijn = "nが不適切": Exit Function
If x < 0 Or x > 200 Then ijn = "xが不適切": Exit Function
If x = 0 Then If ni = 0 Then ijn = 1 Else ijn = 0
Exit Function
Else
End If
i = 2 * Int(0.6 * x + ni+ 5): a = 2 * i / x: b = 1E-50: w = (1 + 1 / a ^ 2) * b: k = 0
While i > 0
a = 2 * i / x + 1 / a: b = b * a: w = w + k * b: k = 1 - k: i = i - 1
If i = ni Then c = b
Wend
a = c / (2 * w - b): a = a * (Exp(x) +Exp(-x)) / 2: ijn = a
End Function
ガンマ関数 Γ(1/2) = √(π): Γ(1)=1
ガンマ関数の漸化式 Γ(x+1) = x Γ(x)
Γ(n+1)=n! : Γ(n+1/2)=((2n-1)!!/2^n)√(π)
楕円積分 (∫φ_0は0からφまでの定積分を表しています)
第1種楕円積分 E1(φ,k)=∫φ_0 1/√(1-k^2 sin^2φ) dφ=φ+(1/2) k^2 S_2(φ)+(1/2)(3/4)k^4 S_4(φ)+ …
第2種楕円積分 E2(φ,k)=∫φ_0 √(1-k^2 sin^2φ) dφ=φ+(-1/2) k^2 S_2(φ)+(-1/2)(3/4) k^4 S_4(φ)+(-1/2)(3/4)(5/6) k^6 S_6(φ)+ …
(ここで∫φ_0 sin^nφ dφ=S_n(φ)と略記)
S_n(φ)は漸化式 S_n(φ) =((n-1)S_n-2(φ) -sin^(n-1)φcosφ)/nで計算出来る
S_0(φ) =φ; S_2(φ) =(1/2)S_0(φ)-sinφcosφ/2=φ/2-sinφcosφ/2
S_4(φ) =(3/4)S_2(φ)-sin^3φcosφ/4=(3/4)(φ/2-sinφcosφ/2)-sin^3φcosφ/4=(3/8)φ-(1/8)sinφcosφ(2sin^2φ+3)
…
第1種楕円積分はkを小さくする変数変換（ランデン変換ではなくガウス変換を使う。ランデン変換は変数の変域が
2倍になってしまうのに対してガウス変換は変域の巾は変わらずただ曲線の形が左に寄るだけなので変換を繰り返して行うに都合が良い。
ガウスのほうがランデンよりも少し後の人なんですね。）
Φ= arcsin ((1+k) sinφ/(1+k sin^2φ)) (Φは古い変数、Kは古い母数K=2√(k)/(1+k)、φは新しい変数、kは変換後の母数)
∫Φ_0 1/√(1-(2√(k)/(1+k))^2 sin^2 Φ) dΦ=(1+k)∫φ_0 1/√(1-k^2 sin^2φ) dφ
をkが実用上0と見なせるまで数回繰り返して行えば被積分関数は実用上定数となり簡単に求められる。(^^
VBによるプログラム例
Function E1(fa, k) '第1種楕円積分,ガウス変換による
kr = k: ph = 1.5707963269
If fa < 0 Or fa > ph Then E1 ="faが不適切": Exit Function
If k < 0 Or k > 1 Then E1 = "kが不適切": Exit Function
s = Sin(fa)
If s = 1 And k > 0.9999999999 Then E1 ="kとfaが不適切": Exit Function
If k > 0.9999999999 Then E1 = Log((1 +Sin(fa)) / (1 - Sin(fa))) / 2: Exit Function
a = 1
Do
k2= k^2: If k2= 0 Then kt=0 else kt = (1 - Sqr(1 - k2))^ 2 / k2
If kt < 0.0000000001 Then Exit Do
k= kt: b = 1 + k: a = a * b: If s <> 0 Then s= (b- Sqr(b ^ 2- 4 * k * s ^ 2))/ 2/ k /s
Loop
If s < 1 Then E1= Atn(s / Sqr(1 - s ^ 2)) * a Else E1= ph * a
k = kr
End Function
ガウス変換を使った第1種楕円積分の計算プログラムを作る過程で副産物として得られた
ヤコビの楕円関数の計算プログラムをここに記載します。
ヤコビの楕円関数は第1種楕円積分z=∫x_0 dx/√((1-k^2 x^2)(1-x^2))の逆関数として定義される関数で
x= sn(z, k) と表現される。
ヤコビの振幅関数は第1種楕円積分z=∫φ_0 dφ/√(1-k^2 sin^2φ) の逆関数として定義される関数でφ= am(z, k)と表現される。
(x= sinφでありヤコビの楕円関数はヤコビの振幅関数から計算出来るのでヤコビの振幅関数のプログラムのみ記載。(^^
Function am(z, k) 'ヤコビの振幅関数
z= z* 1.0
If k < 0 Or k > 1 Then am = "kが不適切": Exit Function
If z < 0 Then am = "zが不適切": Exit Function
kr = k: ph = 1.5707963268: a = 1: le = 1
Do
k2= k^2: If k2= 0 Then kt=0 else kt = (1 - Sqr(1 - k2))^ 2 / k2
If kt < 0.0000000001 Then Exit Do
k = kt: a = a * (1 + k): le = le + 1
Loop
le = le - 1: zh = z/a
If zh > ph Then am = "zが不適切": k = kr: Exit Function
s = Sin(zh)
While le > 0
s = (1 + k) * s / (1 + k * s ^ 2): k = 2 * Sqr(k) / (1 + k): le = le - 1
Wend
If s < 1 Then am = Atn(s / Sqr(1 - s ^ 2)) Else am = ph
k = kr
End Function
Function E2(fa, k) '第2種楕円積分,級数展開による 検算用としてお使い下さい。
If fa < 0 Or fa > 1.5707963269 ThenE2 = "faが不適切": Exit Function
If k < 0 Or k > 1 Then E2 = "kが不適切": Exit Function
if k > 0.9999999999 then E2 = Sin(fa) :Exit Function
i = fa: ku = 2: c = Cos(fa): s = Sin(fa): f= 0: a = 1: k2 = k ^ 2
While Abs(i * a) > 0.0000000001
f = i * a + f: i = ((ku - 1) * i - s ^ (ku- 1) * c) / ku: a = a * (ku - 3) / ku * k2
ku = ku + 2
Wend
E2 = f
End Function
制御変数kuはローマ字でカウンターkauntarから取っているのです。
第2種楕円積分は級数展開による計算でも結構速くプログラムも短くてこのままでも良いのですが
ガウス変換を使った楕円積分の計算の解説に書いているように
第2種楕円積分はガウス変換により第2種楕円積分と第1種楕円積分と初等関数に展開出来るので
(21a展開と呼ぼう E2(Φ,K)= (2/(1+k))E2(φ,k)+(k-1)E1(φ,k)+(2k/(1+k)) Ea(φ,k))
これを使えば若干精度は落ちるが計算速度は早くプログラムは短いのでこれをおすすめ項目にします。(^^
Function E2(fa, k) '第2種楕円積分,21a展開による
ph = 1.5707963269
If fa < 0 Or fa > ph Then E2 = "faが不適切": Exit Function
If k < 0 Or k > 1 Then E2 = "kが不適切": Exit Function
If k > 0.9999999999 Then E2 = Sin(fa): Exit Function
kr = k
a2 = 1: a1 = 0: e = 0: s = Sin(fa)
Do
k2= k^2: If k2= 0 Then kt=0 else kt = (1 - Sqr(1 - k2))^ 2 / k2
If kt < 0.0000000001 Then Exit Do
k = kt: b = 1 + k: If s <> 0 Then s = (b - Sqr(b ^ 2 - 4 * k * s ^ 2)) / 2 / k / s
s2 = s ^ 2
a1 = a1 * b + a2 * (k - 1): ac = a2 * 2 * k / b
c = 1 - s2: If c < 0 Then c = 0 Else c = Sqr(c)
e = e + ac * s * c * Sqr(1 - k ^ 2 * s2) / (1 + k * s2): a2 = a2 * 2 / b
Loop
If s < 1 Then tn = s / Sqr(1 - s ^ 2) Else tn = 1E+20
E2 = Atn(tn) * (a2 + a1) + e
k = kr
End Function
第3種楕円積分 E3(φ,k,kc)=∫φ_0 1/((1+kc sin^2φ)(√(1-k^2 sin^2φ)) dφ
=∫φ_0 (1/(1-kd sin^2φ))(1/√(1-k^2 sin^2φ)) dφ
(符号処理を簡単にするため-kc=kdとおく)
abs(kd) <1 のとき
被積分関数の後ろの部分1/√(1-k^2sin^2φ)= 1+(1/2)k^2 sin^2φ+(1/2)(3/4)k^4 sin^4φ+…=a_0+a_2 sin^2φ+a_4 sin^4φ+…
を前の分数の部分1/(1+kc sin^2φ)の分母1-kd sin^2φで割って出来る
被積分関数=1+(1kd+(1/2)k^2) sin^2φ+((1kd+(1/2)k^2)kd+(1/2)(3/4) k^4) sin^4φ+ …を積分して
E3(φ,ｋc,ｋ)=φ+(1kd+(1/2)k^2)S_2(φ)+((1kd+(1/2)k^2)kd+(1/2)(3/4) k^4) S_4(φ)+ …
=φ+(1kd+(1/2)k^2) S_2(φ)+((1kd+(1/2)k^2)kd+(1/2)(3/4) k^4) S_4(φ)+ …= b_0φ+b_2 S_2(φ)+ b_4 S_4(φ)+ …
係数b_nは漸化式b_n=b_(n-2)kd+a_nで計算出来る
abs(kd) >1 のとき
上の形では後ろの項ほどkd が掛かって大きくなり発散してしまうのでkd sin^2φを先頭に立て-kdsin^2φ+1 で割って出来る
被積分関数=(a_0+a_2/kd+a_4/kd^2+…)/(-kd sin^2φ+1)-(a_2/kd+a_4/kd^2+a_6/kd^3+…)-(a_4/kd+a_6/kd^2+a_8/kd^3+…)sin^2φ-…
=b_(-2)/(-kd sin^2φ+1)+b_0+b_2 sin^2φ+b_4 sin^4φ+… (3)
係数b_nは
b_-2=a_0+a_2/kd+a_4/kd^2+… =1+(1/2)k^2/kd+(1/2)(3/4)k^4/kd^2+…=1/√(1-k^2/kd)
b_0=a_2/kd+a_4/kd^2+a_6/kd^3+…=-(b_(-2)-a0)= -b_(-2)+a0
b_2より後は漸化式 b_n=b_(n-2) kd+a_n で計算出来る。
(3)を積分する。積分I_(n-2)=∫φ_0 1/(1-kd sin^2φ) dφ は
kd>1のときsk=√(kd-1); I_(n-2)= arctan(sk tan(fa))/sk
kd<-1のときsk=√(1-kd); I_(n-2)= ln((sk tan(fa)+ 1)/(sk tan(fa)-1))/(2sk)であるから
E3(φ,k,kc)= b_(-2) I_(n-2)+b_0φ+b_2 S_2+b_4 S_2+…
級数展開による第3種楕円積分の計算プログラム例。検算用としてお使い下さい。
Function E3(fa,kc,k) '第3種楕円積分, 次数ごとにまとめた級数展開による
kd = -kc: k2 = k ^ 2: s = Sin(fa): s2 = s ^ 2: sc = s * Cos(fa): a = 1
If Abs(kd) > 1 Then GoTo 10
i = fa: b = 1: bi = b * i: f = bi: ku = 2
While Abs(bi) > 0.0000000001
a = a * (ku - 1) / ku * k2: b = b * kd + a:i = ((ku - 1) * i - sc) / ku: sc = sc * s2: bi = b * i: f = bi + f
ku = ku + 2
Wend
GoTo 20
'abs(kd)>1のとき
10 b = 1 / Sqr(1 - k2 / kd): If kd < -1 Then sk = Sqr(1 - kd): i = Atn(sk * Tan(fa)) / sk: GoTo 14
sk = Sqr(kd - 1): tn = Tan(fa): i = Log((sk * tn + 1) / (sk * tn - 1)) / 2 / sk
14 f = b * i
b = -b + 1: i = fa: bi = b * i: f = bi + f:ku = 2: bm = Abs(b)
15 If Abs(bi) < 0.0000000001 Then GoTo 20
a = a * (ku - 1) / ku * k2: b = b * kd + a
If Abs(b) > bm Then GoTo 20 'bが前より大きくなったので脱出
i = ((ku - 1) * i - sc) / ku: sc = sc * s2:bi = b * i: f = bi + f
ku = ku + 2: bm = Abs(b): GoTo 15
20 E3 = f
End Function
第3種楕円積分のガウス変換を使った計算プログラムをここに記載します。kが1に近いとき強力です。(^^
第3種楕円積分∫Φ_0 dΦ/((1+KC sin^2Φ)√(1-K^2 sin^2Φ)) は
sinΦの分数f=1/(1+KC sin^2Φ)と第1種楕円積分の部分ｇ=dΦ/√(1-K^2 sin^2Φ)の積で表される。
分数1/(1+KC sin^2Φ)をガウス変換すると定数1と新しいパラメタを持つ2つの分数が現れる。
f=1+(KC(1+k)^2/rt)(1/(1+((mo+rt)/2) sin^2φ)-1/(1+((mo-rt)/2) sin^2φ))
さらにガウス変換するとそれぞれの分数から定数1と新しいパラメタを持つ二つの分数が現れるが定数1は打ち消しあって消え分数の数が2倍になる。
ガウス変換を繰り返すとｇは実用上定数aとなるので最初の定数1を積分したものと各分数を積分して加算したもの
に定数aを掛ければ元の積分の答えになる。
Function E3(fa, k, kc) '第3種楕円積分,ガウス変換による
kr = k: kcr = kc: ph = 1.5707963269
If fa < 0 Or fa > ph Then E3 = "faが不適切": Exit Function
If k < 0 Or k > 1 Then E3 = "kが不適切": Exit Function
s = Sin(fa)
If s = 1 And k > 0.9999999999 Then E3 = "kとfaが不適切": Exit Function
If Abs(1 + kc * s ^ 2) < 0.0000000001 Then E3 = "kcとfaが不適切": Exit Function
If k > 0.9999999999 Then
If kc < 0 Then sk = Sqr(-kc): u = sk * s: us = Log((1 - u) / (1 + u)) / 2 Else sk = Sqr(kc): us = Atn(sk * s)
E3 = (Log((1 + s) / (1 - s)) / 2 + sk * us) / (1 + kc): Exit Function
Else: k2 = k ^ 2: If Abs(kc + k2) < 0.0000000001 Then s2 = s ^ 2: E3 = E2(fa, k) / (1 + kc) + kc * s * Sqr(1 - s2) / Sqr(1 + kc * s2): Exit Function
End If
Dim kh(10), kch(10), ckh(10), ach(10), cah(10), fuh(10)
kch(0) = kc: ckh(0) = 0: a = 1: ach(0) = 1: cah(0) = 0: le = 1
Do
k2= k^2: If k2= 0 Then kt=0 else kt = (1 - Sqr(1 - k2))^ 2 / k2
If kt < 0.0000000001 Then Exit Do
k = kt: b = 1 + k: a = a * b: kh(le) = k: fuh(le) = 0:If s <> 0 Then s = (b - Sqr(b ^ 2 - 4 * k * s ^ 2)) / 2 / k / s
le = le + 1
Loop
If s < 1 Then tn =s / Sqr(1 - s ^ 2) Else tn = 1E+20
If kc = 0 Then E3 = Atn(tn) * a: k = kr: Exit Function
ll = le - 1: iw = Atn(tn): le = 1
Do: Do
If fuh(le) = -1 Then Exit Do
If fuh(le) = 1 Then fuh(le) = -1: If le < ll Then For j = le + 1 To ll: fuh(j) = 0: Next
If fuh(le) = 0 Then fuh(le) = 1
k = kh(le): k2 = k ^ 2: b = 1 + k: b2 = b ^ 2: kc = kch(le - 1): ck = ckh(le - 1): ac = ach(le - 1): ca = cah(le - 1)
mo = b2 * kc + 2 * k: om = b2 * ck: ru = mo ^ 2 - om ^ 2 - 4 * k2
ur = 2 * om * mo: w = Sqr(ru ^ 2 + ur ^ 2): rt = Sqr((ru + w) / 2): tr = Sqr((-ru + w) / 2): If ur < 0 Then tr = -tr
jo = ac * kc - ca * ck: oj = ac * ck + ca * kc: w2 = (rt ^ 2 + tr ^ 2) / b2
ac = (jo * rt + oj * tr) / w2: ca = (-jo * tr + oj * rt) / w2
fu = fuh(le): ac = fu * ac: ca = fu * ca: kc = (mo + fu * rt) / 2: ck = (om + fu * tr) / 2
ach(le) = ac: cah(le) = ca: kch(le) = kc: ckh(le) = ck
If le < ll Then
le = le + 1
Else
ru = 1 + kc: ur = ck: w = Sqr(ru ^ 2 + ur ^ 2): rt = Sqr((ru + w) / 2): tr = Sqr((-ru + w) / 2): If ur < 0 Then tr = -tr
au = rt * tn: ua = tr * tn: w2 = au ^ 2 + ua ^ 2: at = Atn(2 * au / (1 - w2)) / 2: If (1 - w2) < 0 Then at = Sgn(au) * ph + at
u = 2 * ua / (1 + w2): ta = Log((1 + u) / (1 - u)) / 4
jo = ac * at - ca * ta: oj = ac * ta + ca * at: w4 = rt ^ 2 + tr ^ 2
se = (jo * rt + oj * tr) / w4: iw = se + iw: End If
Loop
Do:le = le - 1:If fuh(le) > -1 Then Exit Do
Loop
If le = 0 Then Exit Do
Loop
E3 = iw * a: k = kr: kc = kcr
End Function
完全楕円積分
(主要文字をZとしているのはローマ字のkanzenのZを取っているのです。(^^はじめはKZとしようかとも思ったのですが母数のkと紛れ易くなるのでやめた)
(∫π/2_0は0からπ/2までの定積分を表しています。)
第1種完全楕円積分 Z1(k)=∫π/2_0 1/√(1-k^2 sin^2φ) dφ=π/2(1+(1/2)^2 k^2+((1/2)(3/4))^2 k^4+ …)
第2種完全楕円積分 Z2(k)=∫π/2_0 √(1-k^2 sin^2φ) dφ=π/2(1+(-1/2)(1/2) k^2+(-1/2)(1/4)(1/2)(3/4) k^4+ …
第3種完全楕円積分 Z3(k,kc)=∫π/2_0 1/((1+kc sin^2φ) √(1-k^2 sin^2φ)) dφ
VBによるプログラム例
第1種完全楕円積分は前出のガウス変換をkが実用上0と見なせるまで数回繰り返して行うアルゴリズムが簡単で早くおすすめです。
この簡単明瞭さを見よ。
Function Z1(k) '第1種完全楕円積分,ガウス変換による
If k < 0 Or k > 0.9999999999 Then Z1 = "kが不適切": Exit Function
kr = k: ph = 1.5707963268: a = 1
Do
k2= k^2: If k2= 0 Then kt=0 else kt = (1 - Sqr(1 - k2))^ 2 / k2
If kt < 0.0000000001 Then Exit Do
k = kt: a = a *(1 + k)
Loop
Z1 = ph * a: k = kr
End Function
Function Z2(k) '第2種完全楕円積分,級数展開による
If k < 0 Or k > 1 Then Z2 = "kが不適切": Exit Function
if k =1 then Z2 = 1: Exit Function
i = 1.5707963268: ku = 2: f = 0: a = 1: k2 = k ^ 2
While Abs(i * a) > 0.0000000001
f = i * a + f: i = (ku - 1) / ku * i: a = a * (ku - 3) / ku * k2
ku = ku + 2
Wend
Z2 = f
End Function
Function Z2(k) '第2種完全楕円積分,21a展開による
If k < 0 Or k > 1 Then Z2 = "kが不適切": Exit Function
If k = 1 Then Z2 = 1: Exit Function
kr = k
a2 = 1: a1 = 0
Do
k2= k^2: If k2= 0 Then kt=0 else kt = (1 - Sqr(1 - k2))^ 2 / k2
If kt < 0.0000000001 Then Exit Do
k = kt: b = 1 + k: a1 = a1 * (1 + k) + a2 * (k - 1): a2 = a2 * 2 / b
Loop
Z2 = 1.5707963268 * (a2 + a1)
k = kr
End Function
Function Z3(k, kc) '第3種完全楕円積分,次数ごとにまとめた級数展開による,検算用としてお使い下さい
kd = -kc: k2 = k ^ 2: ph = 1.5707963268: a = 1
If Abs(kd) > 1 Then GoTo 10
i = ph: b = 1: bi = b * i: f = bi: ku = 2
While Abs(bi) > 0.0000000001
a = a * (ku - 1) / ku * k2: b = b * kd + a: i = (ku - 1) / ku * i: bi = b * i: f = bi + f
ku = ku + 2
Wend
GoTo 20
'abs(kd)>1のとき
10 b = 1 / Sqr(1 - k2 / kd): If kd < -1 Then sk = Sqr(1 - kd): i = ph / sk: GoTo 14
i = 0
14 f = b * i
b = -b + 1: i = ph: bi = b * i: f = bi + f: ku = 2: bm = Abs(b)
15 If Abs(bi) < 0.0000000001 Then GoTo 20
a = a * (ku - 1) / ku * k2: b = b * kd + a
If Abs(b) > bm Then GoTo 20 'bが前より大きくなったので脱出
i = (ku - 1) / ku * i: bi = b * i: f = bi + f
ku = ku + 2: bm = Abs(b): GoTo 15
20 Z3 = f
End Function
Function Z3(k, kc) '第3種完全楕円積分,ガウス変換による
kr = k: kcr = kc
If k < 0 Or k > 0.9999999999 Then Z3 = "kが不適切": Exit Function
If Abs(1 + kc) < 0.0000000001 Then Z3 = "kcが不適切": Exit Function
If Abs(kc + k ^ 2) < 0.0000000001 Then Z3 = Z2(k) / (1 + kc): Exit Function
Dim kh(10), kch(10), ckh(10), ach(10), cah(10), fuh(10)
ph = 1.5707963268: kch(0) = kc: ckh(0) = 0: a = 1: ach(0) = 1: cah(0) = 0: le = 1
Do
k2= k^2: If k2= 0 Then kt=0 else kt = (1 - Sqr(1 - k2))^ 2 / k2
If kt < 0.0000000001 Then Exit Do
k = kt: b = 1 + k: a = a * b: kh(le) = k: fuh(le) = 0: le = le + 1
Loop
If kc = 0 Then Z3 = ph * a: k = kr: Exit Function
ll = le - 1: iw = ph: le = 1
Do: Do
If fuh(le) = -1 Then Exit Do
If fuh(le) = 1 Then fuh(le) = -1: If le < ll Then For j = le + 1 To ll: fuh(j) = 0: Next
If fuh(le) = 0 Then fuh(le) = 1
k = kh(le): k2 = k ^ 2: b = 1 + k: b2 = b ^ 2: kc = kch(le - 1): ck = ckh(le - 1): ac = ach(le - 1): ca = cah(le - 1)
mo = b2 * kc + 2 * k: om = b2 * ck: ru = mo ^ 2 - om ^ 2 - 4 * k2
ur = 2 * om * mo: w = Sqr(ru ^ 2 + ur ^ 2): rt = Sqr((ru + w) / 2): tr = Sqr((-ru + w) / 2): If ur < 0 Then tr = -tr
jo = ac * kc - ca * ck: oj = ac * ck + ca * kc: w2 = (rt ^ 2 + tr ^ 2) / b2
ac = (jo * rt + oj * tr) / w2: ca = (-jo * tr + oj * rt) / w2
fu = fuh(le): ac = fu * ac: ca = fu * ca: kc = (mo + fu * rt) / 2: ck = (om + fu * tr) / 2
ach(le) = ac: cah(le) = ca: kch(le) = kc: ckh(le) = ck
If le < ll Then
le = le + 1
Else ru = 1 + kc: ur = ck: w = Sqr(ru ^ 2 + ur ^ 2): rt = Sqr((ru + w) / 2): tr = Sqr((-ru + w) / 2): If ur < 0 Then tr = -tr
jo = ac * ph: oj = ca * ph: w4 = rt ^ 2 + tr ^ 2: se = (jo * rt + oj * tr) / w4: iw = se + iw: End If
Loop
Do:le = le - 1:If fuh(le) > -1 Then Exit Do
Loop
If le = 0 Then Exit Do
Loop
Z3 = iw * a: k = kr: kc = kcr
End Function
第3種楕円積分の有理関数部が2乗であるもの(J2積分)の計算 J2(φ,k,kc)=∫φ_0 (1/((1+kc sin^2φ)^2 √(1-k^2sin^2φ))) dφ
初等関数Ea3= sinφcosφ√(1-k^2 sin^2φ)/(1+kc sin^2φ)を微分すると
d/dφ(Ea3)= ((k^2/kc) sin^6φ+3k^2 sin^4φ-(2+2k^2+kc) sin^2φ+1)/((1+kc sin^2φ)^2 √(1-k^2 sin^2φ))=A となる。
これをα,β,γ,δを未定の係数として第2種楕円積分、第1種楕円積分、第3種楕円積分、J2積分の被積分関数の和の形に展開してBにしたい。
B= α√(1-k^2 sin^2φ)+β/√(1-k^2 sin^2φ)+γ/((1+kc sin^2φ) √(1-k^2 sin^2φ))+δ/((1+kc sin^2φ)^2 √(1-k^2sin^2φ))
Bを通分すると
B=((-k^2 kc^2 α) sin^6φ+((-2 k^2 kc+kc^2)α+kc^2 β) sin^4φ+((-k^2+2kc)α+2kc β+kc γ) sin^2φ+(α+β+γ+δ))/((1+kc sin^2φ)^2 √(1-k^2 sin^2φ))
となる。AとBが等しいとき分子のsinφの6、4、2、0次の係数が等しいので次の四つの式が成り立つ。
-k^2 kc^2 α= k^2 kc；(-2 k^2 kc+kc^2)α+kc^2 β= 3k^2；(-k^2+2kc)α+2kc β+kc γ= -(2+2k^2+kc)；α+β+γ+δ= 1
これから係数がα= -1/kc；β=k^2/kc^2+1/kc；γ=-(1+2/kc+2k^2/kc+3k^2/kc^2)；δ=2(1+1/kc+k^2/kc+k^2/kc^2)=2(k^2+kc)(1/kc+1/kc^2) と求められ
d/dφ(Ea3)= (-1/kc)√(1-k^2 sin^2φ)+(k^2/kc^2+1/kc)/√(1-k^2 sin^2φ)-(1+2/kc+2k^2/kc+3k^2/kc^2)/((1+kc sin^2φ) √(1-k^2 sin^2φ))
+2(k^2+kc)(1/kc+1/kc^2)/((1+kc sin^2φ)^2 √(1-k^2sin^2φ)) となる。
両辺を0からφまで積分すると
sinφcosφ√(1-k^2 sin^2φ)/(1+kc sin^2φ)=-(1/kc) E2(φ,k)+(k^2/kc^2+1/kc) E1(φ,k)-(1+2/kc+2k^2/kc+3k^2/kc^2) E3(φ,k,kc)
+2(k^2+kc)(1/kc+1/kc^2) J2(φ,k,kc)
となる。これを方程式とみなしてJ2(φ,kc,k)について解くと
J2(φ,k,kc)= (1/(2(k^2+kc)(1/kc+1/kc^2))(sinφcosφ√(1-k^2 sin^2φ)/(1+kc sin^2φ)+(1/kc) E2(φ,k)-(k^2/kc^2+1/kc) E1(φ,k)+(1+2/kc+2k^2/kc+3k^2/kc^2) E3(φ,kc,kc)))
となる。
J2積分でφ=π/2のとき(J2完全積分)JZ2と表記すると
JZ2(k,kc)=∫π/2_0 1/((1+kc sin^2φ)^2 √(1-k^2 sin^2φ)) dφ
=(1/(2(k^2+kc)(1/kc+1/kc^2))((1/kc) Z2(k)-(k^2/kc^2+1/kc) Z1(k)+(1+2/kc+2k^2/kc+3k^2/kc^2) Z3(k,kc))
楕円積分にcosφを掛けた形の積分可能な関数の積分
第1種形 E1a=∫φ_0 cosφ/√(1-k^2 sin^2φ) dφ= (1/k) arcsin(k sinφ)
第2種形 E2a=∫φ_0 cosφ √(1-k^2 sin^2φ) dφ= (1/2)(sinφ √(1-k^2 sin^2φ)+(1/k)arcsin(k sinφ))
第3種形 E3a=∫φ_0 cosφ/((1+kc sin^2φ)√(1-k^2 sin^2φ)) dφ= (1/√(k^2+kc)) arctan(√(k^2+kc) sinφ/√(1-k^2 sin^2φ))
J2積分形 J2a=∫φ_0 cosφ/((1+kc sin^2φ)^2 √(1-k^2 sin^2φ)) dφ
= (1/(2(k^2/kc+1))(sinφ √(1-k^2 sin^2φ)+((2k^2/kc+1)/√(k^2+kc)) arctan(√(k^2+kc) sinφ/√(1-k^2 sin^2φ))
ベッセル関数の微分 J`_n(x)= -n J_n(x)/x+J_n-1(x) (J`_n(x)はd/dx (J_n(x)) を表しています)
変形されたベッセル関数の微分 I`_n(x)= -n I_n(x)/x+I_n-1(x)
完全楕円積分の微分
φで積分したものをkで微分したものとkで微分したものをφで積分したものは等しいので先にkで微分したものを
積分する事により次の公式が得られる
d/dk (Z1(k))=(1/k)(-Z1(k)+(1/(1-k^2))Z2(k))
d/dk (Z2(k))=(1/k)(-Z1(k)+Z2(k))
第3種完全楕円積分の母数kによる微分
∂/∂k (Z3(k,kc)) = ∫π/2_0 dφ k sin^2φ/((1+kc sin^2φ)(1-k^2 sin^2φ) √(1-k^2 sin^2φ))
= ∫π/2_0 dφ(k/(kc+k^2))(-1/(1+kc sin^2φ)+1/(1-k^2 sin^2φ))/√(1-k^2 sin^2φ))
= (k/(kc+k^2))((1/(1-k^2))Z2(k)-Z3(k,kc))
第3種完全楕円積分のパラメタkcによる微分
∂/∂kc (Z3(k,kc)) =∫π/2_0 dφ(-sin^2φ/((1+kc sin^2φ)^2 √(1-k^2 sin^2φ)))
=∫π/2_0 dφ(1/kc)(-1/(1+kc sin^2φ)+1/(1+kc sin^2φ)^2)/√(1-k^2 sin^2φ)
=(1/kc)(-(1/kc)(1+k^2/kc) Z1(k)+(1/kc) Z2(k)+(-1+k^2/kc^2) Z3(k,kc))/(2(1+(1+k^2)/kc+k^2/kc^2))
積分区間0〜1の定積分
三角関数、双曲線関数、冪 の積
[f`g]1_0=∫1_0 f``g dx+∫1_0 f`g` dx ; [fg`]1_0 =∫1_0 f`g` dx+∫1_0 f g`` dxの関係を利用して計算すると良い
∫1_0 cos lx dx= sin l /l
∫1_0 cosh mx dx=sinh m/m
∫1_0 cos^2 lx dx=(1+(sin l cos l)/l)/2
∫1_0 sin^2 lx dx=(1-(sin l cos l)/l)/2
∫1_0 cosh^2 mx dx=(1+(sinh m cosh m)/m)/2
∫1_0 sinh^2 mx dx=(-1+(sinh m cosh m)/m)/2
∫1_0 cos lx cos mx dx=(l sin l cos m-m cos l sin m)/D (l^2-m^2=Dと略記)
∫1_0 sin lx sin mx dx = (m sin l cos m-l cos l sin m)/D
∫1_0 cos lx cosh mx dx=(l sin l cosh m+m cos l sinh m)/W (l^2+m^2=Wと略記)
∫1_0 sin lx sinh mx dx = (m sin l cosh m-l cos l sinh m)/W
∫1_0 cosh lx cosh mx dx=(l sinh l cosh m-m cosh l sinh m)/D
∫1_0 sinh lx sinh mx dx = (-m sinh l cosh m+l cosh l sinh m)/D
∫1_0 x sin lx dx=sin l/l^2-cos l/l
∫1_0 x sinh mx dx=-sinh m/m^2+cosh m/m
∫1_0 x sin lx cos lx dx =sin l cos l/(4l^2)+(sin^2 l-cos^2 l)/(4l)
∫1_0 x sinh mx cosh mx dx =-sinh m cosh m/(4m^2)+(sinh^2 m+cosh^2 m)/(4m)
∫1_0 x sin lx cos mx dx =-(l cos l cos m+m sin l sin m)/D+(W sin l cos m-2lm cos l sin m)/D^2
∫1_0 x sin lx cosh mx dx =(-l cos l cosh m+m sin l sinh m)/W +(D sin l cosh m+2lm cos l sinh m)/W^2
∫1_0 x cos lx sinh mx dx =(m cos l cosh m+l sin l sinh m)/W+(D cos l sinh m+2lm sin l cosh m)/W^2
∫1_0 x sinh lx cosh mx dx =(l cosh l cosh m-m sinh l sinh m)/D+(-W sinh l cosh m+2lm cosh l sinh m)/D^2
∫1_0 x^2 cos lx dx= sin l (1/l-2/l^3)
∫1_0 x^2 cosh mx dx=sinh m (1/m-2/m^3)
∫1_0 x^2 cos^2 lx dx=1/6+sin l cos l (1/(2l)-1/(4l^3))-(sin^2 l-cos^2 l)/(4l^2)
∫1_0 x^2 sin^2 lx dx=1/6-sin l cos l (1/(2l)-1/(4l^3))+(sin^2 l-cos^2 l)/(4l^2)
∫1_0 x^2 cosh^2 mx dx=1/6+sinh m cosh m (1/(2m)-1/(4m^3))-(sinh^2 m+cosh^2 m)/(4m^2)
∫1_0 x^2 sinh^2 mx dx=-1/6+sinh m cosh m (1/(2m)-1/(4m^3))-(sinh^2 m+cosh^2 m)/(4m^2)
∫1_0 x^2 cos lx cos mx dx=(l sin l cos m-m cos l sin m)/D+2(W cos l cos m+2l m sin l sin m)/D^2-2(l(W+2m^2)sin l cos m-m(W+2l^2)cos l sin m)/D^3
∫1_0 x^2 sin lx sin mx dx = (m sin l cos m-l cos l sin m)/D+2(2l m cos l cos m+ W sin l sin m)/D^2-2(m(W+2l^2) sin l cos m-l(W+2m^2) cos l sin m)/D^3
∫1_0 x^2 cos lx cosh mx dx=(l sin l coshm+m cos l sinh m)/W+2(D cos l cosh m-2l m sin l sinh m)/W^2-2(l(D+2m^2)sin l cosh m-m(D+2l^2)cos l sinh m)/W^3
∫1_0 x^2 sin lx sinh mx dx = (m sin l cosh m-l cos l sinh m)/W+2(2l m cos l cosh m-D sin l sinh m)/W^2+2(-m(D+2l^2) sin l cosh m- l(D+2m^2)cos l sinh m)/W^3
∫1_0 x^2 cosh lx cosh mx dx=(l sinh l cosh m-m cosh l sinh m)/D+2(-W cosh l cosh m+2l m sinh l sinh m)/D^2+2(l(W+2m^2)sinh l cosh m-m(W+2l^2)cosh l sinh m)/D^3
∫1_0 x^2 sinh lx sinh mx dx = (-m sinh l cosh m+l cosh l sinh m)/D+2(2l m cosh l cosh m+W sinh l sinh m)/D^2+2(-m(W+2l^2)sinh l cosh m+l(W+2m^2) cosh l sinh m)/D^3
ベッセル関数、冪 の積
[f`g x]1_0 =∫1_0(f``+f`/x)g x dx +∫1_0 f`g` x dx ; [f g` x]1_0 =∫1_0 f`g` x dx+∫1_0 f(g``+g`/x) x dxの関係を利用して計算すると良い