1 nint6.c

/*
 * nint6.c --- DE 公式
 *
 *         1      2 (1/2)
 *  I1 = ∫   (1-x )     dx = π/2
 *        -1
 *
 *         1      2 (-1/2)
 *  I2 = ∫   (1-x )      dx = π
 *        -1
 *
 *  いずれも端点に特異性があって、古典的な数値積分公式はうまく行かない。
 *  double exponential formula (DE公式) ならばうまく計算できる。
 *
 *  コンパイル: gcc -o nint6 nint6.c -lm
 *
 *  学生 (横山和正君) の指摘により少し修正を加える (2004/1/15)。
 */

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

typedef double rrfunction(double);

double debug = 0;
double pi, halfpi;

/* φ */
double phi(double t)
{
  return tanh(halfpi * sinh(t));
}

/* 2乗 */
double sqr(double x) { return x * x; }

/* φ' */
double dphi(double t)
{
  return halfpi * cosh(t) / sqr(cosh(halfpi * sinh(t)));
}

/* DE 公式による (-1,1) における定積分の計算 */
double de(rrfunction f, double h, double N)
{
  int n;
  double t, S, dS;
  S = f(phi(0.0)) * dphi(0.0);
  for (n = 1; n <= N; n++) {
    t = n * h;
    dS = f(phi(t)) * dphi(t) + f(phi(-t)) * dphi(- t);
    S += dS;
    if (fabs(dS) < 1.0e-16) {
      if (debug)
        printf("\tde(): n=%d, |t|=%g, fabs(dS)=%g\n", n, t, fabs(dS));
      break;
    }
  }
  return h * S;
}

/* テスト用の被積分関数 その1 */
double f1(double x)
{
  return sqrt(1 - x * x);
}

/* テスト用の被積分関数 その2 */
double f2(double x)
{
  if (x >= 1.0 || x <= -1.0) return 0; else return 1 / sqrt(1 - x * x);
}

void test(rrfunction f, double exact)
{
  int m, N;
  double h, IhN;

  /* |t|≦10 まで計算することにする */
  h = 1.0; N = 10;
  /* h を半分, N を倍にして double exponential formula で計算してゆく */
  for (m = 1; m <= 10; m++) {
    IhN = de(f, h, N);
    printf("h=%f, I_h=%25.16f, I_h-I=%e\n", h, IhN, IhN - exact);
    h /= 2; N *= 2;
  }
}

int main(int argc, char **argv)
{
  if (argc >= 2 && strcmp(argv[1], "-d") == 0)
      debug = 1;
  pi = 4.0 * atan(1.0); halfpi = pi / 2;

  printf("test1 (sqrt(1-x^2) の積分)\n");
  test(f1, halfpi);

  printf("test2 (1/sqrt(1-x^2) の積分)\n");
  test(f2, pi);

  return 0;
}