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

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

/*
例題 4-8-(1) 利用 辛普森理則 計算 雙重積分 f(x,y)= x+y
在[x*x , x+2] dy 與 [-1, 2] dx 的定積分
 c=x*x , d=x+2  , a=-1  , b=2
# Based on Simpson's  Rule to  compute the double integral.
*/
#include <stdio.h>
#include <math.h>
#define F(x,y)  (x+y)

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

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

double D(double x)
{
    return x+2;
}

void main()
{

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



輸出畫面
0.0000 1.0000 1.0000 10 -->x[0]=-1.0000 , hy[0]=0.0000
y[0,0]=1.0000
y[0,1]=1.0000
y[0,2]=1.0000
y[0,3]=1.0000
y[0,4]=1.0000
y[0,5]=1.0000
y[0,6]=1.0000
y[0,7]=1.0000
y[0,8]=1.0000
y[0,9]=1.0000
y[0,10]=1.0000
0.0810 1.3000 0.4900 10 -->x[1]=-0.7000 , hy[1]=0.0810
y[1,0]=0.4900
y[1,1]=0.5710
y[1,2]=0.6520
y[1,3]=0.7330
y[1,4]=0.8140
y[1,5]=0.8950
y[1,6]=0.9760
y[1,7]=1.0570
y[1,8]=1.1380
y[1,9]=1.2190
y[1,10]=1.3000
0.1440 1.6000 0.1600 10 -->x[2]=-0.4000 , hy[2]=0.1440
y[2,0]=0.1600
y[2,1]=0.3040
y[2,2]=0.4480
y[2,3]=0.5920
y[2,4]=0.7360
y[2,5]=0.8800
y[2,6]=1.0240
y[2,7]=1.1680
y[2,8]=1.3120
y[2,9]=1.4560
y[2,10]=1.6000
0.1890 1.9000 0.0100 10 -->x[3]=-0.1000 , hy[3]=0.1890
y[3,0]=0.0100
y[3,1]=0.1990
y[3,2]=0.3880
y[3,3]=0.5770
y[3,4]=0.7660
y[3,5]=0.9550
y[3,6]=1.1440
y[3,7]=1.3330
y[3,8]=1.5220
y[3,9]=1.7110
y[3,10]=1.9000
0.2160 2.2000 0.0400 10 -->x[4]=0.2000 , hy[4]=0.2160
y[4,0]=0.0400
y[4,1]=0.2560
y[4,2]=0.4720
y[4,3]=0.6880
y[4,4]=0.9040
y[4,5]=1.1200
y[4,6]=1.3360
y[4,7]=1.5520
y[4,8]=1.7680
y[4,9]=1.9840
y[4,10]=2.2000
0.2250 2.5000 0.2500 10 -->x[5]=0.5000 , hy[5]=0.2250
y[5,0]=0.2500
y[5,1]=0.4750
y[5,2]=0.7000
y[5,3]=0.9250
y[5,4]=1.1500
y[5,5]=1.3750
y[5,6]=1.6000
y[5,7]=1.8250
y[5,8]=2.0500
y[5,9]=2.2750
y[5,10]=2.5000
0.2160 2.8000 0.6400 10 -->x[6]=0.8000 , hy[6]=0.2160
y[6,0]=0.6400
y[6,1]=0.8560
y[6,2]=1.0720
y[6,3]=1.2880
y[6,4]=1.5040
y[6,5]=1.7200
y[6,6]=1.9360
y[6,7]=2.1520
y[6,8]=2.3680
y[6,9]=2.5840
y[6,10]=2.8000
0.1890 3.1000 1.2100 10 -->x[7]=1.1000 , hy[7]=0.1890
y[7,0]=1.2100
y[7,1]=1.3990
y[7,2]=1.5880
y[7,3]=1.7770
y[7,4]=1.9660
y[7,5]=2.1550
y[7,6]=2.3440
y[7,7]=2.5330
y[7,8]=2.7220
y[7,9]=2.9110
y[7,10]=3.1000
0.1440 3.4000 1.9600 10 -->x[8]=1.4000 , hy[8]=0.1440
y[8,0]=1.9600
y[8,1]=2.1040
y[8,2]=2.2480
y[8,3]=2.3920
y[8,4]=2.5360
y[8,5]=2.6800
y[8,6]=2.8240
y[8,7]=2.9680
y[8,8]=3.1120
y[8,9]=3.2560
y[8,10]=3.4000
0.0810 3.7000 2.8900 10 -->x[9]=1.7000 , hy[9]=0.0810
y[9,0]=2.8900
y[9,1]=2.9710
y[9,2]=3.0520
y[9,3]=3.1330
y[9,4]=3.2140
y[9,5]=3.2950
y[9,6]=3.3760
y[9,7]=3.4570
y[9,8]=3.5380
y[9,9]=3.6190
y[9,10]=3.7000
0.0000 4.0000 4.0000 10 -->x[10]=2.0000 , hy[10]=0.0000
y[10,0]=4.0000
y[10,1]=4.0000
y[10,2]=4.0000
y[10,3]=4.0000
y[10,4]=4.0000
y[10,5]=4.0000
y[10,6]=4.0000
y[10,7]=4.0000
y[10,8]=4.0000
y[10,9]=4.0000
y[10,10]=4.0000
0.0000 0.0000 0.0000 0.0000 gg[0]=0.0000
-0.2100 0.6000 1.5600 3.9000 gg[1]=0.1579
-0.2400 1.2000 3.8400 9.6000 gg[2]=0.6912
-0.0900 1.8000 6.8400 17.1000 gg[3]=1.6159
0.2400 2.4000 10.5600 26.4000 gg[4]=2.8512
0.7500 3.0000 15.0000 37.5000 gg[5]=4.2188
1.4400 3.6000 20.1600 50.4000 gg[6]=5.4432
2.3100 4.2000 26.0400 65.1000 gg[7]=6.1520
3.3600 4.8000 32.6400 81.6000 gg[8]=5.8752
4.5900 5.4000 39.9600 99.9000 gg[9]=4.0460
6.0000 6.0000 48.0000 120.0000 gg[10]=0.0000
S10=9.448380

Command exited with non-zero status 13