別に求めたくないですがレポートで求めろと言われてしまったので、求めます。

いろいろな方法があると思います。ググるとガウス=ルジャンドル法とかいろいろ出てきます。

実際ガウス=ルジャンドル法を使うことでとても少ない計算回数で出せるのですが、なぜそれをすることで求まるのか理解するのは難しそうです。「理解できないものを使いたくない！」という気持ちになってしまうので、これより容易に理解できる手法でやりたいです。今回は数値積分をして円周率を求めたいと思います。

関数y=\sqrt{1 – x^2}を0から1の区間で定積分をし、4倍することで円周率が出てきます。

区分求積法
\int_a^b f(x) dx = \lim_{n\to \infty} \frac{b – a}{n} \sum_{k = 1}^{n} f(a + (b – a)\frac{k}{n})

を使います。ですが、n\to \inftyにはできないので適当な大きい数にしておきます。

#include <stdio.h>
#include <math.h>

// 積分する関数 f(x)
double f(double x) {
    return sqrt(1 - x * x);
}

// \int_a^b f(x)dx を計算する
double integral(double a, double b) {
    int n = 1 << 25;
    double S = 0;
    for(int k = 1; k < n; k++) {
        S += (b - a) / n * f(a + (b - a) * k / n);
    }
    return S;
}

int main() {

    printf("%.15f", integral(0, 1) * 4);

}

実行結果：
3.141592593977977

ダメです。10桁足りませんね。

nをさらに大きく、2^{30}にしてみても、

3.141592651726471

となってしまいダメです。ちなみにこの計算に2秒ほどかかっています。このプログラムの計算量はO(n)なためこれ以上nを大きくすると時間がとてもかかってしまいよくないです。

誤差の原因は何でしょうか？

実際に図にしてみるとわかりやすいのですが、はみ出してしまっている部分がどうしても出てきてしまいます。それが誤差を生んでしまっているので、うまくこの誤差を小さくしたいです。

底辺が(b-a)/n、高さがf(a + (b – a) k / n)の長方形をたくさん足し合わせることにより値を出しましたが、上底f(a + (b – a) k/n)、下底f(a + (b – a) (k – 1) / n )、高さ(b – a)/ nの台形をたくさん足し合わせれば、より誤差が減ります。

#include <stdio.h>
#include <math.h>

// 積分する関数 f(x)
double f(double x) {
    return sqrt(1 - x * x);
}

// \int_a^b f(x)dx を計算する
double integral(double a, double b) {
    int n = 1 << 25;
    double S = 0;
    for(int k = 1; k < n; k++) {
        S += ((f(a + (b - a) * k / n) + f(a + (b - a) * (k - 1) / n) ) * (b - a) / n) / 2;
    }
    return S;
}

int main() {

    printf("%.15f", integral(0, 1) * 4);

}

実行結果：
3.141592653569314

小数点第10位まで求めることができました！

今回は台形近似までで指定の精度で求めることができましたが、このほかにもシンプソン則を用いて数値積分を行う手法もあります。

インライン数式がうまく入らなくて悲しいです……