Numerical integration: Difference between revisions

Content added Content deleted

Inline

@@ Line 4,547: / Line 4,547: @@
 {{works with|R|2.11.0}}
-These presume that f can take a vector argument.
+The integ function defined below uses arbitrary abscissae and weights passed as argument (resp. u and v). It assumes that f can take a vector argument.
-<lang R>integrate.rect <- function(f, a, b, n, k=0) {
+<lang rsplus>integ <- function(f, a, b, n, u, v) {
+  h <- (b - a) / n
-  #k = 0 for left, 1 for right, 0.5 for midpoint
-  h <- (b-a)/n
+  s <- 0
-  x <- seq(a, b, len=n+1)
+  for (i in seq(0, n - 1)) {
-  sum(f(x[-1]-h*(1-k)))*h
+    s <- s + sum(v * f(a + i * h + u * h))
+  }
+  s * h
 }
-integrate.trapezoid <- function(f, a, b, n) {
+test <- function(f, a, b, n) {
+  c(rect.left = integ(f, a, b, n, 0, 1),
-  h <- (b-a)/n
-  x <- seq(a, b, len=n+1)
+    rect.right = integ(f, a, b, n, 1, 1),
+    rect.mid = integ(f, a, b, n, 0.5, 1),
-  fx <- f(x)
+    trapezoidal = integ(f, a, b, n, c(0, 1), c(0.5, 0.5)),
-  sum(fx[-1] + fx[-length(x)])*h/2
+    simpson = integ(f, a, b, n, c(0, 0.5, 1), c(1, 4, 1) / 6))
 }
+test(\(x) x^3, 0, 1, 100)
-integrate.simpsons <- function(f, a, b, n) {
+#  rect.left  rect.right    rect.mid trapezoidal     simpson
-  h <- (b-a)/n
+#  0.2450250   0.2550250   0.2499875   0.2500250   0.2500000
-  x <- seq(a, b, len=n+1)
-  fx <- f(x)
-  sum(fx[-length(x)] + 4*f(x[-1]-h/2) + fx[-1]) * h/6
-}
-f1 <- (function(x) {x^3})
-f2 <- (function(x) {1/x})
-f3 <- (function(x) {x})
-f4 <- (function(x) {x})
+test(\(x) 1 / x, 1, 100, 1000)
-integrate.simpsons(f1,0,1,100) #0.25
+#  rect.left  rect.right    rect.mid trapezoidal     simpson
-integrate.simpsons(f2,1,100,1000) # 4.60517
+#   4.654991    4.556981    4.604763    4.605986    4.605170
-integrate.simpsons(f3,0,5000,5000000) # 12500000
-integrate.simpsons(f4,0,6000,6000000) # 1.8e+07
+test(\(x) x, 0, 5000, 5e6)
-integrate.rect(f1,0,1,100,0) #TopLeft 0.245025
+#  rect.left  rect.right    rect.mid trapezoidal     simpson
-integrate.rect(f1,0,1,100,0.5) #Mid 0.2499875
+#   12499998    12500003    12500000    12500000    12500000
-integrate.rect(f1,0,1,100,1) #TopRight 0.255025
+test(\(x) x, 0, 6000, 6e6)
-integrate.trapezoid(f1,0,1,100) # 0.250025</lang>
+#  rect.left  rect.right    rect.mid trapezoidal     simpson
+#    1.8e+07     1.8e+07     1.8e+07     1.8e+07     1.8e+07</lang>
 =={{header|Racket}}==