範例EX4-8 梯行積分法 雙重積分計算 --- C語言

範例EX4-8 梯行積分法 雙重積分計算 --- C語言

/* ex4-8.c based on Trapezoidal Rule to
 * compute the double integral.
 */
#include <stdio.h>
#include <math.h>
#define F(x,y)  (4*x*y)
#define C(x)    (0)
#define D(x)    (1-pow(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;
   scanf("n=%d a=%lf b=%lf",&n,&a,&b);
   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 0;
}

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;
}

輸入資料
n=10  a=0.0  b=1.0

輸出結果
T10=0.331650