C語言 例題4-6 # 利用 辛普森 理則計算 雙重積分 f(x,y)= x^2 y 在[1,2] dx 與 [0,1] dy 的定積分

C語言  例題4-6 # 利用 辛普森 理則計算 雙重積分 f(x,y)= x^2 y 在[1,2] dx 與 [0,1] dy 的定積分

/*
    Simpson's Rule to conpute f(x,y)=x^2*y , a=0 ,b=1
    c= 1 , d=2  double integral .
*/
#include <stdio.h>
#include <math.h>
#define F(x,y)  (x*x*y)

void gy(int);
double x[50],y[50][50],gg[50],hy[50];

double C(double x)
{
    return 0.0;
}

double D(double x)
{
    return 1.0;
}

void main()
{

    int i,j,m,n;
    double a,b,hx,sum1=0.0,sum2=0.0,sn;
    a=1;
    b=2;
    n=2;
 
    m=n/2;
    hx=(b-a)/n;
    for(i=0;i<=n;i++)
    {
        x[i]=a+i*hx;
        hy[i]=( D(x[i] )-C( x[i] )) / n;
        printf("%1.4lf %1.4lf %1.4lf %1d -->", hy[i], D(x[i]), C(x[i]) ,n );
        printf("x[%1d]=%1.4lf , hy[%1d]=%1.4lf \n",i,x[i],i,hy[i]);
     
        for(j=0;j<=n;j++)
        {
            y[i][j]=C(x[i])+ j*hy[i];
            printf("y[%1d,%1d]=%1.4lf \n",i,j,y[i][j]);
        } 
    }
 
    gy(n);
 

    for(i=1;i<=2*m-1;i++)
    {
        if(i%2==0)
            sum2=sum2+gg[i];
        else
            sum1=sum1+gg[i];
    }
    sn=(hx/3.0)*(gg[0]+gg[n]+2.0*sum2+4.0*sum1);
    printf("S%d=%lf\n",n,sn);
    return;
}
//==========================================
void gy(int n)
{
    int i,j,m;
    double sum1=0.0 , sum2=0.0;
    m=n/2;
    for(i=0;i<=n;i++)
    {
     
    for(j=1;j<=2*m-1;j++)
        {
            if(j%2==0)
                sum2=sum2+ F(x[i],y[i][j]);
            else
                sum1=sum1+ F(x[i],y[i][j]);
        }

        gg[i]=(hy[i]/3.0)*(F(x[i],y[i][0])+F(x[i],y[i][n])+2.0*sum2+4.0*sum1);
        printf("%1.4lf %1.4lf %1.4lf %1.4lf \t", F(x[i],y[i][0]) , F(x[i],y[i][n]) ,2.0*sum2 ,4.0*sum1);
     
        sum1=0.0;
        sum2=0.0;
        printf("gg[%1d]=%1.4lf \n",i,gg[i]);
     
    }
   return;
}
//==========================================



輸出畫面
0.5000 1.0000 0.0000 2 -->x[0]=1.0000 , hy[0]=0.5000
y[0,0]=0.0000
y[0,1]=0.5000
y[0,2]=1.0000
0.5000 1.0000 0.0000 2 -->x[1]=1.5000 , hy[1]=0.5000
y[1,0]=0.0000
y[1,1]=0.5000
y[1,2]=1.0000
0.5000 1.0000 0.0000 2 -->x[2]=2.0000 , hy[2]=0.5000
y[2,0]=0.0000
y[2,1]=0.5000
y[2,2]=1.0000
0.0000 1.0000 0.0000 2.0000 gg[0]=0.5000
0.0000 2.2500 0.0000 4.5000 gg[1]=1.1250
0.0000 4.0000 0.0000 8.0000 gg[2]=2.0000
S2=1.166667

Command exited with non-zero status 12