２重積分の計算をシンプソン法を用いて計算する（C++）

２重積分の計算をシンプソン法を用いて計算する（１次元をシンプソン法はこちら）。
※参考にしたページhttp://www5.airnet.ne.jp/tomy/cpro/science.htm

C++ソース

/*
２次元等間隔シンプソン法
例：f(x,y) = sin(x)*sin(y) [0:PI]
【参考】http://www5.airnet.ne.jp/tomy/cpro/science.htm
*/
#include <math.h>
#include <stdio.h>
double PI = acos(-1.0);

double F(double x, double y);
double simpe2(double **f, const int m, const int n, double h1, double h2);

int main(void){
  double h1, h2, x,y,s,exa ;
  int i, j;

  const int m = 1000; //x軸の刻み数
  const int n = 1000; //y軸の刻み数
  double x0 = 0.0 , x1 = PI ; //xの積分範囲
  double y0 = 0.0 , y1 = PI ; //yの積分範囲
  double **f= new double*[m+1];  // double型 m+1 個分の領域を動的確保
  for (i= 0; i<=m; ++i) f[i]= new double[n+1];  // double型 n+1 個分の領域を動的確保

  for(i = 0; i <= m; i++){
    for (j = 0; j <= n; j++)  {
      x = x0 + (x1-x0)/double(m) * double(i);
      y = y0 + (y1-y0)/double(n) * double(j);
      f[i][j] = F(x,y);
    }
  }
  h1 = (x1-x0) / double(m);
  h2 = (y1-y0) / double(n);
  exa = 4.0; // 解析値
  s = simpe2(f, m, n, h1, h2); // 計算

  for (i= 0; i<=m; ++i) delete[] f[i]; // 動的に確保した領域をそれぞれ解放
  delete[] f;

  //結果の書き出し
  printf("  刻み数 : %d × %d\n", m,n);
  printf("  積分値 : %22.15e\n", s);
  printf("  誤　差 : %22.15e\n", s - exa);

  return 1;
}
double F(double x, double y) //被積分関数
{
  return sin(x)*sin(y);
}

double simpe2(double **f, const int m, const int n, double h1, double h2)
{
  int i, j;
  double v;
  double *temp;
  temp = new double[m+1];  //double型 m+1 個分の領域を動的確保

  for(i = 0; i <= m; i++){
    v = - f[i][0] + f[i][n];
    for(j = 0 ; j < n - 1; j += 2)
      v += (2 * f[i][j] + 4 * f[i][j + 1]);
    temp[i] = v;
  }
  v = - temp[0] + temp[m];
  for(i = 0; i < m - 1; i += 2)  v += (2 * temp[i] + 4 * temp[i + 1]);
  delete [] temp;
  return v * h1 * h2 / 9.0;
}

実行結果