C語言 例題 4-8-(2) 利用 辛普森理則 計算 雙重積分 f(x,y)= (sin(x+y)) 在[ln(x) , 3+exp(x/5)] dy 與 [1,3] dx 的定積分

C語言 例題 4-8-(2) 利用 辛普森理則 計算 雙重積分 f(x,y)= (sin(x+y))  在[ln(x) , 3+exp(x/5)] dy 與 [1,3] dx 的定積分

/*
例題 4-8-(2) 利用 辛普森理則 計算 雙重積分 f(x,y)= (sin(x+y))
在[ln(x) , 3+exp(x/5)] dy 與 [1,3] dx 的定積分
 c=ln(x) , d=3+exp(x/5)  , a=1  , b=3
 C 語言 math.h 的函數log() 回傳參數以常數e 為底數所計算出的對數值。
# Based on Simpson's  Rule to  compute the double integral.
*/
#include <stdio.h>
#include <math.h>
#define F(x,y)  sin(x+y)

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

double C(double x)
{
    return log(x);
}

double D(double x)
{
    return 3+exp(x/5);
}

void main()
{

    int i,j,m,n;
    double a,b,hx,sum1=0.0,sum2=0.0,sn;
    a=1;
    b=3;
    n=10;
   
    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;
}
//==========================================