Integrasi

# FUNGSI BUATAN
trapezoidal <- function(f, a, b, n) {
  # f : fungsi yang akan diintegralkan
  # a : batas bawah
  # b : batas atas
  # n : banyaknya partisi (subinterval)
  
  h <- (b - a) / n                 # panjang tiap subinterval
  x <- seq(a, b, by = h)           # titik-titik partisi
  y <- f(x)                        # nilai fungsi di titik-titik partisi
  
  integral <- (h/2) * (y[1] + 2*sum(y[2:n]) + y[n+1])
  return(integral)
}

# Ilustrasi
f <- function(x) x^3
a <- 0
b <- 1
n <- 1000

hasil <- trapezoidal(f, a, b, n)
print(hasil)

# FUNGSI BAWAAN
library(pracma)
x <- seq(0, 1, length.out = 1000)
y <- x^3
trapz(x, y)

# Bandingkan dengan nilai sejati
print(integrate(f, lower = 0, upper = 1)$value)

# FUNGSI BUATAN
midpoint_integral <- function(f, a, b, n) {
  h <- (b - a) / n
  sum(f(a + ( (1:n) - 0.5 ) * h)) * h
}

# Ilustrasi
f <- function(x) x^3
midpoint_integral(f, 0, 1, 1000)

# FUNGSI BUATAN
simpson13 <- function(f, a, b, n) {
  if (n %% 2 != 0) {
    stop("Jumlah selang (n) harus genap untuk metode Simpson 1/3")
  }
  
  h <- (b - a) / n
  x <- seq(a, b, by = h)
  y <- f(x)
  
  # rumus Simpson 1/3
  hasil <- (h/3) * (y[1] + 
                     4 * sum(y[seq(2, n, by = 2)]) + 
                     2 * sum(y[seq(3, n-1, by = 2)]) + 
                     y[n+1])
  return(hasil)
}

# Ilustrasi
f <- function(x) x^3
simpson13(f, 0, 1, 4)

# Bandingkan dengan nilai sejati
print(integrate(f, lower = 0, upper = 1)$value)

# FUNGSI BUATAN
simpson38 <- function(f, a, b, n) {
  if (n %% 3 != 0) {
    stop("Jumlah selang (n) harus kelipatan 3 untuk metode Simpson 3/8")
  }
  h <- (b - a) / n
  x <- seq(a, b, by = h)
  y <- f(x)
  # rumus Simpson 3/8
  hasil <- (3*h/8) * (y[1] + 
                      3 * sum(y[seq(2, n, by = 3)]) + 
                      3 * sum(y[seq(3, n, by = 3)]) + 
                      2 * sum(y[seq(4, n-2, by = 3)]) +
                      y[n+1])
  return(hasil)
}

# Ilustrasi
f <- function(x) x^3
simpson38(f, 0, 1, 6)

# FUNGSI BUATAN
mromberg <- function(f, a, b, nmax = 10, tol = 1e-8) {
  R <- matrix(0, nrow = nmax, ncol = nmax)
  h <- b - a
  R[1,1] <- 0.5 * h * (f(a) + f(b))
  
  for (i in 2:nmax) {
    h <- h / 2
    x <- seq(a + h, b - h, by = 2*h)
    R[i,1] <- 0.5 * R[i-1,1] + h * sum(f(x))
    for (k in 2:i) {
      R[i,k] <- (4^(k-1) * R[i,k-1] - R[i-1,k-1]) / (4^(k-1) - 1)
    }
    if (abs(R[i,i] - R[i-1,i-1]) < tol) {
      return(list(value = R[i,i], table = R[1:i,1:i]))
    }
  }
  return(list(value = R[nmax,nmax], table = R))
}

# ilustrasi
f <- function(x) x^3
mromberg(f, 0, 1)

# FUNGSI BAWAAN
library(pracma)
romberg(f, 0, 1)

#FUNGSI BAWAAN
library(statmod)

gauss_quad <- function(f, a, b, n = 5) {
  # f  : fungsi yang akan diintegralkan
  # a,b: batas bawah dan atas
  # n  : jumlah titik kuadratur
  
  # titik dan bobot Gauss-Legendre pada interval [-1,1]
  quad <- gauss.quad(n, kind = "legendre")
  x <- quad$nodes
  w <- quad$weights
  
  # transformasi ke interval [a,b]
  t <- ((b - a) * x + (b + a)) / 2
  integral <- ((b - a) / 2) * sum(w * f(t))
  
  return(integral)
}

# ilustrasi
f <- function(x) x^3
gauss_quad(f, a=0, b=1)