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

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

/*
例題 4-8 利用 辛普森理則 計算 雙重積分 f(x,y)= 4xy 在[0 , 1-x^2] dy 與 [0,1] dx 的定積分
# Based on Simpson's  Rule to  compute the double integral.
*/
#include <stdio.h>
#include <math.h>
#define F(x,y)  (4*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-x*x;
}

void main()
{

    int i,j,m,n;
    double a,b,hx,sum1=0.0,sum2=0.0,sn;
    a=0;
    b=1;
    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.1000 1.0000 0.0000 10 -->x[0]=0.0000 , hy[0]=0.1000
y[0,0]=0.0000
y[0,1]=0.1000
y[0,2]=0.2000
y[0,3]=0.3000
y[0,4]=0.4000
y[0,5]=0.5000
y[0,6]=0.6000
y[0,7]=0.7000
y[0,8]=0.8000
y[0,9]=0.9000
y[0,10]=1.0000
0.0990 0.9900 0.0000 10 -->x[1]=0.1000 , hy[1]=0.0990
y[1,0]=0.0000
y[1,1]=0.0990
y[1,2]=0.1980
y[1,3]=0.2970
y[1,4]=0.3960
y[1,5]=0.4950
y[1,6]=0.5940
y[1,7]=0.6930
y[1,8]=0.7920
y[1,9]=0.8910
y[1,10]=0.9900
0.0960 0.9600 0.0000 10 -->x[2]=0.2000 , hy[2]=0.0960
y[2,0]=0.0000
y[2,1]=0.0960
y[2,2]=0.1920
y[2,3]=0.2880
y[2,4]=0.3840
y[2,5]=0.4800
y[2,6]=0.5760
y[2,7]=0.6720
y[2,8]=0.7680
y[2,9]=0.8640
y[2,10]=0.9600
0.0910 0.9100 0.0000 10 -->x[3]=0.3000 , hy[3]=0.0910
y[3,0]=0.0000
y[3,1]=0.0910
y[3,2]=0.1820
y[3,3]=0.2730
y[3,4]=0.3640
y[3,5]=0.4550
y[3,6]=0.5460
y[3,7]=0.6370
y[3,8]=0.7280
y[3,9]=0.8190
y[3,10]=0.9100
0.0840 0.8400 0.0000 10 -->x[4]=0.4000 , hy[4]=0.0840
y[4,0]=0.0000
y[4,1]=0.0840
y[4,2]=0.1680
y[4,3]=0.2520
y[4,4]=0.3360
y[4,5]=0.4200
y[4,6]=0.5040
y[4,7]=0.5880
y[4,8]=0.6720
y[4,9]=0.7560
y[4,10]=0.8400
0.0750 0.7500 0.0000 10 -->x[5]=0.5000 , hy[5]=0.0750
y[5,0]=0.0000
y[5,1]=0.0750
y[5,2]=0.1500
y[5,3]=0.2250
y[5,4]=0.3000
y[5,5]=0.3750
y[5,6]=0.4500
y[5,7]=0.5250
y[5,8]=0.6000
y[5,9]=0.6750
y[5,10]=0.7500
0.0640 0.6400 0.0000 10 -->x[6]=0.6000 , hy[6]=0.0640
y[6,0]=0.0000
y[6,1]=0.0640
y[6,2]=0.1280
y[6,3]=0.1920
y[6,4]=0.2560
y[6,5]=0.3200
y[6,6]=0.3840
y[6,7]=0.4480
y[6,8]=0.5120
y[6,9]=0.5760
y[6,10]=0.6400
0.0510 0.5100 0.0000 10 -->x[7]=0.7000 , hy[7]=0.0510
y[7,0]=0.0000
y[7,1]=0.0510
y[7,2]=0.1020
y[7,3]=0.1530
y[7,4]=0.2040
y[7,5]=0.2550
y[7,6]=0.3060
y[7,7]=0.3570
y[7,8]=0.4080
y[7,9]=0.4590
y[7,10]=0.5100
0.0360 0.3600 0.0000 10 -->x[8]=0.8000 , hy[8]=0.0360
y[8,0]=0.0000
y[8,1]=0.0360
y[8,2]=0.0720
y[8,3]=0.1080
y[8,4]=0.1440
y[8,5]=0.1800
y[8,6]=0.2160
y[8,7]=0.2520
y[8,8]=0.2880
y[8,9]=0.3240
y[8,10]=0.3600
0.0190 0.1900 0.0000 10 -->x[9]=0.9000 , hy[9]=0.0190
y[9,0]=0.0000
y[9,1]=0.0190
y[9,2]=0.0380
y[9,3]=0.0570
y[9,4]=0.0760
y[9,5]=0.0950
y[9,6]=0.1140
y[9,7]=0.1330
y[9,8]=0.1520
y[9,9]=0.1710
y[9,10]=0.1900
0.0000 0.0000 0.0000 10 -->x[10]=1.0000 , hy[10]=0.0000
y[10,0]=0.0000
y[10,1]=0.0000
y[10,2]=0.0000
y[10,3]=0.0000
y[10,4]=0.0000
y[10,5]=0.0000
y[10,6]=0.0000
y[10,7]=0.0000
y[10,8]=0.0000
y[10,9]=0.0000
y[10,10]=0.0000
0.0000 0.0000 0.0000 0.0000 gg[0]=0.0000
0.0000 0.3960 1.5840 3.9600 gg[1]=0.1960
0.0000 0.7680 3.0720 7.6800 gg[2]=0.3686
0.0000 1.0920 4.3680 10.9200 gg[3]=0.4969
0.0000 1.3440 5.3760 13.4400 gg[4]=0.5645
0.0000 1.5000 6.0000 15.0000 gg[5]=0.5625
0.0000 1.5360 6.1440 15.3600 gg[6]=0.4915
0.0000 1.4280 5.7120 14.2800 gg[7]=0.3641
0.0000 1.1520 4.6080 11.5200 gg[8]=0.2074
0.0000 0.6840 2.7360 6.8400 gg[9]=0.0650
0.0000 0.0000 0.0000 0.0000 gg[10]=0.0000
S10=0.333400

Command exited with non-zero status 13