楕円弧の長さは級数解や台形公式による数値積分によっても計算できますが 処理速度を考えると ガウス－ルジャンドルあるいは春日屋の方法が適しています。
楕円弧の長さは
で計算できます。
楕円弧の長さは 楕円の長径×第２種楕円積分 で与えられることがわかります。第２種楕円積分は楕円弧の円周率πのようなものだと考えれば受け入れやすいと思います。
wikipedia 楕円積分によれば、ルジャンドル標準形の楕円積分になります。数式処理ソフト maxima の組込関数はこの形式が採用されています。
※参考文献
○楕円弧の長さの計算式は、応用数学 第４巻 建築構造講座 コロナ社
○Gaussの求積法は、解析概論[改訂第三版] 岩波書店
○Legendreの多項式は、自然科学者のための数学概論[増訂版] 岩波書店
○平均値法は、次元解析・最小２乗法と実験式 応用数学講座 第５巻 コロナ社 本間 仁・春日屋伸昌 共著
今回は 楕円弧の長さの計算式を 級数解で解く やり方を紹介します。
つぎのプログラムを 200-3-1.bat で保存して実行すると
:楕円弧の長さ(級数解)
@echo off
path C:\maxima-5.47.0\bin;%path%
set RUBYLIB=\jww\Lite\pro\ruby\lib
ruby -x %~f0
pause
goto:eof
#!ruby -Ku -rjw
f = open("| maxima --very-quiet", "w")
f.print <<~Maxima
u : taylor(sqrt(1 - t), t, 0, 100)$
e : ratsubst(m * sin(p)^2, t, u)$
[a, b] : [1, 0.5]$
m : 1 - (b / a)^2$
print("[a, b, m] = ", [a, b, m])$
s : integrate(e, p, 0, %pi/2)$
print("s =", float(ev(s)))$
[p1, p2] : [30, 120]$
p1 : (90 - p1) * %pi / 180$
p2 : (p2 - 90) * %pi / 180$
s1 : integrate(e, p, 0, p1)$
s2 : integrate(e, p, 0, p2)$
print("s =", float(ev(s1+s2)))$
quit()$
Maxima
f.close
# JW_CAD 外部変形ライブラリを利用して計算(検算用)
#$n = 50
#puts "s = %s" % "0 0 1 0 90 0.5 0".ci.lr
#puts "s = %s" % "0 0 1 30 120 0.5 0".ci.lr
__END__
計算結果は
となります。級数の項数を 100 としているせいか、処理時間がかかります。
参考までに つぎの ruby プログラムを 200-3-2.bat で保存して実行すると
:楕円弧の長さ(級数解 参考資料)
@echo off
ruby -x %~f0
pause
goto:eof
#!ruby
#ここから
#ここまで
$n = 100 #級数の項数
a, b, p1, p2 = 1, 0.5, 0, 90
puts "a b p1 p2 = %s %s %s %s" % [a, b, p1, p2]
puts "s = %s" % s(a, p1, p2, b/a)
a, b, p1, p2 = 1, 0.5, 30, 120
puts "a b p1 p2 = %s %s %s %s" % [a, b, p1, p2]
puts "s = %s" % s(a, p1, p2, b/a)
__END__
計算結果は
となります。処理速度の問題は改善されています。
有効桁数15 程度の演算で 級数解は、極端に平べったい楕円で誤差がやや大きくなります。
ネット上で完全楕円積分の級数解は見つけだすことできるのですが、不完全楕円積分の級数解を目にすることはあまりありません。
[ 級数展開の補足 ]
楕円弧の長さの計算式を maxima で素直に級数展開するプログラムを書くと以下のようになります。
:楕円弧の長さ(級数展開-2)
@echo off
path C:\maxima-5.47.0\bin;%path%
ruby -x %~f0
pause
goto:eof
#!ruby
f = open("| maxima --very-quiet", "w")
f.print <<~Maxima
e : taylor(sqrt(1 - m * sin(p)^2), p, 0, 8)$
print("sqrt(1 - m * sin(p)^2) = ", expand(e))$
[a, b] : [1, 0.5]$
m : 1 - (b / a)^2$
print("[a, b, m] = ", [a, b, m])$
s : integrate(e, p, 0, %pi/2)$
print("s =", float(ev(s)))$
[p1, p2] : [30, 120]$
p1 : (90 - p1) * %pi / 180$
p2 : (p2 - 90) * %pi / 180$
print("[p1, p2] = ", [p1, p2])$
s1 : integrate(e, p, 0, p1)$
s2 : integrate(e, p, 0, p2)$
print("s =", float(ev(s1+s2)))$
quit()$
Maxima
f.close
__END__
計算結果は
となります。それらしくありませんが、問題ないようです。数式は、答えがひとつではないので迷います。思い込みや先入観は守る側から見れば貴重ですが発想の芽を摘んでいるのも事実です。・・・明治維新から続く親方日の丸の意識はいつまで続くのでしょうか。欧米に追い付け追い越せで、本当に追いついているのか心配です・・・
このプログラムの級数解は、以下のページを参考にしました。
◯楕円積分 数値計算ノート(現在は接続できません)
http://www.ab.cyberhome.ne.jp/~naka-kevin/Elliptic-Integral.html
#ここから
#ここまで
を下記に差し替えてください
# p=φ の補正°
def apa(o, p)
return (90.0 * (2 * ( (o + 1) / 2).to_i - 1) - p) * (-1) ** (o - 1)
end
# 級数解による第２種楕円積分の計算
def E2(p, k = nil)
# E2(k) 第２種完全楕円積分
# E2(p, k) 第２種不完全楕円積分
if k == nil
k = p
p = PI / 2
end
if k < 0 or 1 < k then return print "h#E2 p, %s 計算できませんでした\n" % k end
if k > 1 - 1e-15 then return Math.sin(p) end
n, e, a, s_n = 2.0, 0.0, 1.0, p
c, s = Math.cos(p), Math.sin(p)
j = 0
while (a * s_n).abs > 1e-20
e += a * s_n
a *= (n - 3) / n * k * k
s_n = ( (n - 1) * s_n - s ** (n - 1) * c) / n
n += 2
j += 1
if j > $n then break end
end
return e
end
# 楕円弧の長さの計算
def s(r, p1, p2, w)
# s(r, p1, p2, w)
pi = Math::PI
if w > 1
r *= w
w = 1.0 / w
p1 += 90.0
p2 += 90.0
end
ek = E2(pi / 2, k = Math.sqrt(1 - w * w))
p2 += 360 if p1 == p2
lp = ( (p2 - p1).abs / 360.0).to_i * 4 * ek
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 += ek * ( (n1 + 1) % 2) + E2(apa(n1, p1) * pi / 180, k) * (-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 += ek * ( (n3 - n1) % 4) + E2(apa(n2, p2) * pi / 180, k) * (-1.0) ** n2
lp += 4 * ek if lp < 0.0
return r * lp
end
end
