楕円弧の長さは級数解や台形公式による数値積分によっても計算できますが 処理速度を考えると ガウス－ルジャンドルあるいは春日屋の方法が適しています。
楕円弧の長さは
で計算できます。
楕円弧の長さは 楕円の長径×第２種楕円積分 で与えられることがわかります。第２種楕円積分は楕円弧の円周率πのようなものだと考えれば受け入れやすいと思います。
wikipedia 楕円積分によれば、ルジャンドル標準形の楕円積分になります。数式処理ソフト maxima の組込関数はこの形式が採用されています。
※参考文献
○楕円弧の長さの計算式は、応用数学 第４巻 建築構造講座 コロナ社
○Gaussの求積法は、解析概論[改訂第三版] 岩波書店
○Legendreの多項式は、自然科学者のための数学概論[増訂版] 岩波書店
○平均値法は、次元解析・最小２乗法と実験式 応用数学講座 第５巻 コロナ社 本間 仁・春日屋伸昌 共著
ガウス－ルジャンドルの方法の説明は省略します。参考文献を参照ください。
つぎのプログラムを 200-7-1.bat で保存して実行すると
:楕円弧の長さ(ガウス－ルジャンドルの方法)
@echo off
ruby -x %~f0
pause
goto:eof
#!ruby
def gauss(n)
pi = Math::PI
ee = Array.new(1, 0.0)
hh = Array.new(1, 0.0)
m = (n + n % 2) / 2
e1 = n * (n + 1.0)
# ここから Fortran Program をルビーに移植した
for i in 1..m do
t = (4.0 * i - 1.0) * pi / (4.0 * n + 2.0)
x0 = (1.0 - (1.0 - 1.0 / n) / (8.0 * n * n)) * Math.cos(t)
pkm1= 1.0
pk = x0
for k in 2..n do
t1 = x0 * pk
pkp1 = t1 - pkm1 - (t1 - pkm1) / k + t1
pkm1 = pk
pk = pkp1
end
den = 1.0 - x0 * x0
d1 = n * (pkm1 - x0 * pk)
dpn = d1 / den
d2pn = (2.0 * x0 * dpn - e1 * pk) / den
d3pn = (4.0 * x0 * d2pn + (2.0 - e1) * dpn) / den
d4pn = (6.0 * x0 * d3pn + (6.0 - e1) * d2pn) / den
u = pk / dpn
v = d2pn / dpn
h = -1.0 * u * (1.0 + 0.5 * u * (v + u * (v * v - d3pn / (3.0 * dpn))))
pp = pk + h * (dpn + 0.5 * h * (d2pn + h / 3.0 * (d3pn + 0.25 * h * d4pn)))
dp = dpn + h * (d2pn + 0.5 * h * (d3pn + h * d4pn / 3.0))
h -= pp / dp
xi = x0 + h
fx = d1 - h * e1 * (pk + 0.5 * h * (dpn + h / 3.0 * (d2pn + 0.25 * h * (d3pn + 0.2 * h * d4pn))))
wi = 2.0 * (1.0 - xi * xi) / (fx * fx)
ee[i] = xi
hh[i] = wi
end
return ee, hh
end
def f(x, k)
return Math.sqrt(1.0 - (k * Math.sin(x)) ** 2)
end
def apa(o, p)
return (90.0 * (2 * ( (o + 1) / 2).to_i - 1) - p) * (-1) ** (o - 1)
end
def E2(p, k, ee, hh)
n = ($n == nil ? 60 : $n)
m = (n + n % 2) / 2
if k == 0 || n < 1
return p
end
l = 0.0
for i in 1..n / 2
l += hh[i] * f(ee[i] * p, k)
end
if m == n / 2 + 1 #奇数個の場合
l += hh[m] * f(ee[m] * p, k) / 2.0
end
return l * p
end
def s(a, b, p1, p2)
pi = Math::PI
w = b / a
r = a
if w > 1
w = 1.0 / w
r = r / w
p1 = p1 - 90.0
p2 = p2 - 90.0
end
k = Math.sqrt(1.0 - w * w)
n = ($n == nil ? 60 : $n)
ee, hh = gauss(n)
ec = E2(pi / 2.0, k, ee, hh)
lp = ( (p2 - p1).abs / 360.0).to_i * 4 * ec
p1 += 360 while p1 < 0
p1 -= 360 while p1 >= 360
p2 += 360 while p2 < 0
p2 -= 360 while p2 >= 360
if p1 == p2
return r * lp
else
n1 = 1 + (p1 / 90.0).to_i
n1 = 4 if n1 > 4
lp += ec * ( (n1 + 1) % 2) + E2(apa(n1, p1) * pi / 180, k, ee, hh) * (-1.0) ** (n1 - 1)
n2 = 1 + (p2 / 90.0).to_i
n2 = 4 if n2 > 4
n1 += 4 if n1 == n2
n3 = 5 + 2 * ( (n2 - 1) / 2.0).to_i
lp += ec * ( (n3 - n1) % 4) + E2(apa(n2, p2) * pi / 180, k, ee, hh) * (-1.0) ** n2
lp += 4 * ec if lp < 0.0
return r * lp
end
end
puts "Gauss-Legendre's method"
$n = 30 #積分点数
puts "node n = %s" % $n
puts
a, b, p1, p2 = 1.0, 0.5, 0, 90
puts "a b p1 p2 = %s %s %s %s" % [a, b, p1, p2]
puts "s = %s" % s(a, b, p1, p2)
puts "true value"
puts "s = %s" % 1.2110560275684596
puts
a, b, p1, p2 = 1.0, 0.5, -30, 120
puts "a b p1 p2 = %s %s %s %s" % [a, b, p1, p2]
puts "s1 + s2 = %s" % s(a, b, p1, p2)
puts "true value"
puts "s1 + s2 = %s" % 2.00981083328632
計算結果は
となります。
ガウスールジャンドルの方法は楕円弧の計算に適しています。問題があるとすれば、積分点の位置(プログラムでee)と重み(プログラムでhh)があらかじめ必要なことです。探せば 20点 くらいまでは見つかるのですが、それ以上は計算しなければならないのです。
この方法は有限要素法でよく使われているせいか、当方は、幸いにして、Fortran Program にあるのを発見しました。それまでは 積分点を 150 に決め打ちしていたので円に近いものは不経済となり直線に近いものは誤差が出るといった始末でした。積分点の位置と重みを手軽に計算できるサブルーチンを開発した方には深くお礼申し上げます。
次回は 計算結果のまとめ を取り上げます。
