数値積分における台形公式の高精度化
数値積分におけるシンプソンの公式やブールの公式（ニュートン・コーツの公式）を，台形公式の様に変形した，台形公式より収束の
速いアルゴリズムを考えまし た。
特徴は，分割数による場合分けが不要であり，台形公式の端点の処理を変えただけで，滑らかな関数に対してシンプソンの公式やブールの公式と同程度に高速に
収束することです。特に，端点の処理だけで収束が速くなるのが不思議であったので，高次の場合にも同様であるか試してみたところ，10次式での近似まで実
際に成り立つことを確かめました。
ただし，高次式での近似が常に速く収束するとは限りません。f(x)=1/SQR(1-(k*SIN(x))^2 )
(k2≪1)を0からπ/2まで積分するのような場合は，変形版は元の台形公式やシンプソンの公式より遅く収束します。また，例えば
f(x)=SQR(1-x^2)
を-1から1まで積分してπ/2を求めるような場合は，端点が微分可能でないため，シンプソンの公式やブールの公式を用いても，高次収束することはありま
せん。
個々の関数に対して最適なアルゴリズムを用いれば良いのですが，高次の公式は低次の公式より，平均すると速く収束します。
さらに，和を求めるのに，順番に足していくのではなく，二分木構造により和を求めることにより，丸め誤差の集積が起こらないようにしました。
解説 Newton-Cotes.pdf
台形公式のプログラムは次の様です。これと，他のプログラムを比較してみて下さい。
LET h=(b-a)/n
LET s=(f(a)+f(b))/2
FOR i=1 TO n-1
LET s=s+f(a+h*i)
NEXT i
LET daikei=s*h
2次式で近似するシンプソンの公式は，次のように変形でき，nは3以上の任意の自然数で計算できます。
LET h=(b-a)/n
LET s=(f(a)+f(b))*3
LET s=s+(f(a+h/2)+f(b-h/2))*4
LET s=s+(f(a+h)+f(b-h))*11
LET s=s/12
FOR i=2 TO n-2
LET s=s+f(a+h*i)
NEXT i
LET simpson=s*h
4次式で近似するブールの公式は，次のように変形でき，nは7以上の任意の自然数で計算できます。
LET h=(b-a)/n
LET s=(f(a)+f(b))*70
LET s=s+(f(a+h/4)+f(b-h/4))*32
LET s=s+(f(a+h/2)+f(b-h/2))*76
LET s=s+(f(a+h*3/4)+f(b-h*3/4))*128
LET s=s+(f(a+h)+f(b-h))*187
LET s=s+(f(a+h*3/2)+f(b-h*3/2))*100
LET s=s+(f(a+h*2)+f(b-h*2))*218
LET s=s+(f(a+h*9/4)+f(b-h*9/4))*96
LET s=s+(f(a+h*3)+f(b-h*3))*353
LET s=s/360
FOR i=4 TO n-4
LET s=s+f(a+h*i)
NEXT i
LET boole=s*h
c言語での1次，2次，4次，6次，8次，10次関数近似公式のコード
#include <stdio.h>
#define f(x) 1/(x*x)
double daikei(double a, double b, int n){
int i;
double h,s,x;
h=(b-a)/n;
s=((f(a))+(f(b)))/2.0;
for(i=1;i<=n-1;i++){
x=a+h*i;
s=s+(f(x));
}
return s*h;
}
double simpson(double a, double b, int n){//nは3以上
int i;
double h,s,x,y;
h=(b-a)/n;
s=((f(a))+(f(b)))/4.0;
x=a+h*0.5;
y=b-h*0.5;
s=s+((f(x))+(f(y)))/3.0;
x=a+h;
y=b-h;
s=s+((f(x))+(f(y)))*11.0/12;
for(i=2;i<=n-2;i++){
x=a+h*i;
s=s+(f(x));
}
return s*h;
}
double boole(double a, double b, int n){//nは7以上
int i;
double h,s,dx,x,y;
double
iti[8]={1.0/4,1.0/2,3.0/4,1.0,3.0/2,2.0,9.0/4,3.0};
double keisuu[8]={32,76,128,187,100,218,96,353};
h=(b-a)/n;
s=((f(a))+(f(b)))*70;
for(i=0;i<=7;i++){
dx=h*iti[i];
x=a+dx;
y=b-dx;
s=s+((f(x))+(f(y)))*keisuu[i];
}
s=s/360;
for(i=4;i<=n-4;i++){
x=a+h*i;
s=s+(f(x));
}
return s*h;
}
double int6(double a, double b, int n){//nは11以上
int i;
double h,s,dx,x,y;
double
iti[17]={1.0/6,1.0/3,1.0/2,2.0/3,5.0/6,1.0,4.0/3,3.0/2,5.0/3,2.0,5.0/2,8.0/3,3.0,10.0/3,4.0,25.0/6,5.0};
double
keisuu[17]={216,459,920,945,1296,2208,162,816,567,2955,2008,108,3459,999,3662,1080,4999};
h=(b-a)/n;
s=((f(a))+(f(b)))*861;
for(i=0;i<=16;i++){
dx=h*iti[i];
x=a+dx;
y=b-dx;
s=s+((f(x))+(f(y)))*keisuu[i];
}
s=s/5040;
for(i=6;i<=n-6;i++){
x=a+h*i;
s=s+(f(x));
}
return s*h;
}
double int8(double a, double b, int n){//nは15以上
int i;
double h,s,dx,x,y;
double
iti[29]={1.0/8,1.0/4,3.0/8,1.0/2,5.0/8,3.0/4,7.0/8,1.0,9.0/8,5.0/4,3.0/2,7.0/4,15.0/8,2.0
,9.0/4,5.0/2,21.0/8,3.0,25.0/8,7.0/2,15.0/4,4.0,35.0/8,9.0/2,5.0,21.0/4,6.0,49.0/8,7.0};
double
keisuu[29]={5888,10848,28160,17156,39936,52608,47104,43213,31488,16352,20940,5280,83968,31410
,60192,19284,91136,103575,52480,-8228,58336,99196,102912,-5568,184153,28832,177718,41216,225811};
h=(b-a)/n;
s=((f(a))+(f(b)))*35604;
for(i=0;i<=28;i++){
dx=h*iti[i];
x=a+dx;
y=b-dx;
s=s+((f(x))+(f(y)))*keisuu[i];
}
s=s/226800;
for(i=8;i<=n-8;i++){
x=a+h*i;
s=s+(f(x));
}
return s*h;
}
double int10(double a, double b, int n){//nは19以上
int i;
double h,s,dx,x,y;
double
iti[41]={1.0/10,1.0/5,3.0/10,2.0/5,1.0/2,3.0/5,7.0/10,4.0/5,9.0/10
,1.0,6.0/5,7.0/5,3.0/2,8.0/5,9.0/5
,2.0,21.0/10,12.0/5,5.0/2,27.0/10,14.0/5
,3.0,16.0/5,7.0/2,18.0/5
,4.0,21.0/5,9.0/2,24.0/5,49.0/10
,5.0,27.0/5,28.0/5
,6.0,63.0/10,32.0/5
,7.0,36.0/5
,8.0,81.0/10
,9.0};
double
keisuu[41]={106300,164075,591300,67600,958868,776475,1016500,86675,1880200
,1851848,-504300,205125,2644104,-1527450,628625
,1177276,2724000,-571875,2136840,2770500,-734250
,4772079,-2278500,4353576,-3483050
,4097507,-189450,4377812,-2375550,1906800
,5210935,-1707150,1839525
,2621502,3195700,-388200
,5361569,413675
,4892386,956700
,5971453};
h=(b-a)/n;
s=((f(a))+(f(b)))*883685;
for(i=0;i<=40;i++){
dx=h*iti[i];
x=a+dx;
y=b-dx;
s=s+((f(x))+(f(y)))*keisuu[i];
}
s=s/5987520;
for(i=10;i<=n-10;i++){
x=a+h*i;
s=s+(f(x));
}
return s*h;
}
int main(){
int i,n=20;
double a=0.5,b=1.5,s=4.0/3;
for(i=0;i<=10;i++){
printf("分割数=%5d 誤差=%24.20f
%24.20f %24.20f %24.20f %24.20f
%24.20f\n",n,daikei(a,b,n)-s,simpson(a,b,n)-s,boole(a,b,n)-s,int6(a,b,n)-s,int8(a,b,n)-s,int10(a,b,n)-s);
n=n*2;
}
}
c言語での10次関数近似公式の通常のコード
#include <stdio.h>
#define f(x) 1/(x*x)
double int10(double a, double b, int n){
double h,s,t,x,y;
int i;
h=(b-a)/n;
s=(f(a))+(f(b));
t=0;
for(i=10;i<n;i=i+10){
x=a+h*i;
t=t+(f(x));
}
s=(s+t*2)*16067;
t=0;
for(i=1;i<n;i=i+10){
x=a+h*i;
y=x+h*8;
t=t+(f(x))+(f(y));
}
s=s+t*106300;
t=0;
for(i=2;i<n;i=i+10){
x=a+h*i;
y=x+h*6;
t=t+(f(x))+(f(y));
}
s=s+t*(-48525);
t=0;
for(i=3;i<n;i=i+10){
x=a+h*i;
y=x+h*4;
t=t+(f(x))+(f(y));
}
s=s+t*272400;
t=0;
for(i=4;i<n;i=i+10){
x=a+h*i;
y=x+h*2;
t=t+(f(x))+(f(y));
}
s=s+t*(-260550);
t=0;
for(i=5;i<n;i=i+10){
x=a+h*i;
t=t+(f(x));
}
s=s+t*427368;
return s*h*5.0/299376;
}
int main(){
int i,n=20;//nは10の倍数
for(i=1;i<=10;i++){
printf("分割数=%5d
誤差=%24.20f\n",n,int10(0.5,1.5,n)-4.0/3);
n=n*2;
}
}
数値積分における丸め誤差の集積の解消法
数値積分において，分割数nを大きくすると，丸め誤差のために，かえって誤差が増加することが知られています。しかし，それは，数を次々に加算することに
よって起こることであり，和を取る順序を変えれれば丸め誤差は解消することができます。それを示すのが以下のプログラムであり，二分木構造により和を求め
るものです。ただし，処理速度が遅いのが玉に瑕です。
ポイントは，次々に加算すると，巨大な数と小さな数を足すことになり，このとき，小さな数の有効数字が小さなものとして足され，丸め誤差が大きく作用しま
す。これを解消するには，常に同じくらいの大きさの数同士の和を取れば，丸め誤差は小さくなります。それを実現するのが二分木構造により2つの数の和を求
めることであり，このとき，和を取る2つの数の大きさは常に同じオーダーになります。
さらに，数値積分においては，最後に和を取った個数nで割るため，丸め誤差の集積は起こらなくなります。そのため，sin(x)/(cos(x*x)+1
+1.0/1024) のような特異な関数を0から10まで積分するような場合でも，分割数を十分に大きくできるので，積分値が高精度に求まります。
#include <stdio.h>
#include <math.h>
#define f(x) 4/(1+x*x)
double daikei0(double a, double b, int n){
int i;
double h,s,x;
h=(b-a)/n;
s=((f(a))+(f(b)))/2.0;
for(i=1;i<=n-1;i++){
x=a+h*i;
s=s+(f(x));
}
return s*h;
}
double wa(double a,double h,int k, int n){
if (n>1){
int m=(n+1)>>1;
return
wa(a,h,k,m)+wa(a,h,k+m,n>>1);
}else{
double x=a+h*k;
return f(x);
}
}
double daikei(double a,double b,int n){
double s,h=(b-a)/n;
s=(f(a)+f(b))/2;
s=s+wa(a+h,h,0,n-1);
return s*h;
}
int main(){
double a=0,b=1,t;
int n=1000000000;
t=daikei0(a,b,n);
printf("誤差=%23.20f\n",t-M_PI);
t=daikei(a,b,n);
printf("誤差=%23.20f\n",t-M_PI);
}
高精度に出力する実用コード
以
下のコードは，sqrt(x)，sqrt(1-x*x)，1/sqrt(x)，1/sqrt(1-x*x))などの関数も積分できるように，置換積分を用
いています。区間[0，1]における sin(x)/x，sqrt(sin(x))*cos(x)，log(sin(x))*cos(x)
などの積分や，区間[0，10]における sin(x)/(cos(x*x)+1+1.0/1024) の積分も，高精度かつ高速に行います。
ただし，1/sqrt(1-x*x)の積分は，グラフをx軸方向に1だけ平行移動して，1/sqrt(x*(2-x))とし，0から1まで積分します。す
なわち，±∞に発散する特異点は，積分の分割，平行移動により，積分区間の端点，かつ，x=0の位置にしかないようにします。また，1/sqrt(x)
より発散が速い関数には対応していません。
#include <stdio.h>
#include <math.h>
double g(double x){
return 1/sqrt(x*(2-x));
}
double f(double a,double ab,double t){
double y,tt;
tt=t*t;
y=a+(((t-70.0/20)*t+84.0/20)*t-35.0/20)*tt*tt*ab;
tt=tt-t;
return g(y)*tt*tt*tt;
}
double tanten(double a,double b,double ab,double h){
int i;
double s=0,t;
static double iti[41]={1.0/10,1.0/5,3.0/10,2.0/5,1.0/2,3.0/5,7.0/10,4.0/5,9.0/10
,1.0,6.0/5,7.0/5,3.0/2,8.0/5,9.0/5
,2.0,21.0/10,12.0/5,5.0/2,27.0/10,14.0/5
,3.0,16.0/5,7.0/2,18.0/5
,4.0,21.0/5,9.0/2,24.0/5,49.0/10
,5.0,27.0/5,28.0/5
,6.0,63.0/10,32.0/5
,7.0,36.0/5
,8.0,81.0/10
,9.0};
static double keisuu[41]={106300,164075,591300,67600,958868,776475,1016500,86675,1880200
,1851848,-504300,205125,2644104,-1527450,628625
,1177276,2724000,-571875,2136840,2770500,-734250
,4772079,-2278500,4353576,-3483050
,4097507,-189450,4377812,-2375550,1906800
,5210935,-1707150,1839525
,2621502,3195700,-388200
,5361569,413675
,4892386,956700
,5971453};
//s=(f(a,ab,0)+f(b,-ab,0))*883685;
for(i=0;i<=40;i++){
t=h*iti[i];
s=s+(f(a,ab,t)+f(b,-ab,t))*keisuu[i];
}
return s/5987520;
}
double wa(double a,double ab,double h,int k, int n){
if (n>1){
int m=(n+1)>>1;
return wa(a,ab,h,k,m)+wa(a,ab,h,k+m,n>>1);
}else{
return f(a,ab,h*k);
}
}
double int10(double a, double b, int n){
double h=1.0/n,ab=(b-a)*(-20),s;
int m=(n-19)/2;
s=wa(a,ab,h,10,m)+wa(b,-ab,h,10,n-19-m)+tanten(a,b,ab,h);
return s*(b-a)*(-140)*h;
}
double integral(double a,double b){
int i,n;
double s1,s2,r;
n=80;
s1=int10(a,b,n);
for(i=1;i<=24;i++){
//printf("n=%10d 積分値= %4.20f\n",n,s1);
n=n*2;
s2=int10(a,b,n);
r=(s2-s1)/s1;
if (r<0)r=-r;
if (r<1e-14) break;
s1=s2;
}
return s2;
}
int main(){
double a,b,s;
a=0;
b=1;
s=integral(a,b);
printf("積分値=%.20f\n",s);
}
二重指数関数型数値積分公式(DE公式)によるコード (有限区間)
比較のために書いてみました。発散が速いため上記のコードでは計算が難しい，1/cbrt(x*x)の0から1の積分も計算できます。しかし，sinh(x)，cosh
(x)，tanh(x)を計算するため速度は遅くなります。
ただし，機能が単純であるためか http://www.kurims.kyoto-u.ac.jp/~ooura/intde-j.html
にあるパッケージよりもかなり速いようです。
#include <stdio.h>
#include <math.h>
#define RANGE 4.6
double g(double x){
return 1/sqrt(x*(2-x));
}
double f(double a,double ab,double t){
double x,sh,ey;
sh=sinh(t);
ey=exp(M_PI*sh);
x=a+ab/(ey+1);
return g(x)*sqrt(1+sh*sh)/(ey+1/ey+2);
}
double wa(double a,double ab,double h,int k, int n){
if (n>1){
int m=(n+1)>>1;
return wa(a,ab,h,k,m)+wa(a,ab,h,k-m,n>>1);
}else{
return f(a,ab,h*k);
}
}
double int10(double a, double b, int n){
double h=RANGE*2/n,ab=b-a,s;
int m=n/2;
s=wa(a,ab,h,m,m+1)+wa(b,-ab,h,m,n-m);
return s*ab*M_PI*h;
}
double integral(double a,double b){
int i,n;
double s1,s2,r;
n=40;
s1=int10(a,b,n);
for(i=1;i<=25;i++){
//printf("n=%10d 積分値= %4.20f\n",n,s1);
n=n*2;
s2=int10(a,b,n);
r=(s2-s1)/s1;
if (r<0)r=-r;
if (r<1e-14) break;
s1=s2;
}
return s2;
}
int main(){
double a,b,s;
a=0;
b=1;
s=integral(a,b);
printf("積分値=%.20f\n",s);
}
二重指数関数型数値積分公式によるコード (aから∞の区間)
#include <stdio.h>
#include <math.h>
#define RANGE 4.0
double g(double x){
return 1/(sqrt(x)*(1+x));
}
double f(double a,double t){
double x,sh,y;
sh=sinh(t);
y=exp(M_PI*sh);
return g(a+y)*y*sqrt(1+sh*sh);
}
double wa(double a,double h,int k, int n){
if (n>1){
int m=(n+1)>>1;
return wa(a,h,k,m)+wa(a,h,k+m,n>>1);
}else{
return f(a,h*k);
}
}
double int10(double a,int n){
double h=RANGE*2/n;
return wa(a,h,-n/2,n)*h*M_PI;
}
double integral(double a){
int i,n;
double s1,s2,r;
n=40;
s1=int10(a,n);
for(i=1;i<=25;i++){
//printf("n=%10d 積分値= %4.20f\n",n,s1);
n=n*2;
s2=int10(a,n);
r=(s2-s1)/s1;
if (r<0)r=-r;
if (r<1e-14) break;
s1=s2;
}
return s2;
}
int main(){
double a,s;
a=0;
s=integral(a);
printf("積分値=%.20f\n",s);
}
二重指数関数型数値積分公式によるコード (－∞から∞の区間)
#include <stdio.h>
#include <math.h>
#define RANGE 4.6
double g(double x){
return 1/(1+x*x);
}
double f(double t){
double x,sh,y;
sh=sinh(t);
y=sinh(M_PI_2*sh);
return g(y)*sqrt(1+y*y)*sqrt(1+sh*sh);
}
double wa(double h,int k, int n){
if (n>1){
int m=(n+1)>>1;
return wa(h,k,m)+wa(h,k+m,n>>1);
}else{
return f(h*k);
}
}
double int10(int n){
double h=RANGE*2/n;
return wa(h,-n/2,n)*M_PI_2*h;
}
double integral(){
int i,n;
double s1,s2,r;
n=40;
s1=int10(n);
for(i=1;i<=25;i++){
//printf("n=%10d 積分値= %4.20f\n",n,s1);
n=n*2;
s2=int10(n);
r=(s2-s1)/s1;
if (r<0)r=-r;
if (r<1e-14) break;
s1=s2;
}
return s2;
}
int main(){
double s;
s=integral();
printf("積分値=%.20f\n",s);
}