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

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

// C program to calculate double integral value

/*
    Simpson's Rule to conpute f(x,y)=x*exp(y) , a=0 ,b=1
    c= 0 , d=x   double integral .
*/
#include <stdio.h>
#include <math.h>
#define F(x,y)  (x*exp(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 x;
}

void main()
{

    int i,j,m,n;
    double a,b,hx,sum1=0.0,sum2=0.0,sn;
    a=0;
    b=1.0;
    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.0000 0.0000 0.0000 2 -->x[0]=0.0000 , hy[0]=0.0000
y[0,0]=0.0000
y[0,1]=0.0000
y[0,2]=0.0000
0.2500 0.5000 0.0000 2 -->x[1]=0.5000 , hy[1]=0.2500
y[1,0]=0.0000
y[1,1]=0.2500
y[1,2]=0.5000
0.5000 1.0000 0.0000 2 -->x[2]=1.0000 , hy[2]=0.5000
y[2,0]=0.0000
y[2,1]=0.5000
y[2,2]=1.0000
0.0000 0.0000 0.0000 0.0000 gg[0]=0.0000
0.5000 0.8244 0.0000 2.5681 gg[1]=0.3244
1.0000 2.7183 0.0000 6.5949 gg[2]=1.7189
S2=0.502722