R에서 함수곡선 아래의 면적을 구하려면 어떻게 할까? 예를 들어 y=sin(x) 함수 곡선에서 \(\frac{\pi}{4}\)부터 \(\frac{3}{4}\times\pi\) 까지의 면적은 얼마일까?

이 면적을 수동으로 계산하지 않더라도 R의 내장함수인 integrate()를 사용하면 쉽게 구할 수 있다.

integrate(sin,pi/4,pi*3/4)
## 1.414214 with absolute error < 1.6e-14

이번 장에서는 integrate()함수를 쓰지 않고 Newton-Cotes Formula를 이용하여 면적을 구하는 방법을 다룬다. 이 글의 전반부는 Hadley Wickam의 Advanced-R에 나온 내용에서 인용하였다.

10.5 Case study: numerical integration

Calculation with midpoint and trapezoid

값을 구하고자 하는 범위의 시작을 a, 끝을 b라고 하자.

  1. 먼저 midpoint는 a와 b의 중간값인 \(\frac{a+b}{2}\)에 대한 함수 값을 높이로 하고 b-a를 밑변으로 하는 사각형의 면적으로 구하는 방법이다.
midpoint <- function(f, a, b) {
  (b - a) * f((a + b) / 2)
}

(result=midpoint(sin,pi/4,pi*3/4))
[1] 1.570796

  1. trapezoid는 a의 함수값 f(a)와 b의 함수값 f(b)의 평균을 높이로 하는 사각형의 높이로 구하는 방법이다.
trapezoid1<-function(f,a,b){
    (b-a)*(f(a)+f(b))/2
}

(result=trapezoid1(sin,pi/4,pi*3/4))
## [1] 1.110721

  1. 보다 정밀한 방법으로 a부터 b까지를 몇 개의 구간으로 나누고 각 구간에 해당하는 사각형의 면적을 합해서 근사치를 만드는 방법을 생각해 볼 수 있다. 여러 가지 방법들이 알려져 있는데 이러한 방법들을 Issac Newton과 Roger Cotes의 이름을 따서 Newton-Cotes rules라고 이야기 하며 널리 알려진 formula는 다음과 같은 것들이 있다.

Newton-Cotes Rules

(https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas)

Simpson의 방법은 전체구간을 3개의 사각형으로 나누어 계산하는 방법이다.

Boole의 방법은 전체를 90개의 구간으로 나누고 5개의 사각형으로 면적을 측정한다.

traepzoid방법도 전체를 10개의 구간으로 나누어 적용한다면 보다 정밀해진다.

trapezoid_composite<-function(f,a,b,n=10){
    
    points <- seq(a,b,length=n+1)
    width <- (b-a)/n
    
    area<-0
    for(i in seq_len(n)){
        area<-area+width*f((points[i]+points[i+1])/2)
    }
    area
}

(result=trapezoid_composite(sin,0,pi/2))
## [1] 1.001029

Function for general Newton-Cotes Rule:

Hadley가 쓴 advanced-R 책에는 이와 같은 Newton-Cotes Rule을 일반화하는 함수를 소개하고 있다.

composite=function(f,a,b,n=10,rule){
    points<-seq(a,b,length=n+1)
    area<-0
    result=list()
    for(i in seq_len(n)){
        result[[i]]<-rule(f,points[i],points[i+1])
        area<-area+rule(f,points[i],points[i+1])
    }
    area
}

newton_cotes<-function(coef,open=FALSE){
    n<-length(coef)+open
    
    result=function(f,a,b){
        pos<- function(i) a+i*(b-a)/n
        points<-pos(seq.int(0,length(coef)-1))
        
        (b-a)/sum(coef)*sum(f(points)*coef)
        
    }
}

trapezoid2<-newton_cotes(c(1,1))

Recalcuation with generalized function

이 함수를 이용하여 0부터 pi까지의 값을 계산해보면 다음과 같고 이 값은 trapezoid1(sin,0,pi)의 값은 같아야 한다.

# It should be equal to trapezoid(sin,pi/4,pi*3/4)
composite(sin,0,pi,n=1,rule=trapezoid2)
[1] 1.570796
trapezoid1(sin,0,pi)
[1] 1.923671e-16

왜 값이 다를까?

값이 왜 다른지, 어디에서 에러가 나는지 보기위해 위키의 내용을 다시 리뷰해보았다.

(https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas)

Newton-Cote rule에서 자유도 n 은 계수의 갯수-1로 계산해야 한다. 하지만 Hadley도 신이 아니라 사람인지라 실수한 것으로 생각된다. 함수에 다음과 같이 고치니 바라는 대로 동작하였다.

newton_cotes1<-function(coef,open=FALSE){
    n<-length(coef)-1+open
    
    result=function(f,a,b){
        pos<- function(i) a+i*(b-a)/n
        points<-pos(seq.int(0,length(coef)-1))
        
        (b-a)/sum(coef)*sum(f(points)*coef)
        
    }
}
trapezoid3<-newton_cotes1(c(1,1))
composite(sin,0,pi,n=1,rule=trapezoid3)
[1] 1.923671e-16

Plot for the numerical integration

이 함수을 이용하여 plot을 그리기 위해 함수의 내용을 약간 변경하여 다음과 같은 함수를 제작하였다.

newton_cotes<-function(coef,open=FALSE){
    n<-length(coef)-1+open
    
    result=function(f,a,b){
        pos<- function(i) a+i*(b-a)/n
        points<-pos(seq.int(0,length(coef)-1))
        
        area=(b-a)/sum(coef)*sum(f(points)*coef)
        list(a=a,b=b,y=f(points),coef=coef,area=area)
    }
    result
}

composite=function(f,a,b,n=10,rule){
    points<-seq(a,b,length=n+1)
    area<-0
    result=list()
    for(i in seq_len(n)){
        result[[i]]<-rule(f,points[i],points[i+1])
        area<-area+rule(f,points[i],points[i+1])$area
    }
    list(result,f=f,area=area)
}
require(ggplot2)
plot_quadrature=function(result,color="red",fill="blue",alpha=0.2,size=0.2){
    test=result$f
    cadd <- function(x) Reduce("+", x, accumulate = TRUE)
    p<-ggplot(data.frame(x=c(0,pi)),aes(x))+stat_function(fun=test)
 
    for(i in 1:length(result[[1]])){
        subdata=result[[1]][[i]]
        coef=subdata$coef
        x=seq(subdata$a,subdata$b,length=(sum(coef)+1))
        h=subdata$y
        start=c(1,cadd(coef)+1)
        for(j in 1:length(coef)){
            p<-p+geom_rect(xmin=x[start[j]],xmax=x[start[j+1]],
                           ymin=0,ymax=h[j],colour=color,fill=fill,alpha=alpha,size=size)
        }
    }
    p
}

Plots for general Newton-Cotes Rule

이 함수를 이용하여 0부터 pi까지 사인함수 곡선의 면적을 측정하는 plot을 그리면 다음과 같다.

boole<-newton_cotes(c(7,32,12,32,7))
milne<-newton_cotes(c(2,-1,2),open=TRUE)
trapezoid<-newton_cotes(c(1,1))
simpson<-newton_cotes(c(1,4,1))
simpson8<-newton_cotes(c(1,3,3,1))

rules=list(trapezoid,simpson,simpson8,boole)

result=lapply(rules,function(x) composite(sin,0,pi,n=10,rule=x))
plot_quadrature(result[[1]])+ggtitle("Rule : Trapezoid, n=10")+
    geom_text(x=pi/2,y=sin(pi/2)/2,label=result[[1]]$area)

plot_quadrature(result[[2]])+ggtitle("Rule : Simpson's , n=10")+
    geom_text(x=pi/2,y=sin(pi/2)/2,label=result[[2]]$area)

plot_quadrature(result[[3]])+ggtitle("Rule : Simpson's 3/8, n=10")+
    geom_text(x=pi/2,y=sin(pi/2)/2,label=result[[3]]$area)

plot_quadrature(result[[4]])+ggtitle("Rule : Boole's, n=10")+
    geom_text(x=pi/2,y=sin(pi/2)/2,label=result[[4]]$area)

이번에는 0부터 pi까지를 100개의 구간으로 나누어 trapezoid rule을 적용해 보았다.

result<-composite(sin,0,pi,n=100,rule=trapezoid)
plot_quadrature(result,color=NA)+ggtitle("Rule : Trapezoid, n=100")+
    geom_text(x=pi/2,y=sin(pi/2)/2,label=result$area)

0부터 pi까지를 100개의 구간으로 나누어 Boole’s rule을 적용해 보니 참값인 2를 얻을 수 있었다.

result1<-composite(sin,0,pi,n=100,rule=boole)
plot_quadrature(result,color=NA)+ggtitle("Rule : Boole's, n=100")+
    geom_text(x=pi/2,y=sin(pi/2)/2,label=result1$area)

결론

Newton-Cotes Rule의 일반화를 위한 함수를 함수형 프로그래밍을 이용해 제작하였고 Hadley Wickham도 실수할 수도 있다는 것을 깨달았다. 이와 같이 책의 error를 발견한 경우 이 책의 github의 issue에 글을 올려 공유하기를 바란다. 오늘 언급한 내용도 advanced-R책의 github에 올려져 있다(https://github.com/hadley/adv-r/issues/783). 이 글을 읽으시는 분들께서도 책을 읽으시면서 함수를 눈으로 보는 것으로 만족하지 말고 직접 제작해보면서 공부할 경우 보다 많은 것을 얻을 수 있고 집단 지성을 발휘해서 책의 내용을 가다듬는데 일조할 수 있다.