楕円弧の長さ(5) －台形積分－

楕円弧の長さは級数解や台形公式による数値積分によっても計算できますが処理速度を考えるとガウス－ルジャンドルあるいは春日屋の方法が適しています。

楕円弧の長さは

で計算できます。

楕円弧の長さは楕円の長径×第２種楕円積分で与えられることがわかります。第２種楕円積分は楕円弧の円周率πのようなものだと考えれば受け入れやすいと思います。
wikipedia 楕円積分によれば、ルジャンドル標準形の楕円積分になります。数式処理ソフト maxima の組込関数はこの形式が採用されています。

※参考文献
○楕円弧の長さの計算式は、応用数学第４巻建築構造講座コロナ社
○Gaussの求積法は、解析概論[改訂第三版] 岩波書店
○Legendreの多項式は、自然科学者のための数学概論[増訂版] 岩波書店
○平均値法は、次元解析・最小２乗法と実験式応用数学講座第５巻コロナ社本間　仁・春日屋伸昌共著

今回は台形積分を紹介します。

台形積分の詳細は wikipedia 台形公式が参考になります。

つぎのプログラムを 200-5-1.bat で保存して実行すると

:楕円弧の長さ(台形積分)
@echo off
ruby -x %~f0
pause
goto:eof

#!ruby
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, a, b, n=$n)
n = 100 if $n == nil
hh = (b - a) / n.to_f
x = 0.0 ; l = 0.0
for i in 1..n
s1 = f(x * p, k)
x += hh
s2 = f(x * p, k)
l += (s1 + s2) * hh / 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)
ec = E2(pi / 2.0, k, 0, 1)
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, 0, 1) * (-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, 0, 1) * (-1.0) ** n2
lp += 4 * ec if lp < 0.0
return r * lp
end
end

puts "trapezoidal rule"
$n = 10000 #積分区間の分割数
puts "division 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