Simpson's rule

[정현성]

[펌] 출처 영남대

5.2 Simpson의 법칙

[원 리]

사다리꼴 적분공식에서는 함수 곡선 위의 점들을 직선으로 연결하여 그 면적을 구하는데

비하여, Simpson의 법칙은 함수 곡선 위의 점들을 2차 곡선으로 근사하여 그 면적을 계산한다.

위 그림의 (a)와 같이 2차 곡선으로 근사하여 그 면적을 계산할 때, 계산을 쉽게 하기

위하여 그림 (b)와 같이 곡선을 평행 이동시켜 2차 곡선식으로부터 면적을 계산한다.

		-- (1)
		-- (2)

식 (2)를 이용하여 적분값을 계산하기 위하여 식 (1)의 계수들을 결정하여야 한다. 식 (1)의

함수는 그림 (b)에서 보는 바와 같이 세 점을 지나므로, 식 (3)과 같이 쓸 수 있다.

-- (3)

식 (3)의 연립방정식에서 식 (4)와 같이 계수들을 구할 수 있다.

-- (4)

[입 력]

(1)	: 적분값을 구할 방정식
(2) 적분 구간	: 적분 구간
(3)	: 적분 구간의 분할 수

[방 법]

1. 분할 간격 를 계산한다.
2. 각 구간에 대하여 식 (4)와 같이 계수들을 계산하고, 식(2)와 같이 면적을 구하여 더한다.

 

/***********************************************************/

/*                                                         */

/*          Chap.5  Numerical Integration                  */

/*             5.2  Simpson's Rule                         */

/*                                                         */

/***********************************************************/

 

#include < stdio.h >

#include < stdlib.h >

#include < math.h >

 

#define   F(x)          (sin(x)+cos(x))     //..Function

#define   A             -3.0                //..Integration start-pos

#define   B              2.0                //..Integration end  -pos

 

void main()

{

    int         i,N;

    double      delX,FI;

    double              y1,y2,y3;

    double              coeffA,coeffB,coeffC;

 

    printf(" number of gap : ");

    scanf("%d",&N);

 

    delX=(B-A)/(double)N;

 

    FI=0.;

    for(i=0;i < N;i++){

        y1=F(A+i*delX);

        y2=F(A+i*delX+delX/2.);

        y3=F(A+(i+1)*delX);

        coeffA=2./pow(delX,2.)*(y1-2.*y2+y3);

        coeffB=1./(delX)*(y3-y1);

        coeffC=y2;

 

        FI+=coeffA/12.*pow(delX,3.)+coeffC*delX;

    }

 

    printf("\nI is %f\n",FI);

}

 

C..........................................................

C           Chap.5  Numerical Integration 

C              5.2  Simpson's Rule

C..........................................................

C

C     ..Function

      F(X)=SIN(X)+COS(X)

      A=-3.0

      B= 2.0

      WRITE(*,*) ' number of gap : '

      READ(*,*) N

      DELX=(B-A)/N

      FI=0.

      DO 100 I=0,N-1

         Y1=F(A+I*DELX)

         Y2=F(A+I*DELX+DELX/2.)

         Y3=F(A+(I+1)*DELX)

         COEFFA=2./DELX**2*(Y1-2.*Y2+Y3)

         COEFFB=1./DELX*(Y3-Y1)

         COEFFC=Y2

         FI=FI+COEFFA/12.*DELX**3+COEFFC*DELX

  100 CONTINUE

      WRITE(*,110) FI

  110 FORMAT('FI is ',E15.7)

      STOP

      END