RosettaCodeData/Task/Numerical-integration/R/numerical-integration.r

33 lines
1.0 KiB
R

integ <- function(f, a, b, n, u, v) {
h <- (b - a) / n
s <- 0
for (i in seq(0, n - 1)) {
s <- s + sum(v * f(a + i * h + u * h))
}
s * h
}
test <- function(f, a, b, n) {
c(rect.left = integ(f, a, b, n, 0, 1),
rect.right = integ(f, a, b, n, 1, 1),
rect.mid = integ(f, a, b, n, 0.5, 1),
trapezoidal = integ(f, a, b, n, c(0, 1), c(0.5, 0.5)),
simpson = integ(f, a, b, n, c(0, 0.5, 1), c(1, 4, 1) / 6))
}
test(\(x) x^3, 0, 1, 100)
# rect.left rect.right rect.mid trapezoidal simpson
# 0.2450250 0.2550250 0.2499875 0.2500250 0.2500000
test(\(x) 1 / x, 1, 100, 1000)
# rect.left rect.right rect.mid trapezoidal simpson
# 4.654991 4.556981 4.604763 4.605986 4.605170
test(\(x) x, 0, 5000, 5e6)
# rect.left rect.right rect.mid trapezoidal simpson
# 12499998 12500003 12500000 12500000 12500000
test(\(x) x, 0, 6000, 6e6)
# rect.left rect.right rect.mid trapezoidal simpson
# 1.8e+07 1.8e+07 1.8e+07 1.8e+07 1.8e+07