数値積分について
ニートンコーツ法、チェビシェフ公式、ガウス公式を用いた計算
<<<program>>>
<stdio.h>
<math.h>
double fun( double x)
{
//return x*exp(x);
return ((x+2)*x+3)*x+4;
}
double scale( double x, double x0, double x1 )
{
return (x0 * (1.0 - x) + x1 * (1.0 + x)) / 2.0;
}
double newton_cote(int n,double xfrom_local,double xto_local,int xdivision){
int xidx;
double x1,x2,x3,x4,sum;
sum=0;
if(n==1){
for (xidx=1; xidx<xdivision; xidx++) {
x1 = (xfrom_local * (xdivision - xidx) + xto_local * xidx) / xdivision;
sum += fun( x1 );
}
sum *= 2;
sum += fun( xfrom_local ) + fun( xto_local );
sum = (xto_local - xfrom_local) * sum / ( 2.0 * xdivision );
}else if(n==2){
if(xdivision%2==0){
for (xidx=0; xidx<=((xdivision/2)-1); xidx++) {
x1 = (xfrom_local * (xdivision - 2*xidx) + xto_local * 2*xidx)
/ xdivision;
x2 = (xfrom_local * (xdivision - (2*xidx+1)) + xto_local * (2*xidx+1))
/ xdivision;
x3 = (xfrom_local * (xdivision - (2*xidx+2)) + xto_local * (2*xidx+2))
/ xdivision;
sum += (fun(x1)+4*fun(x2)+fun(x3));
}
sum = (xto_local - xfrom_local) * sum / ( 3.0 * xdivision);
}else{
sum = -1;//error
}
}else if(n==3){
if(xdivision%3==0){
for (xidx=0; xidx<=((xdivision/3)-1); xidx++) {
x1 = (xfrom_local * (xdivision - 3*xidx) + xto_local * 3*xidx)
/ xdivision;
x2 = (xfrom_local * (xdivision - (3*xidx+1)) + xto_local * (3*xidx+1))
/ xdivision;
x3 = (xfrom_local * (xdivision - (3*xidx+2)) + xto_local * (3*xidx+2))
/ xdivision;
x4 = (xfrom_local * (xdivision - (3*xidx+3)) + xto_local * (3*xidx+3))
/ xdivision;
sum += (fun(x1)+3*fun(x2)+3*fun(x3)+fun(x4));
}
sum = 3.0 * (xto_local - xfrom_local) * sum / ( 8.0 * xdivision);
}else{
sum = -1;//error
}
}
return sum;
}
double chebyshev(int n,double xfrom_local,double xto_local,int xdivision){
double x0, x1, y;
int xidx;
xdivision=(int)(xdivision/3);//hikakunotame
y = 0;
x0 = xfrom_local;
for( xidx=1; xidx<=xdivision; xidx++ ) {
x1 = (xfrom_local * (xdivision - xidx) + xto_local * xidx) / xdivision;
y += ( fun( scale( -1.0/sqrt(2.0), x0, x1 ) )
+ fun( scale( 0.0,x0, x1 ) )
+ fun( scale( +1.0/sqrt(2.0), x0, x1 ) )
) * (x1 - x0) / 3.0;
x0 = x1;
}
return y;
}
double gauss(int n,double xfrom_local,double xto_local,int xdivision){
double x0, x1, y;
int xidx;
xdivision=(int)(xdivision/3);//hikakunotame
y = 0;
x0 = xfrom_local;
for( xidx=1; xidx<=xdivision; xidx++ ) {
x1 = (xfrom_local * (xdivision - xidx) + xto_local * xidx) / xdivision;
y += ( (5.0/9.0) * fun( scale( -sqrt(3.0/5.0), x0, x1 ) )
+ (8.0/9.0) * fun( scale(0.0, x0, x1 ) )
+ (5.0/9.0) * fun( scale( +sqrt(3.0/5.0), x0, x1 ) )
) * (x1 - x0) / 2.0;
x0 = x1;
}
return y;
}
int main(void){
double ans;
double xfrom, xto,true_value;
int xdivision;
//init
//true_value=2/exp(1);
true_value=28.0/3.0;
xfrom=-1;
xto=1;
printf("division number(xdivision) > ");
scanf( "%d",&xdivision );
do{
printf( "true value:\t%20.15lf\n", true_value );
ans=newton_cote(1,xfrom,xto,xdivision);
printf( "newton_cote(1):\t%20.15lf\n", ans );
ans=newton_cote(2,xfrom,xto,xdivision);
printf( "newton_cote(2):\t%20.15lf\n", ans );
ans=newton_cote(3,xfrom,xto,xdivision);
printf( "newton_cote(3):\t%20.15lf\n", ans );
ans=chebyshev(3,xfrom,xto,xdivision);
printf( "chebyshev(3):\t%20.15lf\n", ans );
ans=gauss(3,xfrom,xto,xdivision);
printf( "gauss(3):\t%20.15lf\n", ans );
printf("division number(xdivision) (if(xdivision<=0){END}) > ");
scanf( "%d",&xdivision );
}while(xdivision>0);
return 0;
}
<<<results-1(x*exp(x))>>>与えられた課題のsimpsonの公式のみを出力したもの
division number(xdivision) > 2
newton_cote(2): 0.783467462429201
division number(xdivision) (if(xdivision<=0){END}) > 4
newton_cote(2): 0.739130601543765
division number(xdivision) (if(xdivision<=0){END}) > 6
newton_cote(2): 0.736440917298669
division number(xdivision) (if(xdivision<=0){END}) > 8
newton_cote(2): 0.735976501958868
division number(xdivision) (if(xdivision<=0){END}) > 10
newton_cote(2): 0.735848367669737
division number(xdivision) (if(xdivision<=0){END}) > 20
newton_cote(2): 0.735764504414181
division number(xdivision) (if(xdivision<=0){END}) > 40
newton_cote(2): 0.735759234181559
division number(xdivision) (if(xdivision<=0){END}) > 80
newton_cote(2): 0.735758904339986
division number(xdivision) (if(xdivision<=0){END}) > 200
newton_cote(2): 0.735758882906062
division number(xdivision) (if(xdivision<=0){END}) > 500
newton_cote(2): 0.735758882357302
division number(xdivision) (if(xdivision<=0){END}) > 1000
newton_cote(2): 0.735758882343786
division number(xdivision) (if(xdivision<=0){END}) > 2000
newton_cote(2): 0.735758882342941
division number(xdivision) (if(xdivision<=0){END}) > 5000
newton_cote(2): 0.735758882342885
division number(xdivision) (if(xdivision<=0){END}) > 10000
newton_cote(2): 0.735758882342883
division number(xdivision) (if(xdivision<=0){END}) > 30000
newton_cote(2): 0.735758882342885
division number(xdivision) (if(xdivision<=0){END}) > 80000
newton_cote(2): 0.735758882342884
division number(xdivision) (if(xdivision<=0){END}) > 100000
newton_cote(2): 0.735758882342886
division number(xdivision) (if(xdivision<=0){END}) >0
<<<results-2(x*exp(x))>>>真の値、台形公式、simpsonの公式、simpsonの3/8公式、チェビシェフ、ガウスの順に計算したもの
division number(xdivision) > 6
true value: 0.735758882342885
newton_cote(1): 0.785924282165319
newton_cote(2): 0.736440917298669
newton_cote(3): 0.737264668230509
chebyshev(3): 0.734912297016692
gauss(3): 0.735751841168723
division number(xdivision) (if(xdivision<=0){END}) > 12
true value: 0.735758882342885
newton_cote(1): 0.748332666995611
newton_cote(2): 0.735802128605708
newton_cote(3): 0.735855719173966
chebyshev(3): 0.735704418254167
gauss(3): 0.735758768617193
division number(xdivision) (if(xdivision<=0){END}) > 30
true value: 0.735758882342885
newton_cote(1): 0.737772146364007
newton_cote(2): 0.735759993950277
newton_cote(3): 0.735761381524993
chebyshev(3): 0.735757476580191
gauss(3): 0.735758881872675
division number(xdivision) (if(xdivision<=0){END}) > 60
true value: 0.735758882342885
newton_cote(1): 0.736262250485027
newton_cote(2): 0.735758951858701
newton_cote(3): 0.735759038723198
chebyshev(3): 0.735758794379384
gauss(3): 0.735758882335528
division number(xdivision) (if(xdivision<=0){END}) > 300
true value: 0.735758882342885
newton_cote(1): 0.735779017736025
newton_cote(2): 0.735758882454131
newton_cote(3): 0.735758882593186
chebyshev(3): 0.735758882202090
gauss(3): 0.735758882342884
division number(xdivision) (if(xdivision<=0){END}) > 600
true value: 0.735758882342885
newton_cote(1): 0.735763916196385
newton_cote(2): 0.735758882349838
newton_cote(3): 0.735758882358529
chebyshev(3): 0.735758882334085
gauss(3): 0.735758882342885
division number(xdivision) (if(xdivision<=0){END}) > 3000
true value: 0.735758882342885
newton_cote(1): 0.735759083697091
newton_cote(2): 0.735758882342896
newton_cote(3): 0.735758882342909
chebyshev(3): 0.735758882342871
gauss(3): 0.735758882342885
division number(xdivision) (if(xdivision<=0){END}) > 6000
true value: 0.735758882342885
newton_cote(1): 0.735758932681438
newton_cote(2): 0.735758882342885
newton_cote(3): 0.735758882342886
chebyshev(3): 0.735758882342885
gauss(3): 0.735758882342886
division number(xdivision) (if(xdivision<=0){END}) > 0
<<<results-3(((x+2)*x+3)*x+4)>>>自分で設定した被積分関数について
真の値、台形公式、simpsonの公式、simpsonの3/8公式、チェビシェフ、ガウスの順に計算したもの
division number(xdivision) > 6
true value: 9.333333333333334
newton_cote(1): 9.407407407407407
newton_cote(2): 9.333333333333334
newton_cote(3): 9.333333333333334
chebyshev(3): 9.333333333333334
gauss(3): 9.333333333333332
division number(xdivision) (if(xdivision<=0){END}) > 12
true value: 9.333333333333334
newton_cote(1): 9.351851851851849
newton_cote(2): 9.333333333333334
newton_cote(3): 9.333333333333332
chebyshev(3): 9.333333333333334
gauss(3): 9.333333333333332
division number(xdivision) (if(xdivision<=0){END}) > 30
true value: 9.333333333333334
newton_cote(1): 9.336296296296297
newton_cote(2): 9.333333333333334
newton_cote(3): 9.333333333333334
chebyshev(3): 9.333333333333334
gauss(3): 9.333333333333334
division number(xdivision) (if(xdivision<=0){END}) > 60
true value: 9.333333333333334
newton_cote(1): 9.334074074074074
newton_cote(2): 9.333333333333334
newton_cote(3): 9.333333333333334
chebyshev(3): 9.333333333333336
gauss(3): 9.333333333333334
division number(xdivision) (if(xdivision<=0){END}) > 300
true value: 9.333333333333334
newton_cote(1): 9.333362962962966
newton_cote(2): 9.333333333333334
newton_cote(3): 9.333333333333330
chebyshev(3): 9.333333333333334
gauss(3): 9.333333333333332
division number(xdivision) (if(xdivision<=0){END}) > 600
true value: 9.333333333333334
newton_cote(1): 9.333340740740741
newton_cote(2): 9.333333333333327
newton_cote(3): 9.333333333333336
chebyshev(3): 9.333333333333334
gauss(3): 9.333333333333334
division number(xdivision) (if(xdivision<=0){END}) > 3000
true value: 9.333333333333334
newton_cote(1): 9.333333629629646
newton_cote(2): 9.333333333333323
newton_cote(3): 9.333333333333343
chebyshev(3): 9.333333333333330
gauss(3): 9.333333333333330
division number(xdivision) (if(xdivision<=0){END}) > 6000
true value: 9.333333333333334
newton_cote(1): 9.333333407407418
newton_cote(2): 9.333333333333341
newton_cote(3): 9.333333333333332
chebyshev(3): 9.333333333333327
gauss(3): 9.333333333333327
division number(xdivision) (if(xdivision<=0){END}) > 0
<<<考察１（精度について考える）>>>
シンプソンの公式で得られた結果をまとめると以下のようになる。
newton_cote(2): 0.783467462429201[2]
newton_cote(2): 0.739130601543765[4]
newton_cote(2): 0.736440917298669[6]
newton_cote(2): 0.735976501958868[8]
newton_cote(2): 0.735848367669737[10]
newton_cote(2): 0.735764504414181[20]
newton_cote(2): 0.735759234181559[40]
newton_cote(2): 0.735758904339986[80]
newton_cote(2): 0.735758882906062[200]
newton_cote(2): 0.735758882357302[500]
newton_cote(2): 0.735758882343786[1000]
newton_cote(2): 0.735758882342941[2000]
newton_cote(2): 0.735758882342885[5000]
newton_cote(2): 0.735758882342883[10000]
newton_cote(2): 0.735758882342885[30000]
newton_cote(2): 0.735758882342884[80000]
newton_cote(2): 0.735758882342886[100000]
true value: 0.735758882342885
正しい桁数まで表示してみると・・・
newton_cote(2): 0.7[2]
newton_cote(2): 0.73[4]
newton_cote(2): 0.73[6]
newton_cote(2): 0.735[8]
newton_cote(2): 0.735[10]
newton_cote(2): 0.7357[20]
newton_cote(2): 0.73575[40]
newton_cote(2): 0.735758[80]
newton_cote(2): 0.735758882[200]
newton_cote(2): 0.7357588823[500]
newton_cote(2): 0.73575888234[1000]
newton_cote(2): 0.735758882342[2000]
newton_cote(2): 0.735758882342885[5000]
newton_cote(2): 0.73575888234288[10000]
newton_cote(2): 0.735758882342885[30000]
newton_cote(2): 0.73575888234288[80000]
newton_cote(2): 0.73575888234288[100000]
true value: 0.735758882342885
これより、全桁正しいのは5000分割したときなのだが、それ以降増やしていくと最後の桁がずれていることに気づく。求めたい精度の桁数に対して、適切な分割数があると考えられる。
また、分割数を１０倍したもののみ見てみると・・・
newton_cote(2): 0.7[2] 有効数字１桁
newton_cote(2): 0.7357[20] 有効数字４桁
newton_cote(2): 0.735758882[200] 有効数字９桁（８桁）
newton_cote(2): 0.735758882342[2000] 有効数字１２桁
分割数200のとき最後の有効が偶然合ったとして有効数字８桁とすると・・・
分割数を10倍するごとに有効数字が４桁多くなっていることがわかる。（だいたい）
文献によると・・シンプソンの公式（閉じた式）の誤差項は-h^5*f(4)/90とのことなので、誤差項のオーダーがh^5によると簡易化して考えてみる。
h=(xto-xfrom)/xdivisionよりhは分割数に反比例するので、分割数が10倍になると、オーダーが5桁上がるはずだが、実際のところ４桁しか上がっていない。
これの原因として、上の誤差項は閉じた式の誤差項なので、実際のプログラムではその誤差項をfor文のなかで、((xdivision/2)-1)回繰り返していることから、誤差が((xdivision/2)-1)回たまるとすると（h^5*xdivisionと簡易化して）、オーダーはだいたいh^4となっているのではないか。
<<<考察２（関数が違った場合について考える）>>>
被積分関数が違った場合・・・（今回は以下の場合）
return x*exp(x);
return ((x+2)*x+3)*x+4;
このとき、分割数を６にしたときから、台形公式以外で求めた積分値はほぼ真の値を示している。これは被積分関数が、ｘの多項式であるからかもしれない（実験不足で自信がありません）アルゴリズム的にラグランジュ補間式を用いているからであるかもしれない。
<<<考察３おまけ（区間毎に Chebyshev 公式や Gauss 公式を用いた結果とも比較 ）>>>
結果2を必要な公式で計算したもののみ取り出し、真の値と同じ桁のみ残したものが以下である。プログラム中のscale関数でチェビシェフとガウスは分割しているとみなして、同じ分割数にするため、xdivision =xdivision/3としている。これで比較してみると・・・
division number(xdivision) > 6
newton_cote(3): 0.73
chebyshev(3): 0.73
gauss(3): 0.73575
division number(xdivision) (if(xdivision<=0){END}) > 12
newton_cote(3): 0.735
chebyshev(3): 0.7357
gauss(3): 0.735758
division number(xdivision) (if(xdivision<=0){END}) > 30
newton_cote(3): 0.7357
chebyshev(3): 0.73575
gauss(3): 0.73575888
division number(xdivision) (if(xdivision<=0){END}) > 60
newton_cote(3): 0.73575
chebyshev(3): 0.735758
gauss(3): 0.7357588823
division number(xdivision) (if(xdivision<=0){END}) > 300
newton_cote(3): 0.735758882
chebyshev(3): 0.735758882
gauss(3): 0.73575888234288
division number(xdivision) (if(xdivision<=0){END}) > 600
newton_cote(3): 0.7357588823
chebyshev(3): 0.7357588823
gauss(3): 0.735758882342885
division number(xdivision) (if(xdivision<=0){END}) > 3000
newton_cote(3): 0.735758882342
chebyshev(3): 0.7357588823428
gauss(3): 0.735758882342885
division number(xdivision) (if(xdivision<=0){END}) > 6000
newton_cote(3): 0.73575888234288
chebyshev(3): 0.73575888234288
gauss(3): 0.73575888234288
ガウス公式は収束がはやい。
シンプソンの3/8公式とチェビシェフ公式はほぼ同じくらいで若干チェビシェフ公式が早い。
（シンプソンの3/8公式よりこれだけ収束が早いということは何か比較を間違えているかもしれない・・・あるいは、たまたま早く求まる関数だったのかもしれない。）
最終更新：2011年01月18日 16:04
