「advanced R」10.5数值积分

数值积分的思想很简单：通过用简单的成分金斯曲线来寻找曲线下面积。两个最简单的做法便是中点法则和梯形法则，中点法则利用的是矩阵近似曲线，梯形法则利用的是梯形近似。每个法则都要以要积分的函数f和积分范围（a_{b）为参数进行积分。例如，尝试对sinx(0}pi)积分。

1. 中点法则：

midpoint <- function(f,a,b){
  
  h=f((a+b)/2)
  s=h*(b-a)
  return(s)
  
}


midpoint(sin,0,2*pi)
trapzoid(sin,0,2*pi)

2.梯形法则

trapzoid <- function(f,a,b){
  
  h=b-a
  s=(f(a)+f(b))/2*h
  return(s)
  
}

为了使答案更加准确，我们可以将整个积分区间划分为多个小区间，然后对每个小区间进行积分，这称为组合积分，使用下面两个新函数来实现它

# 区间划分 --------------------------------------------------------------------

midpoint_com <- function(f,a,b,n=10){
  points = seq(a,b,length.out = n+1)
  h = (b-a)/n
  sum=0
  for (i in seq_len(n)) {
    area=f((points[i]+points[i+1])/2)*h
    sum=sum+area
  }
  
  return(sum)
  
}


trapezoid_com <- function(f,a,b,n=10){
  points=seq(a,b,length.out = n+1)
  h=(b-a)/n
  sum=0
  for (i in seq_len(n)) {
    area=(f(points[i])+f(points[i+1]))*h/2
    sum=area+sum
  }
  return(sum)
}

trapezoid(sin,0,pi,100) 
midpoint_com(sin,0,pi,100)
midpoint(sin,0,pi,100)

composite <- function(f,a,b,n=10,rule){
  points=seq(a,b,length.out = n+1)
  h=(b-a)/n
  sum=0
  for (i in seq_len(n)) {
    area=rule(f,points[i],points[i+1])
    sum=area+sum
  }
  return(sum)
}
}

composite(sin,0,pi,n=10,rule =midpoint )

composite(sin,0,pi,n=10,rule =trapzoid )