| 場所 | 正 | 誤 |
|---|---|---|
| p. 62左 下から7行目 | p_0((a+b)/2) | p_0((z+b)/2) |
| p. 64右 下から13行目 | 削除 | および(a,b) |
| p. 65左 下から11行目 | f(tanh | (tanh |
表2の出力に用いたCプログラム(文字コードはUTF-8)
/* Double exponential formula */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX 100
#define EPS1 1.0e-16
#define EPS2 1.0e-12
#define PI 3.141592653589793238
double interval(double *xa, double *xb);
double trapezoidal(double h, double *xa, double *xb, double s);
double func(double x);
int main(void)
{
double *xa, *xb, s = 0.0, h = 1.0, s0 = 0.0, I;
xa = malloc(sizeof(double));
xb = malloc(sizeof(double));
I = 4.442882938158366101788487867452204227448; /* B(1/4,3/4) */
s = interval(xa,xb);
printf("[%5.3e,%5.3e]\n",*xa,*xb);
printf("%15.8e %20.15f %15.8e %15.8e\n",h,s,s-I,s-s0);
while(fabs(s-s0) > EPS2){
s0 = s;
h *= 0.5;
s = trapezoidal(h,xa,xb,s);
printf("%15.8e %20.15f %15.8e %15.8e\n",h,s,s-I,s-s0);
}
free(xa);
free(xb);
return 0;
}
double interval(double *xa, double *xb)
{
int m, n;
double x, fx, fx0, fxp, w, h = 1.0, s;
x = 0.0;
fx = func(x);
fx0 = fx;
fxp = fx;
x += h;
fx = func(x);
n = 1;
s = h * (fxp + fx);
w = fabs(fxp) + fabs(fx);
/* determination of upper boundary */
while(w >EPS1 && n < MAX){
fxp = fx;
x += h;
fx = func(x);
// printf("%10.3f %10.3e %20.15g\n",x,fx,1.0-tanh((PI/2.0)*sinh(x)));
n++;
s += h * fx;
w = fabs(fxp) + fabs(fx);
}
*xb = h * (n - 1);
fxp = fx0;
x = -h;
fx = func(x);
m = 1;
s += h * fx;
w = fabs(fxp) + fabs(fx);
/* determination of lower boundary */
while(w >EPS1 && m < MAX){
x -= h;
fxp = fx;
fx = func(x);
// printf("%10.3f %10.3e %20.15g\n",x,fx,tanh((PI/2.0)*sinh(x))+1.0);
m++;
s += h * fx;
w = fabs(fxp) + fabs(fx);
}
*xa = -h * (m - 1);
return s;
}
double trapezoidal(double h, double *xa, double *xb, double s)
{
double x, w, fx;
x = *xa + h;
w = 0.0;
/* midpoint rule */
while(x < *xb){
fx = func(x);
w += fx;
x += 2 * h;
}
s = 0.5 * s + h * w;
return s;
}
double func(double x)
{
return (PI/2.0)*cosh(x)/(sqrt(exp((PI/2.0)*sinh(x)))*cosh((PI/2.0)*sinh(x)));
}