# Memanggil package yang digunakan
library(deSolve); library(lpSolve); library(GA); library(Matrix); library(numDeriv)
## Warning: package 'deSolve' was built under R version 4.5.3
## Warning: package 'GA' was built under R version 4.5.3
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 4.5.3
## Loading required package: iterators
## Warning: package 'iterators' was built under R version 4.5.3
## 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
## Warning: package 'Matrix' was built under R version 4.5.3
Sub Bab 1 Galat dan Aritmetika Komputer
# Membandingkan bilangan desimal
0.1 + 0.2 == 0.3
## [1] FALSE
# Membandingkan dengan toleransi
isTRUE(all.equal(0.1 + 0.2, 0.3))
## [1] TRUE
# Epsilon mesin
.Machine$double.eps
## [1] 2.220446e-16
# Nilai eksak dan hampiran
x <- pi
xh <- 22/7
# Menghitung galat
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
# Fungsi hampiran deret Taylor untuk e^x
taylor_exp <- function(x, n) {
sum(x^(0:n) / factorial(0:n))
}
# Fungsi hampiran deret Taylor untuk e^x
taylor_exp <- function(x, n) {
sum(x^(0:n) / factorial(0:n))
}
# Menghitung galat untuk beberapa jumlah suku
hasil_taylor <- sapply(1:8, function(n) {
abs(exp(1) - taylor_exp(1, n))
})
# Menampilkan hasil
data.frame(
jumlah_suku = 1:8,
galat = hasil_taylor
)
## jumlah_suku galat
## 1 1 7.182818e-01
## 2 2 2.182818e-01
## 3 3 5.161516e-02
## 4 4 9.948495e-03
## 5 5 1.615162e-03
## 6 6 2.262729e-04
## 7 7 2.786021e-05
## 8 8 3.058618e-06
f1 <- function(x) sqrt(x + 1) - sqrt(x)
f2 <- function(x) {
1 / (sqrt(x + 1) + sqrt(x))
}
x <- 1e15
print(c(
bentuk_langsung = f1(x),
bentuk_stabil = f2(x)
), digits = 15)
## bentuk_langsung bentuk_stabil
## 1.86264514923096e-08 1.58113883008419e-08
Sub Bab 2 Akar Persamaan Nonlinear
# Fungsi yang akan dicari akarnya
f <- function(x) {
x^3 - x - 2
}
# 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)
}
# Menjalankan metode Bisection
bisection(f, 1, 2)
## $akar
## [1] 1.52138
##
## $iter
## [1] 27
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
# Turunan pertama f(x)
df <- function(x) {
3 * x^2 - 1
}
# 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)
}
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'))
# Menampilkan hasil turunan saat x = 2
df2(2)
## [1] 11
# Menampilkan hasil turunan saat x = 3
df2(3)
## [1] 26
# 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
# 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
# Hasil masing-masing metode
hasil_bis <- bisection(f, 1, 2, 1e-10)
hasil_rf <- regula_falsi(f, 1, 2, 1e-10)
hasil_new <- newton(f, df, 1.5)
hasil_sec <- secant(f, 1, 2)
# Tabel perbandingan
perbandingan <- data.frame(
metode = c(
"Bisection",
"Regula Falsi",
"Newton-Raphson",
"Secant"
),
akar = c(
hasil_bis$akar,
hasil_rf$akar,
hasil_new$akar,
hasil_sec$akar
),
iterasi = c(
hasil_bis$iter,
hasil_rf$iter,
hasil_new$iter,
hasil_sec$iter
)
)
perbandingan
## metode akar iterasi
## 1 Bisection 1.52138 33
## 2 Regula Falsi 1.52138 20
## 3 Newton-Raphson 1.52138 4
## 4 Secant 1.52138 8
# Verifikasi dengan fungsi bawaan R
uniroot(f, interval = c(1, 2), tol = 1e-10)$root
## [1] 1.52138
Sub Bab 3 Sistem Persamaan Linear (SPL)
# Membuat matriks koefisien
A <- matrix(
c(4, -1, 0,
-1, 4, -1,
0, -1, 4),
nrow = 3,
byrow = TRUE
)
# Membuat vektor konstanta
b <- c(15, 10, 10)
# Menampilkan matriks dan vektor
print(A)
## [,1] [,2] [,3]
## [1,] 4 -1 0
## [2,] -1 4 -1
## [3,] 0 -1 4
print(b)
## [1] 15 10 10
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
}
# Menjalankan eliminasi Gauss
# Hasil eliminasi Gauss
hasil_gauss <- gauss(A, b)
print(hasil_gauss)
## [1] 4.910714 4.642857 3.660714
# Hasil menggunakan fungsi solve()
hasil_solve <- solve(A, b)
print(hasil_solve)
## [1] 4.910714 4.642857 3.660714
# Dekomposisi LU
lu_res <- expand(lu(Matrix(A)))
lu_res$L
## 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
lu_res$U
## 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
# Dekomposisi Cholesky
chol(A)
## [,1] [,2] [,3]
## [1,] 2 -0.500000 0.0000000
## [2,] 0 1.936492 -0.5163978
## [3,] 0 0.000000 1.9321836
# Penyelesaian menggunakan dekomposisi QR
qr.solve(A, b)
## [1] 4.910714 4.642857 3.660714
# metode jacobi
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 <- as.vector((b - R %*% x) / D)
if (max(abs(xn - x)) < tol) {
return(list(x = xn, iter = k))
}
x <- xn
}
list(x = as.vector(x), iter = maxit)
}
# Menjalankan metode Jacobi
jacobi(A, b)
## $x
## [1] 4.910714 4.642857 3.660714
##
## $iter
## [1] 21
# metode gauss seidel
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)
}
# Menjalankan metode Gauss-Seidel
gauss_seidel(A, b)
## $x
## [1] 4.910714 4.642857 3.660714
##
## $iter
## [1] 12
# Matriks mirip Hilbert
H <- 1 / outer(1:6, 1:6, "+")
# Bilangan kondisi matriks
kappa(H)
## [1] 61143998
kappa(A)
## [1] 2.308061
Sub Bab 4 Interpolasi dan Pencocokan Kurva
# Data yang diketahui
x <- c(0, 1, 2, 3)
y <- c(1, 2, 5, 10)
# Fungsi 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]))
}))
})
}
# Menghitung nilai interpolasi
lagrange(x, y, c(0.5, 1.5, 2.5))
## [1] 1.25 3.25 7.25
# Menghitung koefisien Newton
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
}
# Mengevaluasi polinom Newton
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
})
}
# Menghitung koefisien dan hasil interpolasi
cf <- newton_dd(x, y)
print(cf)
## [1] 1 1 1 0
newton_eval(x, cf, c(0.5, 1.5, 2.5))
## [1] 1.25 3.25 7.25
# Membuat data
xs <- seq(-5, 5, length.out = 11)
ys <- 1 / (1 + xs^2)
# Membuat titik untuk kurva yang lebih halus
xx <- seq(-5, 5, length.out = 200)
# Menggambar fungsi asli
plot(
xx, 1 / (1 + xx^2),
type = "l",
lwd = 2,
main = "Perbandingan Fungsi, Lagrange, dan Spline",
xlab = "x",
ylab = "y"
)
# Menambahkan polinom Lagrange
lines(xx, lagrange(xs, ys, xx), col = "red")
# Menambahkan spline kubik
lines(spline(xs, ys, n = 200), col = "blue", lwd = 2)
# Menampilkan titik data
points(xs, ys, pch = 19)
# Menambahkan legenda
legend(
"topright",
legend = c("Fungsi asli", "Lagrange", "Spline"),
col = c("black", "red", "blue"),
lty = 1,
lwd = 2
)
approx(xs, ys, xout = 0.5) # interpolasi linear
## $x
## [1] 0.5
##
## $y
## [1] 0.75
# Membuat data simulasi
set.seed(1)
x <- 1:20
y <- 3 + 2 * x + 0.1 * x^2 + rnorm(20, sd = 2)
# Membentuk matriks desain
X <- cbind(1, x, x^2)
# Menghitung koefisien dengan persamaan normal
beta <- solve(t(X) %*% X, t(X) %*% y)
beta
## [,1]
## 3.0382456
## x 2.0130475
## 0.1014342
# Membandingkan dengan fungsi lm()
coef(lm(y ~ x + I(x^2)))
## (Intercept) x I(x^2)
## 3.0382456 2.0130475 0.1014342
# Regresi nonlinear: y = a * exp(b * x)
# Membuat data simulasi
xd <- 1:10
yd <- 2 * exp(0.3 * xd) * (1 + rnorm(10, sd = 0.02))
# Estimasi model nonlinear
fit <- nls(
yd ~ a * exp(b * xd),
start = list(a = 1, b = 0.1)
)
# Melihat hasil estimasi
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