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

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

/* 習題 8-1.c based on Trapezoidal Rule to
 * compute the double integral.
 */
#include <stdio.h>
#include <math.h>
#define F(x,y)  (x+y)
#define C(x)    (pow(x,2))
#define D(x)    (x+2)
void gy(int);
double x[50],y[50][50],gg[50],hy[50];
void main()
{
   int i,j,n;
   double a,b,hx,sum=0.0,ts;
   n=10;
   a=-1.0;
   b=2.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=9.315405