コンピュータによる定積分の精度の検証

Last updated at 2025-12-02Posted at 2025-12-01

プログラムで定積分を計算するときは、細長い長方形の合計で近似的に計算する(いわゆる区分求積法)。ただ近似なので、例えば手計算で原始関数を求めて具体的な値を代入した値とは僅かにずれが生じるのはいわずもがな。
では両者の誤差がどれほどなのか実際に計算して観察してみる。
分かりやすさのために、以下の条件で検証する。

関数f(x)は常に正になるものにする。
f(x)が正になったり負になったりすると、定積分が足されるときと引かるときがあり、誤差が分かりづらくなるため。
f(x)は大きくなりすぎないものにする。
定積分の範囲を大きくすると値が大きくなりすぎてオーバーフローする可能性があるため。例えばf(x)=x^2とすると、積分値がとても大きくなる場合がある。とはいってもよほど範囲を広げない限り大丈夫だと思うが。

上記2つの条件を満たす関数として、f(x)=1-cos(x)とした。

以下の手順で比較する。

積分範囲は0~x_maxとする。
x_maxは適宜設定。
高さf(x)=1-cos(x)と横dxの長方形を足し合わせた値と、f(x)の原始関数F(x)=x-sin(x)に直接積分範囲の0とx_maxを代入した値(=F(x)-F(0))を比較する(細かい話だが、今回の計算には影響がないため、原始関数において積分定数は省略する)。

    int N;  // 0 <= x <= x_max におけるサンプル数
    int n = 100;    // 1あたりの刻み数
    int x_max = 100;

    // S1 : 足し合わせによる積分値(区分求積法), S2 : 原始関数から算出する積分値
    float S1 = 0.0, S2 = 0.0;
    // 2つの手法で算出した定積分の差分
    float diffS;

    dx = 1.0/(float)n;
    N = n * x_max;

nは範囲1あたりのサンプル数(例えば0<=x<1の範囲でのサンプル数)。Nは積分範囲のサンプル数。上記で書いてある通りN=n*x_maxである。Nの値が中途半端な値にすること誤差が生じることを防ぐために、x_maxは整数に設定する。

以下、区分求積法による計算

    // 細かい長方形の足し合わせで 0 <~ x <= x_max の範囲の定積分を求める
    for(i = 0; i < N; i++){
        x = (float)i * dx;
       S1 = S1 + function(x) * dx;
    }

// 関数f(x) = 1 - cos(x)
float function(float x){
    return(1.0 - cos(x));
}

原始関数による計算。やっていることは手計算で定積分を求めるのと同じ。下記では紛らわしいがxの値はx_maxになる。x_maxはint型なので、float型であるxにx_maxを同じ値を代入した。

    // 原始関数により 0 <= x <= x_max の範囲の定積分を求める
    x = (float)N * dx;
    S2 = primitiveFunction(x) - primitiveFunction(0.0);

// 原始関数F(x) = x - sin(x)
float primitiveFunction(float x){
    return(x - sin(x));
}

上記はfloat型でプログラムを書いたが、double型でも同じ計算をして区分求積法と原始関数による値の誤差の比較をした。より高精度なdouble型のほうが誤差は小さくなると予想できる。
n, N, x_maxを変化させたときの定積分をfloat型、double型のそれぞれで計算する。
どのn, N, x_maxの値で計算したか、floatかdoubleで計算したかは、下記のコンソールでの表示結果を見ればわかると思う。

検証結果

== n = 10, x_max = 100 ==
-- integralFloat.c
N = 1000, n = 10
x_max=100 (N*dx=1000*0.100000=100.000000), dx=0.1000000015
Integral( 0.0000000000 <= x <= 100.0000000000 ) S1 = 100.4990234375
Integral( 0.0000000000 <= x <= 100.0000000000 ) S2 = 100.5063629150
difference between S1 and S2 : S1 - S2 = -0.0073394775

-- integralDouble.c
N = 1000, n = 10
x_max=100 (N*dx=1000*0.100000=100.000000), dx=0.1000000000
Integral( 0.0000000000 <= x <= 100.0000000000 ) S1 = 100.4990595430
Integral( 0.0000000000 <= x <= 100.0000000000 ) S2 = 100.5063656411
difference between S1 and S2 : S1 - S2 = -0.0073060981

== n = 100, x_max = 100 ==
-- integralFloat.c
N = 10000, n = 100
x_max=100 (N*dx=10000*0.010000=100.000000), dx=0.0099999998
Integral( 0.0000000000 <= x <= 100.0000000000 ) S1 = 100.5055618286
Integral( 0.0000000000 <= x <= 100.0000000000 ) S2 = 100.5063629150
difference between S1 and S2 : S1 - S2 = -0.0008010864

-- integralDouble.c
N = 10000, n = 100
x_max=100 (N*dx=10000*0.010000=100.000000), dx=0.0100000000
Integral( 0.0000000000 <= x <= 100.0000000000 ) S1 = 100.5056730158
Integral( 0.0000000000 <= x <= 100.0000000000 ) S2 = 100.5063656411
difference between S1 and S2 : S1 - S2 = -0.0006926254

== n = 1000, x_max = 100 ==
-- integralFloat.c
N = 100000, n = 1000
x_max=100 (N*dx=100000*0.001000=100.000008), dx=0.0010000000
Integral( 0.0000000000 <= x <= 100.0000076294 ) S1 = 100.5052642822
Integral( 0.0000000000 <= x <= 100.0000076294 ) S2 = 100.5063629150
difference between S1 and S2 : S1 - S2 = -0.0010986328

-- integralDouble.c
N = 100000, n = 1000
x_max=100 (N*dx=100000*0.001000=100.000000), dx=0.0010000000
Integral( 0.0000000000 <= x <= 100.0000000000 ) S1 = 100.5062967583
Integral( 0.0000000000 <= x <= 100.0000000000 ) S2 = 100.5063656411
difference between S1 and S2 : S1 - S2 = -0.0000688828

== n = 10, x_max = 1000 ==
-- integralFloat.c
N = 10000, n = 10
x_max=1000 (N*dx=10000*0.100000=1000.000000), dx=0.1000000015
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S1 = 999.1522827148
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S2 = 999.1730957031
difference between S1 and S2 : S1 - S2 = -0.0208129883

-- integralDouble.c
N = 10000, n = 10
x_max=1000 (N*dx=10000*0.100000=1000.000000), dx=0.1000000000
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S1 = 999.1519285944
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S2 = 999.1731204595
difference between S1 and S2 : S1 - S2 = -0.0211918650

== n = 100, x_max = 1000 ==
-- integralFloat.c
N = 100000, n = 100
x_max=1000 (N*dx=100000*0.010000=1000.000000), dx=0.0099999998
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S1 = 999.1686401367
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S2 = 999.1730957031
difference between S1 and S2 : S1 - S2 = -0.0044555664

-- integralDouble.c
N = 100000, n = 100
x_max=1000 (N*dx=100000*0.010000=1000.000000), dx=0.0100000000
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S1 = 999.1709392455
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S2 = 999.1731204595
difference between S1 and S2 : S1 - S2 = -0.0021812139

== n = 1000, x_max = 1000 ==
-- integralFloat.c
N = 1000000, n = 1000
x_max=1000 (N*dx=1000000*0.001000=1000.000061), dx=0.0010000000
Integral( 0.0000000000 <= x <= 1000.0000610352 ) S1 = 999.2051391602
Integral( 0.0000000000 <= x <= 1000.0000610352 ) S2 = 999.1731567383
difference between S1 and S2 : S1 - S2 = 0.0319824219

-- integralDouble.c
N = 1000000, n = 1000
x_max=1000 (N*dx=1000000*0.001000=1000.000000), dx=0.0010000000
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S1 = 999.1729017179
Integral( 0.0000000000 <= x <= 1000.0000000000 ) S2 = 999.1731204595
difference between S1 and S2 : S1 - S2 = -0.0002187415

== n = 10, x_max = 10000 ==
-- integralFloat.c
N = 100000, n = 10
x_max=10000 (N*dx=100000*0.100000=10000.000000), dx=0.1000000015
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S1 = 10000.2119140625
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S2 = 10000.3056640625
difference between S1 and S2 : S1 - S2 = -0.0937500000

-- integralDouble.c
N = 100000, n = 10
x_max=10000 (N*dx=100000*0.100000=10000.000000), dx=0.1000000000
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S1 = 10000.2077518994
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S2 = 10000.3056143889
difference between S1 and S2 : S1 - S2 = -0.0978624895

== n = 100, x_max = 10000 ==
-- integralFloat.c
N = 1000000, n = 100
x_max=10000 (N*dx=1000000*0.010000=10000.000000), dx=0.0099999998
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S1 = 9994.2968750000
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S2 = 10000.3056640625
difference between S1 and S2 : S1 - S2 = -6.0087890625

-- integralDouble.c
N = 1000000, n = 100
x_max=10000 (N*dx=1000000*0.010000=10000.000000), dx=0.0100000000
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S1 = 10000.2958510654
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S2 = 10000.3056143889
difference between S1 and S2 : S1 - S2 = -0.0097633234

== n = 1000, x_max = 10000 ==
-- integralFloat.c
N = 10000000, n = 1000
x_max=10000 (N*dx=10000000*0.001000=10000.000000), dx=0.0010000000
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S1 = 9954.6572265625
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S2 = 10000.3056640625
difference between S1 and S2 : S1 - S2 = -45.6484375000

-- integralDouble.c
N = 10000000, n = 1000
x_max=10000 (N*dx=10000000*0.001000=10000.000000), dx=0.0010000000
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S1 = 10000.3046382862
Integral( 0.0000000000 <= x <= 10000.0000000000 ) S2 = 10000.3056143889
difference between S1 and S2 : S1 - S2 = -0.0009761027

やはりdouble型で計算したほうがfloat型より誤差が小さくなる。
x_max=10000かつ、n=100とn=1000のときでfloat型では大きな差が出ている。より細かく刻んでいるn=1000のほうが誤差が大きいのは謎だが。加えて、fooat型では積分範囲自体にわずかに誤差が生じているケースもある。
一方でdouble型では誤差は非常に小さい。しかし、これが研究などのシミュレーションにおいてどれだけの影響があるのかは分からない。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up