2.3 数値例

以下にあげる例は、サンプル・プログラム (C言語で記述) を用意してある。 現象数理学科 Mac であれば、ターミナルで以下のコマンドを実行すれば動くはず。

curl -O http://nalab.mind.meiji.ac.jp/~mk/complex2/prog20180625.tar.gz
tar xzf prog20180625.tar.gz
cd prog20180625
make

中点公式、台形公式、Simpson公式は次のコードで計算できる (以下の数値実験例のプログラムの多くで、 この nc.c をインクルードして使っている)。

nc.c

/*
 * nc.c --- Newton-Cotes の積分公式: 複合中点公式, 複合台形公式, 複合Simpson公式
 */

typedef double ddfunction(double);

double midpoint(ddfunction, double, double, int);
double trapezoidal(ddfunction, double, double, int);
double simpson(ddfunction, double, double, int);

/* 関数 f の [a,b] における積分の複合中点公式による数値積分 M_N */
double midpoint(ddfunction f, double a, double b, int N)
{
  int j;
  double h, M;
  h = (b - a) / N;
  M = 0.0;
  for (j = 1; j <= N; j++) M += f(a + (j - 0.5) * h);
  M *= h;
  return M;
}

/* 関数 f の [a,b] における積分の複合台形公式による数値積分 T_N */
double trapezoidal(ddfunction f, double a, double b, int N)
{
  int j;
  double h, T;
  h = (b - a)/ N;
  T = (f(a) + f(b)) / 2;
  for (j = 1; j < N; j++) T += f(a + j * h);
  T *= h;
  return T;
}

/* 関数 f の [a,b] における積分の複合 Simpson 公式による数値積分 S_{N} */
double simpson(ddfunction f, double a, double b, int N)
{
  int m = N / 2;
  return (trapezoidal(f, a, b, m) + 2 * midpoint(f, a, b, m)) / 3;
}