1.2 Praktikum Galat dan Aritmetika Komputer

0.1 + 0.2 == 0.3                      # FALSE
## [1] FALSE
isTRUE(all.equal(0.1 + 0.2, 0.3))     # TRUE, dibandingkan dengan toleransi
## [1] TRUE
.Machine$double.eps                   # epsilon mesin sekitar 2.22e-16
## [1] 2.220446e-16

Galat hampiran pi oleh 22/7

x <- pi; xh <- 22/7
c(mutlak = abs(x - xh), relatif = abs(x - xh) / abs(x))
##       mutlak      relatif 
## 0.0012644893 0.0004024994

Galat pemotongan deret Taylor e^x di x = 1

taylor_exp <- function(x, n) sum(x^(0:n) / factorial(0:n))
sapply(1:8, function(n) abs(exp(1) - taylor_exp(1, n)))
## [1] 7.182818e-01 2.182818e-01 5.161516e-02 9.948495e-03 1.615162e-03
## [6] 2.262729e-04 2.786021e-05 3.058618e-06

Cancellation: pengurangan dua bilangan hampir sama

f1 <- function(x) sqrt(x + 1) - sqrt(x)
f2 <- function(x) 1 / (sqrt(x + 1) + sqrt(x))   # bentuk ekuivalen yang stabil
x <- 1e15
print(c(f1(x), f2(x)), digits = 15)
## [1] 1.86264514923096e-08 1.58113883008419e-08

2. Akar Persamaan Nonlinear

Contoh: \(f(x) = x^3 - x - 2\) dengan akar di sekitar 1,5214.

2.1 Metode Bisection

bisection <- function(f, a, b, tol = 1e-8, maxit = 100) {
  if (f(a) * f(b) > 0) stop('f(a) dan f(b) harus berbeda tanda')
  for (i in 1:maxit) {
    c <- (a + b) / 2
    if (abs(f(c)) < tol || (b - a) / 2 < tol) return(list(akar = c, iter = i))
    if (f(a) * f(c) < 0) b <- c else a <- c
  }
  list(akar = c, iter = maxit)
}

f <- function(x) x^3 - x - 2
bisection(f, 1, 2)
## $akar
## [1] 1.52138
## 
## $iter
## [1] 27

2.2 Metode Regula Falsi

regula_falsi <- function(f, a, b, tol = 1e-8, maxit = 100) {
  for (i in 1:maxit) {
    c <- b - f(b) * (b - a) / (f(b) - f(a))
    if (abs(f(c)) < tol) return(list(akar = c, iter = i))
    if (f(a) * f(c) < 0) b <- c else a <- c
  }
  list(akar = c, iter = maxit)
}

regula_falsi(f, 1, 2)
## $akar
## [1] 1.52138
## 
## $iter
## [1] 17

2.3 Metode Newton-Raphson

newton <- function(f, df, x0, tol = 1e-10, maxit = 100) {
  x <- x0
  for (i in 1:maxit) {
    xn <- x - f(x) / df(x)
    if (abs(xn - x) < tol) return(list(akar = xn, iter = i))
    x <- xn
  }
  warning('belum konvergen')
  list(akar = x, iter = maxit)
}

df <- function(x) 3 * x^2 - 1
newton(f, df, x0 = 1.5)
## $akar
## [1] 1.52138
## 
## $iter
## [1] 4
# Turunan simbolik otomatis dengan D()
ex <- expression(x^3 - x - 2)
df2 <- function(x) eval(D(ex[[1]], 'x'))
df2(1.5)
## [1] 5.75

2.4 Metode Secant

secant <- function(f, x0, x1, tol = 1e-10, maxit = 100) {
  for (i in 1:maxit) {
    x2 <- x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
    if (abs(x2 - x1) < tol) return(list(akar = x2, iter = i))
    x0 <- x1; x1 <- x2
  }
  list(akar = x1, iter = maxit)
}

secant(f, 1, 2)
## $akar
## [1] 1.52138
## 
## $iter
## [1] 8

2.5 Iterasi Titik Tetap

fixed_point <- function(g, x0, tol = 1e-10, maxit = 1000) {
  x <- x0
  for (i in 1:maxit) {
    xn <- g(x)
    if (abs(xn - x) < tol) return(list(akar = xn, iter = i))
    x <- xn
  }
  list(akar = x, iter = maxit)
}

fixed_point(function(x) (x + 2)^(1/3), 1.5)   # dari x^3 = x + 2
## $akar
## [1] 1.52138
## 
## $iter
## [1] 11

2.6 Fungsi bawaan R dan perbandingan

uniroot(f, interval = c(1, 2), tol = 1e-10)$root
## [1] 1.52138
polyroot(c(-2, -1, 0, 1))   # koefisien urut naik
## [1]  1.5213797+1.392310e-20i -0.7606899+8.578736e-01i -0.7606899-8.578736e-01i
# Perbandingan banyak iterasi
data.frame(
  metode  = c('Bisection', 'Regula Falsi', 'Newton', 'Secant'),
  iterasi = c(bisection(f, 1, 2, 1e-10)$iter,
              regula_falsi(f, 1, 2, 1e-10)$iter,
              newton(f, df, 1.5)$iter,
              secant(f, 1, 2)$iter)
)
##         metode iterasi
## 1    Bisection      33
## 2 Regula Falsi      20
## 3       Newton       4
## 4       Secant       8

3. Sistem Persamaan Linear (SPL)

3.1 Eliminasi Gauss dengan pivoting parsial

gauss <- function(A, b) {
  n <- nrow(A); Ab <- cbind(A, b)
  for (k in 1:(n - 1)) {
    p <- which.max(abs(Ab[k:n, k])) + k - 1
    if (p != k) Ab[c(k, p), ] <- Ab[c(p, k), ]
    for (i in (k + 1):n) {
      m <- Ab[i, k] / Ab[k, k]
      Ab[i, ] <- Ab[i, ] - m * Ab[k, ]
    }
  }
  x <- numeric(n)
  for (i in n:1) {
    s <- if (i < n) sum(Ab[i, (i + 1):n] * x[(i + 1):n]) else 0
    x[i] <- (Ab[i, n + 1] - s) / Ab[i, i]
  }
  x
}

A <- matrix(c(4, -1, 0, -1, 4, -1, 0, -1, 4), 3, byrow = TRUE)
b <- c(15, 10, 10)
gauss(A, b)
## [1] 4.910714 4.642857 3.660714
solve(A, b)   # pembanding bawaan R
## [1] 4.910714 4.642857 3.660714

3.2 Dekomposisi LU, Cholesky, QR

lu_res <- expand(lu(Matrix(A)))   # paket Matrix: A = P L U
lu_res$L; lu_res$U
## 3 x 3 Matrix of class "dtrMatrix" (unitriangular)
##      [,1]       [,2]       [,3]      
## [1,]  1.0000000          .          .
## [2,] -0.2500000  1.0000000          .
## [3,]  0.0000000 -0.2666667  1.0000000
## 3 x 3 Matrix of class "dtrMatrix"
##      [,1]      [,2]      [,3]     
## [1,]  4.000000 -1.000000  0.000000
## [2,]         .  3.750000 -1.000000
## [3,]         .         .  3.733333
chol(A)         # A simetrik definit positif: A = R'R
##      [,1]      [,2]       [,3]
## [1,]    2 -0.500000  0.0000000
## [2,]    0  1.936492 -0.5163978
## [3,]    0  0.000000  1.9321836
qr.solve(A, b)  # via dekomposisi QR
## [1] 4.910714 4.642857 3.660714

3.3 Metode iteratif: Jacobi dan Gauss-Seidel

jacobi <- function(A, b, x0 = rep(0, length(b)), tol = 1e-8, maxit = 1000) {
  D <- diag(A); R <- A - diag(D); x <- x0
  for (k in 1:maxit) {
    xn <- (b - R %*% x) / D
    if (max(abs(xn - x)) < tol) return(list(x = as.vector(xn), iter = k))
    x <- xn
  }
  list(x = as.vector(x), iter = maxit)
}

gauss_seidel <- function(A, b, x0 = rep(0, length(b)), tol = 1e-8, maxit = 1000) {
  n <- length(b); x <- x0
  for (k in 1:maxit) {
    xold <- x
    for (i in 1:n) x[i] <- (b[i] - sum(A[i, -i] * x[-i])) / A[i, i]
    if (max(abs(x - xold)) < tol) return(list(x = x, iter = k))
  }
  list(x = x, iter = maxit)
}

jacobi(A, b)
## $x
## [1] 4.910714 4.642857 3.660714
## 
## $iter
## [1] 21
gauss_seidel(A, b)
## $x
## [1] 4.910714 4.642857 3.660714
## 
## $iter
## [1] 12

3.4 Kondisi matriks

H <- 1 / outer(1:6, 1:6, '+')   # matriks mirip Hilbert
kappa(H); kappa(A)
## [1] 61143998
## [1] 2.308061

4. Interpolasi dan Pencocokan Kurva

4.1 Interpolasi Lagrange

lagrange <- function(x, y, xp) {
  sapply(xp, function(t) sum(sapply(seq_along(x), function(i)
    y[i] * prod((t - x[-i]) / (x[i] - x[-i])))))
}

x <- c(0, 1, 2, 3); y <- c(1, 2, 5, 10)   # y = x^2 + 1
lagrange(x, y, c(0.5, 1.5, 2.5))
## [1] 1.25 3.25 7.25

4.2 Interpolasi Newton (selisih terbagi)

newton_dd <- function(x, y) {
  n <- length(x); d <- y; coef <- numeric(n); coef[1] <- d[1]
  for (j in 2:n) {
    d <- (d[-1] - d[-length(d)]) / (x[j:n] - x[1:(n - j + 1)])
    coef[j] <- d[1]
  }
  coef
}

newton_eval <- function(x, coef, xp) {
  n <- length(coef)
  sapply(xp, function(t) {
    r <- coef[n]
    if (n > 1) for (i in (n - 1):1) r <- r * (t - x[i]) + coef[i]
    r
  })
}

cf <- newton_dd(x, y)
newton_eval(x, cf, c(0.5, 1.5, 2.5))
## [1] 1.25 3.25 7.25

4.3 Spline dan interpolasi linear bawaan

xs <- seq(-5, 5, length.out = 11); ys <- 1 / (1 + xs^2)
xx <- seq(-5, 5, length.out = 200)

plot(xx, 1 / (1 + xx^2), type = 'l', lwd = 2, main = 'Runge vs Spline')
lines(xx, lagrange(xs, ys, xx), col = 'red')    # berosilasi di tepi
lines(spline(xs, ys, n = 200), col = 'blue', lwd = 2)   # spline kubik
points(xs, ys, pch = 19)

approx(xs, ys, xout = 0.5)   # interpolasi linear
## $x
## [1] 0.5
## 
## $y
## [1] 0.75

4.4 Regresi kuadrat terkecil

set.seed(1)
x <- 1:20; y <- 3 + 2 * x + 0.1 * x^2 + rnorm(20, sd = 2)
X <- cbind(1, x, x^2)
beta <- solve(t(X) %*% X, t(X) %*% y)   # persamaan normal
beta
##        [,1]
##   3.0382456
## x 2.0130475
##   0.1014342
coef(lm(y ~ x + I(x^2)))   # pembanding
## (Intercept)           x      I(x^2) 
##   3.0382456   2.0130475   0.1014342
# Regresi nonlinear: y = a * exp(b * x)
xd <- 1:10; yd <- 2 * exp(0.3 * xd) * (1 + rnorm(10, sd = 0.02))
fit <- nls(yd ~ a * exp(b * xd), start = list(a = 1, b = 0.1))
summary(fit)
## 
## Formula: yd ~ a * exp(b * xd)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 1.933451   0.050086    38.6 2.23e-10 ***
## b 0.303501   0.002891   105.0 7.57e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2667 on 8 degrees of freedom
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 9.979e-06