楕円は円錐を平面で切ったときにできる図形のうちの1つです。今回は楕円の周の長さCを数値計算によって求めます。長軸2a、単軸2bのx2/a2+y2/b2=1で表される楕円の面積はπabと簡単に計算できますが、楕円の周の長さを普通の関数で表すことはできません。ゆえに数値計算が必要になります。周の長さは以下の第2種完全楕円方程式により計算されます。
区分求積法
解析解を求められない積分値を計算する場合には区分求積法を用いて解くことができます。下のコードでは楕円の周長を区分求積法によりn個の区間に分割して求め、その値をdouble型で返す関数を作成しています。この関数は楕円のパラメータa、bの他に精度を決定するnを引数に取ります。
#include <stdio.h>
#include <math.h>
#define PI 3.141592
double calculateC(double a, double b, int n){
double result=0.0;
for(int i=0; i<n; i++){
double t=(PI*i)/(2*n);
result+=pow(a*a*cos(t)*cos(t)+b*b*sin(t)*sin(t),0.5);
}
result*=(PI*0.5/n);
result*=4;
return result;
}
int main(){
printf("%f\n", calculateC(6,4,10000) );
return 0;
}
a=6,b=4,n=10000のとき、31.731502という結果が得られました。このような数値計算の結果が正しいのか確かめるには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
