楕円は円錐を平面で切ったときにできる図形のうちの1つです。今回は楕円の周の長さCを数値計算によって求めます。以下の式で表現される、長軸2a、単軸2bの平面上の楕円を考えます。
\frac{x^2}{a^2}+\frac{y^2}{b^2}=1楕円の面積はπabと簡単に計算できますが、楕円の周の長さを普通の関数で表すことはできません。ゆえに数値計算が必要になります。周の長さは以下の第2種完全楕円方程式により計算されます。
L = 4 \int_0^{\frac{\pi}{2}} \sqrt{a^2 \cos^2 t + b^2 \sin^2 t} dt区分求積法
解析解を求められない積分値を計算する場合には区分求積法を用いて解くことができます。下のコードでは楕円の周長を区分求積法によりn個の区間に分割して求め、その値を返す関数を作成しています。この関数は楕円のパラメータa、bの他に精度を決定するnを引数に取ります。
import math
def calculateC(a, b, n):
result=0.0
for i in range(n):
t=(math.pi*i)/(2*n)
result+=math.pow(a*a*math.cos(t)*math.cos(t)+b*b*math.sin(t)*math.sin(t),0.5)
result*=(math.pi*0.5/n)
result*=4
return result
print(calculateC(6,4,10000))a=6,b=4,n=10000のとき、31.731507という結果が得られました。このような数値計算の結果が正しいのか確かめるにはWolfram Alphaが便利です。Wolfram Alphaの計算結果は31.7309ですので、およそ合致していることがわかります。さらに精度が必要であればnの値を大きくすればよいだけです。
Integrate[4(6^2cos(t)cos(t)+4^2sin(t)*sin(t))^0.5,{t,0,pi/2}] - Wolfram|Alpha
