習題 EX4-8-2 梯行法雙重積分---C語言

//習題 EX4-8-2 梯行法雙重積分

/* 習題 8-2.c based on Trapezoidal Rule to
 * compute the double integral.
 */
#include <stdio.h>
#include <math.h>
#define F(x,y)  ( sin (x+y))
#define C(x)    (log(x))
#define D(x)    (3+ exp(0.2*x))
void gy(int);
double x[50],y[50][50],gg[50],hy[50];
void main()
{
   int i,j,n;
   double a,b,hx=0.0,sum=0.0,ts=0.0;
   n=10;
   a=1.0;
   b=3.0;
 
   hx=(b-a)/n;
   for(i=0;i<=n;i++)
   {
      x[i]=a+i*hx;
      hy[i]=(D(x[i])-C(x[i]))/n;
      for(j=0;j<=n;j++)
         y[i][j]=C(x[i])+j*hy[i];
   }
   gy(n);
   for(i=1;i<=n-1;i++)
      sum=sum+gg[i];
   ts=(hx/2.0)*(gg[0]+gg[n]+2*sum);
   printf("T%d=%.6lf\n",n,ts);
   return;
}

void gy(int n)
{
   int i,j;
   double sum=0.0;
   for(i=0;i<=n;i++)
   {
      for(j=1;j<=n-1;j++)
      {
  sum=sum+F(x[i],y[i][j]);
      }
      gg[i]=(hy[i]/2.0)*(F(x[i],y[i][0])+
     F(x[i],y[i][n])+2*sum);
      sum=0.0;
   }
   return;
}

輸出結果
T10=-2.562917