Persiapan

knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  fig.align = "center", fig.width = 8, fig.height = 5, comment = "#>"
)
options(scipen = 999)

paket <- c("deSolve", "lpSolve", "GA")
belum_terpasang <- paket[
  !vapply(paket, requireNamespace, logical(1), quietly = TRUE)
]
if (length(belum_terpasang) > 0) {
  install.packages(belum_terpasang, repos = "https://cloud.r-project.org")
}
library(deSolve)
## Warning: package 'deSolve' was built under R version 4.5.3
library(lpSolve)
## Warning: package 'lpSolve' was built under R version 4.5.2
library(GA)
## Warning: package 'GA' was built under R version 4.5.3
## Loading required package: foreach
## Loading required package: iterators
## Package 'GA' version 3.2.5
## Type 'citation("GA")' for citing this R package in publications.
## 
## Attaching package: 'GA'
## The following object is masked from 'package:utils':
## 
##     de

Latihan 1 Akar Persamaan Nonlinear

Tentukan akar \(x^3-2x-5=0\) dengan bisection, Newton-Raphson, dan secant pada toleransi \(10^{-10}\), lalu verifikasi menggunakan uniroot().

f_akar <- function(x) x^3 - 2 * x - 5
df_akar <- function(x) 3 * x^2 - 2

bisection <- function(fun, a, b, tol = 1e-10, maxit = 1000) {
  fa <- fun(a); fb <- fun(b)
  if (abs(fa) <= tol) return(list(root = a, iter = 0L, converged = TRUE))
  if (abs(fb) <= tol) return(list(root = b, iter = 0L, converged = TRUE))
  if (fa * fb > 0) stop("f(a) dan f(b) harus memiliki tanda berbeda.")

  for (i in seq_len(maxit)) {
    tengah <- (a + b) / 2
    ft <- fun(tengah)
    if (abs(ft) <= tol || (b - a) / 2 <= tol) {
      return(list(root = tengah, iter = i, converged = TRUE))
    }
    if (fa * ft < 0) {
      b <- tengah; fb <- ft
    } else {
      a <- tengah; fa <- ft
    }
  }
  list(root = (a + b) / 2, iter = maxit, converged = FALSE)
}

newton_raphson <- function(fun, dfun, x0, tol = 1e-10, maxit = 1000) {
  x <- x0
  for (i in seq_len(maxit)) {
    fx <- fun(x); dfx <- dfun(x)
    if (!is.finite(dfx) || abs(dfx) < .Machine$double.eps) {
      stop("Turunan terlalu dekat dengan nol.")
    }
    x_baru <- x - fx / dfx
    if (abs(x_baru - x) <= tol || abs(fun(x_baru)) <= tol) {
      return(list(root = x_baru, iter = i, converged = TRUE))
    }
    x <- x_baru
  }
  list(root = x, iter = maxit, converged = FALSE)
}

secant <- function(fun, x0, x1, tol = 1e-10, maxit = 1000) {
  f0 <- fun(x0); f1 <- fun(x1)
  for (i in seq_len(maxit)) {
    penyebut <- f1 - f0
    if (!is.finite(penyebut) || abs(penyebut) < .Machine$double.eps) {
      stop("Penyebut metode secant terlalu dekat dengan nol.")
    }
    x2 <- x1 - f1 * (x1 - x0) / penyebut
    if (abs(x2 - x1) <= tol || abs(fun(x2)) <= tol) {
      return(list(root = x2, iter = i, converged = TRUE))
    }
    x0 <- x1; f0 <- f1
    x1 <- x2; f1 <- fun(x1)
  }
  list(root = x1, iter = maxit, converged = FALSE)
}
tol_akar <- 1e-10
r_bis <- bisection(f_akar, 2, 3, tol_akar)
r_new <- newton_raphson(f_akar, df_akar, 2.5, tol_akar)
r_sec <- secant(f_akar, 2, 3, tol_akar)
r_uni <- uniroot(f_akar, interval = c(2, 3), tol = tol_akar)

akar <- c(r_bis$root, r_new$root, r_sec$root, r_uni$root)
tabel_akar <- data.frame(
  Metode = c("Bisection", "Newton-Raphson", "Secant", "uniroot (verifikasi)"),
  Akar = akar,
  Iterasi = c(r_bis$iter, r_new$iter, r_sec$iter, r_uni$iter),
  Residual = abs(f_akar(akar)),
  Konvergen = c(r_bis$converged, r_new$converged, r_sec$converged, TRUE)
)
knitr::kable(tabel_akar, digits = 12, caption = "Perbandingan metode pencarian akar")
Perbandingan metode pencarian akar
Metode Akar Iterasi Residual Konvergen
Bisection 2.094551 32 0.000000000033 TRUE
Newton-Raphson 2.094551 4 0.000000000088 TRUE
Secant 2.094551 6 0.000000000000 TRUE
uniroot (verifikasi) 2.094551 7 0.000000000069 TRUE

Akar persamaan berada di sekitar 2,094551. Ketiga metode seharusnya menghasilkan hampiran yang sangat berdekatan. Jumlah iterasi dapat berbeda karena cara memperbarui hampiran berbeda; uniroot() dipakai sebagai verifikasi.

Latihan 2 Sistem Persamaan Linear

Selesaikan SPL tiga variabel dengan gauss(), jacobi(), dan solve(). Setelah itu, perkecil salah satu elemen diagonal dan amati perilaku Jacobi.

gauss <- function(A, b, tol = 1e-12) {
  A <- as.matrix(A); b <- as.numeric(b); n <- nrow(A)
  if (n != ncol(A) || length(b) != n) stop("Dimensi A dan b tidak sesuai.")
  Ab <- cbind(A, b)

  if (n > 1) {
    for (k in seq_len(n - 1)) {
      p <- which.max(abs(Ab[k:n, k])) + k - 1
      if (abs(Ab[p, k]) < tol) stop("Matriks singular atau hampir singular.")
      if (p != k) Ab[c(k, p), ] <- Ab[c(p, k), ]
      for (i in seq.int(k + 1, n)) {
        pengali <- Ab[i, k] / Ab[k, k]
        Ab[i, k:(n + 1)] <- Ab[i, k:(n + 1)] - pengali * Ab[k, k:(n + 1)]
      }
    }
  }

  if (abs(Ab[n, n]) < tol) stop("Matriks singular atau hampir singular.")
  x <- numeric(n)
  for (i in seq.int(n, 1)) {
    jumlah <- if (i < n) sum(Ab[i, (i + 1):n] * x[(i + 1):n]) else 0
    x[i] <- (Ab[i, n + 1] - jumlah) / Ab[i, i]
  }
  x
}

jacobi <- function(A, b, x0 = rep(0, length(b)),
                   tol = 1e-10, maxit = 1000, divergence_limit = 1e10) {
  A <- as.matrix(A); b <- as.numeric(b); D <- diag(A)
  if (any(abs(D) < .Machine$double.eps)) {
    return(list(x = rep(NA_real_, length(b)), iter = 0L,
                converged = FALSE, reason = "Diagonal nol atau terlalu kecil."))
  }

  R <- A - diag(D); x <- as.numeric(x0)
  for (k in seq_len(maxit)) {
    x_baru <- as.numeric((b - R %*% x) / D)
    if (any(!is.finite(x_baru)) || max(abs(x_baru)) > divergence_limit) {
      return(list(x = x_baru, iter = k, converged = FALSE,
                  reason = "Iterasi divergen atau nilainya membesar."))
    }
    if (max(abs(x_baru - x)) < tol) {
      return(list(x = x_baru, iter = k, converged = TRUE, reason = "Konvergen."))
    }
    x <- x_baru
  }
  list(x = x, iter = maxit, converged = FALSE,
       reason = "Belum konvergen sampai iterasi maksimum.")
}
A_spl <- matrix(c(10, -1, 2, -1, 11, -1, 2, -1, 10), nrow = 3, byrow = TRUE)
x_diketahui <- c(1, 2, -1)
b_spl <- as.vector(A_spl %*% x_diketahui)

x_gauss <- gauss(A_spl, b_spl)
r_jacobi <- jacobi(A_spl, b_spl)
x_solve <- solve(A_spl, b_spl)

knitr::kable(
  data.frame(Variabel = c("x1", "x2", "x3"),
             Gauss = x_gauss, Jacobi = r_jacobi$x, solve = x_solve),
  digits = 10, caption = "Perbandingan solusi SPL"
)
Perbandingan solusi SPL
Variabel Gauss Jacobi solve
x1 1 1 1
x2 2 2 2
x3 -1 -1 -1 
knitr::kable(
  data.frame(
    Metode = c("Gauss", "Jacobi", "solve"),
    Norma_Residu = c(
      sqrt(sum((A_spl %*% x_gauss - b_spl)^2)),
      sqrt(sum((A_spl %*% r_jacobi$x - b_spl)^2)),
      sqrt(sum((A_spl %*% x_solve - b_spl)^2))
    )
  ),
  digits = 12, caption = "Norma residu setiap metode"
)
Norma residu setiap metode
Metode Norma_Residu
Gauss 0.000000000000
Jacobi 0.000000000163
solve 0.000000000000 
A_kecil <- A_spl
A_kecil[1, 1] <- 0.1

knitr::kable(
  data.frame(
    Baris = 1:3,
    Nilai_Diagonal = abs(diag(A_kecil)),
    Jumlah_Nondiagonal = rowSums(abs(A_kecil)) - abs(diag(A_kecil))
  ),
  digits = 4, caption = "Dominasi diagonal setelah elemen diagonal diperkecil"
)
Dominasi diagonal setelah elemen diagonal diperkecil
Baris Nilai_Diagonal Jumlah_Nondiagonal
1 0.1 3
2 11.0 2
3 10.0 3 
r_jacobi_kecil <- jacobi(
  A_kecil, b_spl, tol = 1e-10, maxit = 100, divergence_limit = 1e8
)
knitr::kable(
  data.frame(Iterasi = r_jacobi_kecil$iter,
             Konvergen = r_jacobi_kecil$converged,
             Keterangan = r_jacobi_kecil$reason),
  caption = "Perilaku Jacobi setelah perubahan diagonal"
)
Perilaku Jacobi setelah perubahan diagonal
Iterasi Konvergen Keterangan
19 FALSE Iterasi divergen atau nilainya membesar.

Dominasi diagonal pada sistem awal mendukung konvergensi Jacobi. Ketika elemen diagonal pertama diperkecil, syarat tersebut tidak lagi terpenuhi pada baris pertama, sehingga Jacobi dapat menjadi tidak stabil atau divergen.

Latihan 3 Interpolasi Lagrange dan Spline

Interpolasikan \(f(x)=1/(1+25x^2)\) pada \([-1,1]\) menggunakan 11 titik berjarak sama. Bandingkan interpolasi Lagrange dan spline kubik.

f_runge <- function(x) 1 / (1 + 25 * x^2)

lagrange_eval <- function(x, y, xout) {
  n <- length(x)
  w <- vapply(seq_len(n), function(i) 1 / prod(x[i] - x[-i]), numeric(1))

  vapply(xout, function(t) {
    titik_sama <- which(abs(t - x) < 1e-12)
    if (length(titik_sama) > 0) return(y[titik_sama[1]])
    q <- w / (t - x)
    sum(q * y) / sum(q)
  }, numeric(1))
}

x_titik <- seq(-1, 1, length.out = 11)
y_titik <- f_runge(x_titik)
x_plot <- seq(-1, 1, length.out = 500)
y_asli <- f_runge(x_plot)
y_lagrange <- lagrange_eval(x_titik, y_titik, x_plot)
y_spline <- splinefun(x_titik, y_titik, method = "natural")(x_plot)

knitr::kable(
  data.frame(
    Metode = c("Lagrange", "Spline kubik natural"),
    Galat_Maksimum_Absolut = c(
      max(abs(y_lagrange - y_asli)),
      max(abs(y_spline - y_asli))
    )
  ),
  digits = 6, caption = "Galat maksimum terhadap fungsi asli"
)
Galat maksimum terhadap fungsi asli
Metode Galat_Maksimum_Absolut
Lagrange 1.915621
Spline kubik natural 0.021974 
plot(x_plot, y_asli, type = "l", lwd = 2, lty = 1,
     xlab = "x", ylab = "f(x)",
     main = "Interpolasi Lagrange dan Spline Kubik",
     ylim = range(c(y_asli, y_lagrange, y_spline)))
lines(x_plot, y_lagrange, lwd = 2, lty = 2)
lines(x_plot, y_spline, lwd = 2, lty = 3)
points(x_titik, y_titik, pch = 19)
legend("top", legend = c("Fungsi asli", "Lagrange derajat 10",
                         "Spline kubik", "Titik data"),
       lty = c(1, 2, 3, NA), pch = c(NA, NA, NA, 19),
       lwd = c(2, 2, 2, NA), bty = "n")

Polinom Lagrange berderajat tinggi dapat mengalami osilasi dekat batas interval pada fungsi Runge. Spline kubik menggunakan polinom berderajat rendah pada setiap subinterval dan biasanya menghasilkan kurva yang lebih stabil. Tabel galat maksimum memberi ukuran kuantitatif untuk membandingkan keduanya.

Latihan 4 Integrasi Numerik

Hitung \(\int_{0.001}^{2}\sin(x)/x\,dx\) dengan aturan trapesium, Simpson 1/3, dan integrate(). Gunakan hasil integrate() sebagai nilai acuan dan plot galat terhadap \(n\) pada skala logaritmik.

f_integral <- function(x) sin(x) / x

trapesium <- function(fun, a, b, n) {
  x <- seq(a, b, length.out = n + 1)
  h <- (b - a) / n
  h * (sum(fun(x)) - (fun(a) + fun(b)) / 2)
}

simpson <- function(fun, a, b, n) {
  if (n %% 2 != 0) stop("Simpson 1/3 membutuhkan n genap.")
  x <- seq(a, b, length.out = n + 1)
  h <- (b - a) / n
  h / 3 * (
    fun(x[1]) + fun(x[n + 1]) +
      4 * sum(fun(x[seq(2, n, by = 2)])) +
      2 * sum(fun(x[seq(3, n - 1, by = 2)]))
  )
}
batas_bawah <- 0.001
batas_atas <- 2
hasil_integrate <- integrate(
  f_integral, lower = batas_bawah, upper = batas_atas,
  rel.tol = 1e-12, abs.tol = 1e-14
)
nilai_acuan <- hasil_integrate$value
n_laporan <- 1000

hasil_trapesium <- trapesium(f_integral, batas_bawah, batas_atas, n_laporan)
hasil_simpson <- simpson(f_integral, batas_bawah, batas_atas, n_laporan)

knitr::kable(
  data.frame(
    Metode = c(paste0("Trapesium, n = ", n_laporan),
               paste0("Simpson 1/3, n = ", n_laporan), "integrate()"),
    Nilai_Integral = c(hasil_trapesium, hasil_simpson, nilai_acuan),
    Galat_Absolut_vs_integrate = c(
      abs(hasil_trapesium - nilai_acuan),
      abs(hasil_simpson - nilai_acuan),
      NA_real_
    )
  ),
  digits = 14, caption = "Perbandingan hasil integrasi"
)
Perbandingan hasil integrasi
Metode Nilai_Integral Galat_Absolut_vs_integrate
Trapesium, n = 1000 1.604413 0.00000014487650
Simpson 1/3, n = 1000 1.604413 0.00000000000002
integrate() 1.604413 NA 
n_values <- 2^(2:10)
nilai_trapesium <- sapply(
  n_values, function(n) trapesium(f_integral, batas_bawah, batas_atas, n)
)
nilai_simpson <- sapply(
  n_values, function(n) simpson(f_integral, batas_bawah, batas_atas, n)
)
galat_trapesium <- abs(nilai_trapesium - nilai_acuan)
galat_simpson <- abs(nilai_simpson - nilai_acuan)

plot(n_values, galat_trapesium, type = "b", log = "xy", lwd = 2, pch = 16,
     xlab = "Jumlah subinterval (n)",
     ylab = "Galat absolut terhadap integrate()",
     main = "Galat Integrasi pada Skala Log-Log")
lines(n_values, galat_simpson, type = "b", lwd = 2, pch = 17, lty = 2)
legend("bottomleft", legend = c("Trapesium", "Simpson 1/3"),
       lty = c(1, 2), pch = c(16, 17), lwd = 2, bty = "n")

knitr::kable(
  data.frame(n = n_values, Galat_Trapesium = galat_trapesium,
             Galat_Simpson = galat_simpson),
  digits = 12, caption = "Galat untuk berbagai nilai n"
)
Galat untuk berbagai nilai n
n Galat_Trapesium Galat_Simpson
4 0.009075373797 0.000083720013
8 0.002264978410 0.000005153385
16 0.000566003948 0.000000320873
32 0.000141485960 0.000000020036
64 0.000035370551 0.000000001252
128 0.000008842579 0.000000000078
256 0.000002210641 0.000000000005
512 0.000000552660 0.000000000000
1024 0.000000138165 0.000000000000

Aturan trapesium memiliki galat teoritis \(O(h^2)\), sedangkan Simpson 1/3 memiliki galat \(O(h^4)\) untuk fungsi yang cukup mulus. Karena itu, galat Simpson umumnya turun lebih cepat saat jumlah subinterval bertambah. integrate() menggunakan integrasi adaptif dan tidak memerlukan parameter \(n\), sehingga dipakai sebagai nilai acuan.

Latihan 5 Persamaan Diferensial Biasa

Selesaikan \(y'=-2y+t\), \(y(0)=1\), pada \(t\in[0,2]\) dengan Euler dan RK4 untuk \(h=0.2\) dan \(h=0.05\). Bandingkan dengan ode() dari paket deSolve.

f_pdb <- function(t, y) -2 * y + t
solusi_eksak <- function(t) 0.5 * t - 0.25 + 1.25 * exp(-2 * t)

euler_pdb <- function(fun, t0, y0, h, t_akhir) {
  n <- round((t_akhir - t0) / h)
  waktu <- seq(t0, by = h, length.out = n + 1)
  y <- numeric(n + 1)
  y[1] <- y0
  for (i in seq_len(n)) y[i + 1] <- y[i] + h * fun(waktu[i], y[i])
  data.frame(waktu = waktu, y = y)
}

rk4_pdb <- function(fun, t0, y0, h, t_akhir) {
  n <- round((t_akhir - t0) / h)
  waktu <- seq(t0, by = h, length.out = n + 1)
  y <- numeric(n + 1)
  y[1] <- y0
  for (i in seq_len(n)) {
    k1 <- fun(waktu[i], y[i])
    k2 <- fun(waktu[i] + h / 2, y[i] + h * k1 / 2)
    k3 <- fun(waktu[i] + h / 2, y[i] + h * k2 / 2)
    k4 <- fun(waktu[i] + h, y[i] + h * k3)
    y[i + 1] <- y[i] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
  }
  data.frame(waktu = waktu, y = y)
}

ode_acuan <- function(waktu) {
  hasil <- deSolve::ode(
    y = c(y = 1), times = waktu,
    func = function(t, y, parms) list(-2 * y[1] + t),
    parms = NULL
  )
  as.data.frame(hasil)
}
bandingkan_pdb <- function(h) {
  waktu <- seq(0, 2, by = h)
  E <- euler_pdb(f_pdb, 0, 1, h, 2)
  R <- rk4_pdb(f_pdb, 0, 1, h, 2)
  O <- ode_acuan(waktu)
  data.frame(h = h, waktu = waktu, Eksak = solusi_eksak(waktu),
             Euler = E$y, RK4 = R$y, deSolve = O$y)
}

hasil_pdb <- do.call(rbind, lapply(c(0.2, 0.05), bandingkan_pdb))
hasil_pdb$Galat_Euler <- abs(hasil_pdb$Euler - hasil_pdb$deSolve)
hasil_pdb$Galat_RK4 <- abs(hasil_pdb$RK4 - hasil_pdb$deSolve)

knitr::kable(
  aggregate(cbind(Galat_Euler, Galat_RK4) ~ h, data = hasil_pdb, FUN = max),
  digits = 10, caption = "Galat maksimum terhadap hasil deSolve()"
)
Galat maksimum terhadap hasil deSolve()
h Galat_Euler Galat_RK4
0.05 0.02400155 0.0000015320
0.20 0.11166138 0.0001338167 
knitr::kable(head(hasil_pdb, 12), digits = 8,
             caption = "Sebagian hasil numerik untuk h = 0.2")
Sebagian hasil numerik untuk h = 0.2
h waktu Eksak Euler RK4 deSolve Galat_Euler Galat_RK4
0.20 0.0 1.0000000 1.0000000 1.0000000 1.0000000 0.00000000 0.00000000
0.20 0.2 0.6879001 0.6000000 0.6880000 0.6879016 0.08790161 0.00009839
0.20 0.4 0.5116612 0.4000000 0.5117952 0.5116614 0.11166138 0.00013382
0.20 0.6 0.4264928 0.3200000 0.4266275 0.4264950 0.10649496 0.00013254
0.20 0.8 0.4023706 0.3120000 0.4024911 0.4023718 0.09037177 0.00011930
0.20 1.0 0.4191691 0.3472000 0.4192700 0.4191699 0.07196992 0.00010010
0.20 1.2 0.4633974 0.4083200 0.4634786 0.4633983 0.05507826 0.00008036
0.20 1.4 0.5260126 0.4849920 0.5260761 0.5260133 0.04102130 0.00006277
0.20 1.6 0.6009528 0.5709952 0.6010014 0.6009534 0.02995816 0.00004804
0.20 1.8 0.6841546 0.6625971 0.6841913 0.6841554 0.02155831 0.00003590
0.20 2.0 0.7728945 0.7575583 0.7729219 0.7728951 0.01533685 0.00002674
0.05 0.0 1.0000000 1.0000000 1.0000000 1.0000000 0.00000000 0.00000000 
par(mfrow = c(1, 2))
for (h_sekarang in c(0.2, 0.05)) {
  d <- subset(hasil_pdb, h == h_sekarang)
  plot(d$waktu, d$deSolve, type = "l", lwd = 2, lty = 1,
       xlab = "Waktu (t)", ylab = "y(t)",
       main = paste("Perbandingan metode, h =", h_sekarang))
  lines(d$waktu, d$Eksak, lwd = 2, lty = 2)
  lines(d$waktu, d$Euler, lwd = 2, lty = 3)
  lines(d$waktu, d$RK4, lwd = 2, lty = 4)
  legend("topright", legend = c("deSolve", "Eksak", "Euler", "RK4"),
         lty = 1:4, lwd = 2, bty = "n", cex = 0.8)
}

par(mfrow = c(1, 1))

Euler merupakan metode orde satu, sedangkan RK4 merupakan metode orde empat. Umumnya RK4 memberikan hampiran lebih dekat terhadap solusi acuan untuk ukuran langkah yang sama. Pengurangan \(h\) juga diharapkan memperkecil galat Euler.

Latihan 6 Optimasi Tanpa Kendala

Minimumkan \(f(x,y)=x^4+y^4-4xy+1\) dengan gradient descent, Newton multivariabel, dan optim() metode BFGS. Tentukan seluruh titik stasioner dan periksa Hessian-nya.

f_opt <- function(p) {
  x <- p[1]; y <- p[2]
  x^4 + y^4 - 4 * x * y + 1
}

grad_opt <- function(p) {
  x <- p[1]; y <- p[2]
  c(4 * x^3 - 4 * y, 4 * y^3 - 4 * x)
}

hess_opt <- function(p) {
  x <- p[1]; y <- p[2]
  matrix(c(12 * x^2, -4, -4, 12 * y^2), nrow = 2, byrow = TRUE)
}

gradient_descent <- function(start, alpha = 0.01, tol = 1e-8, maxit = 100000) {
  p <- as.numeric(start)
  konvergen <- FALSE
  for (i in seq_len(maxit)) {
    g <- grad_opt(p)
    if (sqrt(sum(g^2)) < tol) {
      konvergen <- TRUE
      break
    }
    p_baru <- p - alpha * g
    if (any(!is.finite(p_baru))) break
    p <- p_baru
  }
  list(par = p, value = f_opt(p), iter = i, converged = konvergen,
       grad_norm = sqrt(sum(grad_opt(p)^2)))
}

newton_optimasi <- function(start, tol = 1e-10, maxit = 100) {
  p <- as.numeric(start)
  konvergen <- FALSE
  for (i in seq_len(maxit)) {
    g <- grad_opt(p)
    if (sqrt(sum(g^2)) < tol) {
      konvergen <- TRUE
      break
    }

    langkah <- tryCatch(solve(hess_opt(p), g), error = function(e) NULL)
    if (is.null(langkah)) break

    # Backtracking sederhana agar langkah tidak menaikkan fungsi objektif.
    skala <- 1
    while (skala > 1e-8 && f_opt(p - skala * langkah) > f_opt(p)) {
      skala <- skala / 2
    }
    if (skala <= 1e-8) break
    p <- p - skala * langkah
  }
  list(par = p, value = f_opt(p), iter = i, converged = konvergen,
       grad_norm = sqrt(sum(grad_opt(p)^2)))
}
bandingkan_optimasi <- function(start, label_start) {
  gd <- gradient_descent(start)
  nw <- newton_optimasi(start)
  bg <- optim(
    par = start, fn = f_opt, gr = grad_opt, method = "BFGS",
    control = list(reltol = 1e-12, maxit = 10000)
  )

  rbind(
    data.frame(Titik_Awal = label_start, Metode = "Gradient descent",
               x = gd$par[1], y = gd$par[2], Nilai_f = gd$value,
               Iterasi = gd$iter, Evaluasi_Fungsi = NA_integer_,
               Konvergen = gd$converged),
    data.frame(Titik_Awal = label_start, Metode = "Newton multivariabel",
               x = nw$par[1], y = nw$par[2], Nilai_f = nw$value,
               Iterasi = nw$iter, Evaluasi_Fungsi = NA_integer_,
               Konvergen = nw$converged),
    data.frame(Titik_Awal = label_start, Metode = "optim BFGS",
               x = bg$par[1], y = bg$par[2], Nilai_f = bg$value,
               Iterasi = NA_integer_, Evaluasi_Fungsi = unname(bg$counts[1]),
               Konvergen = (bg$convergence == 0))
  )
}

hasil_optimasi <- rbind(
  bandingkan_optimasi(c(2, 2), "(2, 2)"),
  bandingkan_optimasi(c(-2, -2), "(-2, -2)")
)
knitr::kable(hasil_optimasi, digits = 10,
             caption = "Hasil tiga metode optimasi dari dua titik awal")
Hasil tiga metode optimasi dari dua titik awal
Titik_Awal Metode x y Nilai_f Iterasi Evaluasi_Fungsi Konvergen
(2, 2) Gradient descent 1.0000000 1.0000000 -1 239 NA TRUE
(2, 2) Newton multivariabel 1.0000000 1.0000000 -1 7 NA TRUE
(2, 2) optim BFGS 0.9999999 0.9999999 -1 NA 10 TRUE
(-2, -2) Gradient descent -1.0000000 -1.0000000 -1 239 NA TRUE
(-2, -2) Newton multivariabel -1.0000000 -1.0000000 -1 7 NA TRUE
(-2, -2) optim BFGS -0.9999999 -0.9999999 -1 NA 10 TRUE

Gradien sama dengan nol jika \(4x^3-4y=0\) dan \(4y^3-4x=0\). Titik stasioner riilnya adalah \((0,0)\), \((1,1)\), dan \((-1,-1)\). Klasifikasi dilakukan berdasarkan nilai eigen Hessian.

titik_stasioner <- data.frame(x = c(0, 1, -1), y = c(0, 1, -1))
titik_stasioner$Nilai_f <- apply(
  titik_stasioner[, c("x", "y")], 1,
  function(p) f_opt(as.numeric(p))
)

nilai_eigen <- t(vapply(
  seq_len(nrow(titik_stasioner)),
  function(i) eigen(
    hess_opt(c(titik_stasioner$x[i], titik_stasioner$y[i])),
    symmetric = TRUE, only.values = TRUE
  )$values,
  numeric(2)
))

titik_stasioner$Eigen_Min <- apply(nilai_eigen, 1, min)
titik_stasioner$Eigen_Maks <- apply(nilai_eigen, 1, max)
titik_stasioner$Klasifikasi <- ifelse(
  titik_stasioner$Eigen_Min > 0, "Minimum lokal",
  ifelse(titik_stasioner$Eigen_Maks < 0, "Maksimum lokal", "Titik pelana")
)

knitr::kable(titik_stasioner, digits = 6,
             caption = "Klasifikasi titik stasioner berdasarkan Hessian")
Klasifikasi titik stasioner berdasarkan Hessian
x y Nilai_f Eigen_Min Eigen_Maks Klasifikasi
0 0 1 -4 4 Titik pelana
1 1 -1 8 16 Minimum lokal
-1 -1 -1 8 16 Minimum lokal

Titik \((1,1)\) dan \((-1,-1)\) memiliki nilai fungsi \(-1\) dan Hessian definit positif, sehingga keduanya merupakan minimum lokal sekaligus minimum global. Titik \((0,0)\) adalah titik pelana karena Hessian memiliki nilai eigen dengan tanda berbeda.

Latihan 7 Optimasi Berkendala dengan Pemrograman Linear

Keuntungan per unit produk \(x\) dan \(y\) adalah 40 dan 30. Kendala: \(2x+y\leq100\), \(x+y\leq80\), \(x\leq40\), dan \(x,y\geq0\). Tentukan produksi optimal menggunakan lpSolve dan gambarkan daerah layak.

keuntungan <- c(40, 30)
kendala_lp <- rbind(c(2, 1), c(1, 1), c(1, 0))
arah_kendala <- c("<=", "<=", "<=")
ruas_kanan <- c(100, 80, 40)

hasil_lp <- lpSolve::lp(
  direction = "max",
  objective.in = keuntungan,
  const.mat = kendala_lp,
  const.dir = arah_kendala,
  const.rhs = ruas_kanan
)

solusi_x <- hasil_lp$solution[1]
solusi_y <- hasil_lp$solution[2]

knitr::kable(
  data.frame(
    Status = ifelse(hasil_lp$status == 0, "Optimal", "Periksa status solver"),
    Produk_X = solusi_x,
    Produk_Y = solusi_y,
    Keuntungan_Maksimum = hasil_lp$objval
  ),
  digits = 6, caption = "Solusi optimal pemrograman linear"
)
Solusi optimal pemrograman linear
Status Produk_X Produk_Y Keuntungan_Maksimum
Optimal 20 60 2600 
pemakaian <- as.vector(kendala_lp %*% hasil_lp$solution)
knitr::kable(
  data.frame(
    Kendala = c("2x + y", "x + y", "x"),
    Nilai_Kiri = pemakaian,
    Batas = ruas_kanan,
    Sisa_Kapasitas = ruas_kanan - pemakaian
  ),
  digits = 6, caption = "Pemeriksaan kendala pada solusi optimal"
)
Pemeriksaan kendala pada solusi optimal
Kendala Nilai_Kiri Batas Sisa_Kapasitas
2x + y 100 100 0
x + y 80 80 0
x 20 40 20 
plot(NA, xlim = c(0, 45), ylim = c(0, 90),
     xlab = "Jumlah produk X", ylab = "Jumlah produk Y",
     main = "Daerah Layak Pemrograman Linear")

# Titik sudut daerah layak: (0,0), (40,0), (40,20), (20,60), (0,80).
polygon(
  x = c(0, 40, 40, 20, 0),
  y = c(0, 0, 20, 60, 80),
  col = gray(0.9), border = "gray40"
)

xx <- seq(0, 45, length.out = 300)
lines(xx, 100 - 2 * xx, lty = 1, lwd = 2)
lines(xx, 80 - xx, lty = 2, lwd = 2)
abline(v = 40, lty = 3, lwd = 2)
points(solusi_x, solusi_y, pch = 19, cex = 1.3)
text(solusi_x, solusi_y, labels = "(20, 60)", pos = 4)

legend(
  "topright",
  legend = c("2x + y = 100", "x + y = 80", "x = 40", "Daerah layak", "Optimum"),
  lty = c(1, 2, 3, 0, 0), pch = c(NA, NA, NA, 15, 19),
  lwd = c(2, 2, 2, 0, 0), pt.cex = 1.2, bty = "n"
)

Solusi optimal adalah 20 unit produk X dan 60 unit produk Y, dengan keuntungan maksimum 2600 satuan keuntungan. Solusi tersebut memenuhi semua kendala dan berada pada perpotongan \(2x+y=100\) dan \(x+y=80\).

Latihan 8 Perbandingan Metode Optimasi Global

Bandingkan BFGS, Nelder-Mead, SANN, dan algoritma genetika (GA) pada fungsi Rastrigin lima dimensi. Lakukan 20 ulangan acak dan laporkan rata-rata serta simpangan baku nilai akhir dan waktu komputasi dari system.time().

rastrigin <- function(x) {
  10 * length(x) + sum(x^2 - 10 * cos(2 * pi * x))
}

grad_rastrigin <- function(x) {
  2 * x + 20 * pi * sin(2 * pi * x)
}

# Minimum global fungsi Rastrigin berada di (0, 0, ..., 0), dengan nilai 0.

Menjalankan 20 ulangan setiap metode

Pengaturan dibuat tetap agar hasil dapat direproduksi. BFGS dan Nelder-Mead memakai titik awal acak yang sama pada setiap ulangan. SANN dan GA menggunakan bilangan acak internal; seed diatur pada setiap ulangan.

jumlah_ulangan <- 20
metode_uji <- c("BFGS", "Nelder-Mead", "SANN", "GA")

set.seed(2026)
titik_awal <- replicate(
  jumlah_ulangan,
  runif(5, min = -5.12, max = 5.12),
  simplify = FALSE
)

hasil_benchmark <- vector("list", length(metode_uji) * jumlah_ulangan)
indeks_hasil <- 1

for (nama_metode in metode_uji) {
  for (i in seq_len(jumlah_ulangan)) {
    set.seed(100000 + 1000 * match(nama_metode, metode_uji) + i)
    awal <- titik_awal[[i]]

    if (nama_metode == "BFGS") {
      waktu <- system.time({
        fit <- optim(
          par = awal, fn = rastrigin, gr = grad_rastrigin,
          method = "BFGS",
          control = list(maxit = 2000, reltol = 1e-10)
        )
        nilai_akhir <- fit$value
      })[["elapsed"]]

    } else if (nama_metode == "Nelder-Mead") {
      waktu <- system.time({
        fit <- optim(
          par = awal, fn = rastrigin, method = "Nelder-Mead",
          control = list(maxit = 5000, reltol = 1e-10)
        )
        nilai_akhir <- fit$value
      })[["elapsed"]]

    } else if (nama_metode == "SANN") {
      waktu <- system.time({
        fit <- optim(
          par = awal, fn = rastrigin, method = "SANN",
          control = list(maxit = 20000, temp = 20, tmax = 10, trace = 0)
        )
        nilai_akhir <- fit$value
      })[["elapsed"]]

    } else if (nama_metode == "GA") {
      waktu <- system.time({
        fit_ga <- GA::ga(
          type = "real-valued",
          fitness = function(x) -rastrigin(x),
          lower = rep(-5.12, 5),
          upper = rep(5.12, 5),
          popSize = 50,
          maxiter = 150,
          run = 30,
          monitor = FALSE
        )
        solusi_ga <- as.numeric(fit_ga@solution[1, ])
        nilai_akhir <- rastrigin(solusi_ga)
      })[["elapsed"]]
    }

    hasil_benchmark[[indeks_hasil]] <- data.frame(
      Metode = nama_metode,
      Ulangan = i,
      Nilai_Akhir = as.numeric(nilai_akhir),
      Waktu_Detik = as.numeric(waktu)
    )
    indeks_hasil <- indeks_hasil + 1
  }
}

data_benchmark <- do.call(rbind, hasil_benchmark)
data_benchmark$Metode <- factor(data_benchmark$Metode, levels = metode_uji)

knitr::kable(
  head(data_benchmark, 12), digits = 6,
  caption = "Sebagian hasil dari 20 ulangan setiap metode"
)
Sebagian hasil dari 20 ulangan setiap metode
Metode Ulangan Nilai_Akhir Waktu_Detik
BFGS 1 22.88395 0.03
BFGS 2 50.74253 0.00
BFGS 3 69.64664 0.00
BFGS 4 54.72224 0.00
BFGS 5 62.68200 0.00
BFGS 6 51.73749 0.00
BFGS 7 47.75772 0.00
BFGS 8 19.89912 0.00
BFGS 9 50.74234 0.00
BFGS 10 41.78799 0.00
BFGS 11 52.73252 0.00
BFGS 12 63.67684 0.00

Rata-rata, simpangan baku, dan waktu komputasi

rata_rata <- aggregate(
  cbind(Nilai_Akhir, Waktu_Detik) ~ Metode,
  data = data_benchmark, FUN = mean
)
simpangan_baku <- aggregate(
  cbind(Nilai_Akhir, Waktu_Detik) ~ Metode,
  data = data_benchmark, FUN = sd
)

tabel_benchmark <- data.frame(
  Metode = rata_rata$Metode,
  Rata_Rata_Nilai_Akhir = rata_rata$Nilai_Akhir,
  SD_Nilai_Akhir = simpangan_baku$Nilai_Akhir,
  Rata_Rata_Waktu_Detik = rata_rata$Waktu_Detik,
  SD_Waktu_Detik = simpangan_baku$Waktu_Detik
)

knitr::kable(
  tabel_benchmark, digits = 6,
  caption = "Ringkasan kinerja metode optimasi Rastrigin 5 dimensi"
)
Ringkasan kinerja metode optimasi Rastrigin 5 dimensi
Metode Rata_Rata_Nilai_Akhir SD_Nilai_Akhir Rata_Rata_Waktu_Detik SD_Waktu_Detik
BFGS 45.519026 19.025927 0.0015 0.006708
Nelder-Mead 45.270328 20.463082 0.0020 0.006156
SANN 24.687925 13.702407 0.1260 0.040833
GA 0.559775 0.926565 0.3565 0.125416

Grafik distribusi nilai akhir dan waktu komputasi

par(mfrow = c(1, 2))
boxplot(
  Nilai_Akhir ~ Metode, data = data_benchmark,
  xlab = "Metode", ylab = "Nilai fungsi akhir",
  main = "Distribusi Nilai Akhir", las = 2
)
boxplot(
  Waktu_Detik ~ Metode, data = data_benchmark,
  xlab = "Metode", ylab = "Waktu (detik)",
  main = "Distribusi Waktu Komputasi", las = 2
)

par(mfrow = c(1, 1))

Fungsi Rastrigin mempunyai banyak minimum lokal dan minimum global bernilai 0 pada titik asal. BFGS dan Nelder-Mead dapat berhenti pada minimum lokal bergantung pada titik awal. SANN dan GA dirancang untuk mengeksplorasi ruang pencarian secara lebih global, tetapi hasilnya tetap dipengaruhi parameter dan bilangan acak. Kesimpulan akhir perlu didasarkan pada rata-rata, simpangan baku, dan waktu pada tabel, bukan hanya satu ulangan.

Kesimpulan

  1. Metode pencarian akar dapat dibandingkan melalui residual dan jumlah iterasi.
  2. Konvergensi Jacobi dipengaruhi sifat matriks; dominasi diagonal mendukung konvergensi.
  3. Interpolasi polinom berderajat tinggi dapat berosilasi, sedangkan spline kubik biasanya lebih stabil.
  4. Galat integrasi berkurang saat ukuran langkah diperkecil; Simpson umumnya lebih cepat akurat daripada trapesium untuk fungsi yang mulus.
  5. RK4 umumnya lebih akurat daripada Euler untuk ukuran langkah yang sama.
  6. Hessian membantu membedakan minimum lokal dan titik pelana.
  7. Solusi pemrograman linear perlu diperiksa kembali terhadap seluruh kendala.
  8. Perbandingan metode optimasi global sebaiknya memakai beberapa ulangan dan melaporkan kualitas solusi serta waktu komputasi.

Referensi

Virgania Sari, S.Pd., M.Si. Komputasi Numerik dan Optimasi: Modul Materi dan Praktikum Menggunakan R. Mata Kuliah Komputasi Statistika.

Dokumentasi paket R: deSolve, lpSolve, GA, dan fungsi dasar R seperti optim(), uniroot(), integrate(), serta solve().