수치 적분 신소재공학과 김 혜 민
과제 X=0~0.8 의 구간에서 f(x) 적분하기 Midpoint Rule Trapezoidal Rule Simpson’s Rule Romberg Integration
적분공식 Midpoint Trapezoidal Simpson’s Romberg
Midpoint Rule n=1; while(fabs(sum-A)>EL){ h=(b-a)/n; x=a+h/2; sum=0; for(i=1;i<=n;i++){ sum=sum+f(x); x=x+h; } sum=sum*h; n++; } 적분공식 A= , EL=10 -7
Trapezoidal Rule n=1; while(fabs(sum-A)>EL){ h=(b-a)/n; x=a; sum=f(a); for(i=1;i<=n-1;i++){ x=x+h; sum=sum+2*f(x); } sum=sum+f(b); sum=sum*(b-a)/(2*n); n++; } 적분공식
Simpson’s Rule n=2; while(fabs(sum-A)>EL){ h=(b-a)/n; x=a; sum=f(a); for(i=1;i<=(n/2-1);i++){ x=x+2*h; sum=sum+2*f(x); } x=a+h; for(i=1;i<=(n/2);i++){ sum=sum+4*f(x); x=x+2*h; } sum=sum+f(b); sum=sum*(b-a)/(3*n); n=n+2; } 적분공식 구간의 개수가 짝수일 경우만 계산
Romberg Integration romb[1][1]=(b-a)/2*(f(a)+f(b)); while(1){ sigma=0; for(i=1;i<=pow(2,k-2);i++){ sigma=sigma+f(a+(2*i-1)*(b-a)/pow(2,k-1)); } romb[k][1]=(romb[k-1][1]+(b-a)/pow(2,k-2)*sigma)/2; for(j=2;j<=k;j++){ romb[k][j]=romb[k][j-1]+(romb[k][j-1]-romb[k-1][j-1])/(pow(4,j-1)-1); } if(fabs(romb[k][k]-romb[k][k-1])<=EL) break; k++; } 적분공식 Extrapolation 을 해도 차이가 없을 때까지
실행결과 효율순위 : Simpson’s Rule > Midpoint Rule > Trapezoidal Rule
오차비교
실행결과 (Romberg Integration) 정답 : 와 거의 일치