Pertemuan 3 · Analisis Numerik · Program Studi Statistika · Universitas Matana
Implementasikan fungsi yang mengubah bilangan desimal menjadi tiga komponen 64-bit (sign, exponent field, mantissa) memakai algoritma konversi biner yang sudah dikuasai.
encode64 <- function(d) {
sign <- ifelse(d < 0, 1, 0)
d <- abs(d)
exponent <- floor(log2(d)) # pangkat pada bentuk standar
# 1.b1b2... x 2^exponent
remainder <- d - 2^exponent # sisa yang akan diuraikan jadi 52 bit mantissa
bits <- character(52)
for (i in 1:52) {
p <- exponent - i # posisi bit ke-i bernilai 2^p
if (remainder >= 2^p) {
bits[i] <- '1'
remainder <- remainder - 2^p
} else {
bits[i] <- '0'
}
}
list(
sign = sign,
exponent_field = exponent + 1023,
mantissa = paste0(bits, collapse = '')
)
}
encode64(5.0) # bandingkan dengan Contoh Dasar di kelas## $sign
## [1] 0
##
## $exponent_field
## [1] 1025
##
## $mantissa
## [1] "0100000000000000000000000000000000000000000000000000"encode64(11.5625) # bandingkan dengan Contoh Menengah di kelas## $sign
## [1] 0
##
## $exponent_field
## [1] 1026
##
## $mantissa
## [1] "0111001000000000000000000000000000000000000000000000"encode64(2.5) # coba bilangan barumu sendiri## $sign
## [1] 0
##
## $exponent_field
## [1] 1024
##
## $mantissa
## [1] "0100000000000000000000000000000000000000000000000000"Apakah exponent_field dari
encode64(5.0) dan encode64(11.5625) sesuai
dengan hitungan manual di kelas (1025 dan 1026)?
Ya. Untuk 5.0:
\[ 5=(101)_2=1.01\times2^2 \]
Sehingga:
\[ exponent\ field=2+1023=1025 \]
Untuk 11.5625:
\[ 11.5625=(1011.1001)_2=1.0111001\times2^3 \]
Sehingga:
\[ exponent\ field=3+1023=1026 \]
Jadi hasilnya sesuai dengan hitungan manual di kelas.
Untuk encode64(2.5), tuliskan bentuk standar
binernya secara manual (1.b1b2... × 2^exponent), lalu
cocokkan dengan output fungsi.
\[ 2.5=(10.1)_2 \]
Bentuk standarnya:
\[ 2.5=1.01\times2^1 \]
Jadi exponent = 1 dan:
\[ exponent\ field=1+1023=1024 \]
Mantissanya adalah:
0100000000000000000000000000000000000000000000000000Hasil tersebut sesuai dengan output fungsi
encode64(2.5).
R menyimpan machine epsilon di .Machine$double.eps. Kita
buktikan langsung bahwa gap antar bilangan floating-point
melebar seiring membesarnya bilangan.
eps <- .Machine$double.eps
gap_di <- function(exponent) {
eps * 2^exponent
}
gap0 <- gap_di(0) # gap di sekitar 1
gap10 <- gap_di(10) # gap di sekitar 1024
c(
gap0 = gap0,
gap10 = gap10,
rasio = gap10 / gap0
)## gap0 gap10 rasio
## 2.220446e-16 2.273737e-13 1.024000e+03# Uji langsung: apakah menambahkan SETENGAH gap ke 1024
# cukup untuk mengubah nilainya di floating-point?
(2^10 + gap10 / 2) == 2^10 # apa dugaanmu?## [1] TRUE(2^10 + gap10) == 2^10 # bandingkan dengan baris di atas## [1] FALSEBerapa kali lipat gap10 dibanding
gap0? Apakah cocok dengan hitungan ‘Contoh Menantang’ di
slide (sekitar 1024 kali)?
\[ \frac{gap10}{gap0}=1024 \]
Jadi gap10 adalah 1024 kali lebih besar
dibandingkan gap0.
Hasil ini cocok dengan hitungan pada contoh di kelas. Hal ini menunjukkan bahwa jarak antar bilangan floating-point semakin besar ketika skala bilangannya semakin besar.
Jelaskan hasil dua baris terakhir kode — kenapa salah satu
TRUE dan yang lain FALSE? Kaitkan dengan
definisi machine epsilon sebagai jarak floating-point
terdekat.
Ketika setengah gap ditambahkan ke 1024, hasilnya masih dianggap sama
dengan 1024 dalam floating-point sehingga menghasilkan
TRUE.
Sedangkan ketika satu gap penuh ditambahkan ke 1024, hasilnya sudah
mencapai bilangan floating-point berikutnya. Nilainya berubah sehingga
menghasilkan FALSE.
Hal ini berkaitan dengan machine epsilon sebagai ukuran ketelitian floating-point. Jika perubahan yang diberikan terlalu kecil dibandingkan jarak floating-point, perubahan tersebut tidak selalu mengubah nilai yang tersimpan.
Jalankan ulang dua demo dari kelas, lalu perluas sedikit untuk melihat pola yang lebih lengkap.
h_vals <- 10^-(10:16)
hasil <- 1 + h_vals - 1
data.frame(
h = h_vals,
hasil = hasil,
rel_error = abs(hasil - h_vals) / h_vals
)## h hasil rel_error
## 1 1e-10 1.000000e-10 8.274037e-08
## 2 1e-11 1.000000e-11 8.274037e-08
## 3 1e-12 1.000089e-12 8.890058e-05
## 4 1e-13 9.992007e-14 7.992778e-04
## 5 1e-14 9.992007e-15 7.992778e-04
## 6 1e-15 1.110223e-15 1.102230e-01
## 7 1e-16 0.000000e+00 1.000000e+00y <- sqrt(2)^2 - 2 # secara matematis PERSIS 0 (akar dikuadratkan lagi)
print(y, digits = 17)## [1] 4.4408920985006262e-16y / .Machine$double.eps # y ini berapa kali lipat machine epsilon?## [1] 2Pada tabel (a), di titik h berapa rel_error
mulai melonjak signifikan? Kaitkan dengan machine epsilon dari Aktivitas
2.
rel_error mulai meningkat cukup jelas ketika
h mendekati 10^-15 dan menjadi sangat besar
pada h = 10^-16.
Hal ini berkaitan dengan machine epsilon yang bernilai sekitar
2.22 × 10^-16. Ketika h sangat kecil dan
mendekati batas ketelitian floating-point, penambahan h ke
1 dapat menyebabkan informasi tentang h tidak tersimpan
dengan akurat.
Pada (b), y seharusnya TEPAT 0 secara matematis tapi tidak di R. Dari mana asal ‘digit sampah’ itu, dan kenapa besarnya sudah dalam orde machine epsilon (lihat hasil baris terakhir)? Jelaskan kaitannya dengan Subtracting Nearly Equal.
Secara matematis:
\[ (\sqrt{2})^2-2=0 \]
Namun komputer tidak dapat menyimpan nilai sqrt(2)
secara tepat dalam floating-point karena jumlah digitnya terbatas.
Ketika nilai pendekatan sqrt(2) tersebut dikuadratkan,
hasilnya sedikit berbeda dari 2. Saat dikurangi dengan 2, muncul sisa
kecil atau “digit sampah”.
Nilai tersebut berada pada orde machine epsilon karena kesalahan pembulatannya sangat kecil dan berada di sekitar batas ketelitian bilangan floating-point.
Hal ini merupakan contoh Subtracting Nearly Equal, yaitu ketika dua bilangan yang nilainya sangat dekat dikurangkan, sehingga kesalahan pembulatan dapat menjadi lebih terlihat.
Implementasikan Eliminasi Gauss TANPA pivoting (persis alur ‘apa adanya’ dari kelas), lalu uji pada dua sistem: yang mulus dan yang bermasalah.
gauss_naive <- function(A, b) {
aug <- cbind(A, b)
n <- nrow(aug)
for (k in 1:(n - 1)) {
for (i in (k + 1):n) {
faktor <- aug[i, k] / aug[k, k]
aug[i, ] <- aug[i, ] - faktor * aug[k, ]
}
}
x <- numeric(n)
x[n] <- aug[n, n + 1] / aug[n, n]
for (i in (n - 1):1) {
x[i] <- (aug[i, n + 1] -
sum(aug[i, (i + 1):n] * x[(i + 1):n])) /
aug[i, i]
}
x
}Jawaban seharusnya:
x = 2, y = 3, z = -1
A1 <- matrix(
c(2, 1, -1,
-3, -1, 2,
-2, 1, 2),
nrow = 3,
byrow = TRUE
)
b1 <- c(8, -11, -3)
gauss_naive(A1, b1)## [1] 2 3 -1gauss_naive(A1, b1) menghasilkan
x=2, y=3, z=-1 sesuai kelas?
Ya. Fungsi gauss_naive(A1, b1) menghasilkan:
x = 2
y = 3
z = -1Hasil tersebut sesuai dengan jawaban yang diberikan pada contoh di kelas.
Jawaban seharusnya:
x = 1, y = 2, z = 3
A2 <- matrix(
c(1, 1, 1,
2, 2, 5,
1, -1, 1),
nrow = 3,
byrow = TRUE
)
b2 <- c(6, 21, 2)
gauss_naive(A2, b2) # perhatikan baik-baik apa yang terjadi!## [1] NaN NaN NaNsolve(A2, b2) # bandingkan dengan fungsi bawaan R## [1] 1 2 3Setelah eliminasi kolom x, diperoleh:
\[ \begin{bmatrix} 1 & 1 & 1 & | & 6\\ 0 & 0 & 3 & | & 9\\ 0 & -2 & 0 & | & -4 \end{bmatrix} \]
Pada baris kedua, pivot untuk kolom y bernilai 0.
gauss_naive(A2, b2) semestinya bermasalah (hasil
NaN/Inf atau salah). Jelaskan MENGAPA — kaitkan dengan
‘Contoh Menantang: Ketika Pivot Bernilai 0’ yang dibahas di
kelas.
Setelah eliminasi kolom x, diperoleh:
[1 1 1 | 6]
[0 0 3 | 9]
[0 -2 0 | -4]Pada baris kedua, elemen pivot untuk kolom y bernilai
0.
Fungsi gauss_naive tidak melakukan pertukaran baris atau
pivoting, sehingga tetap menggunakan angka 0 sebagai pivot.
Akibatnya terjadi pembagian dengan 0 ketika menghitung faktor
eliminasi. Hal tersebut dapat menghasilkan Inf,
NaN, atau hasil yang salah.
Jadi masalahnya bukan karena sistem persamaannya tidak memiliki solusi, tetapi karena algoritma eliminasi Gauss tanpa pivoting tidak mampu menangani pivot yang bernilai 0.
solve(A2, b2) tetap berhasil memberi jawaban
yang benar walau gauss_naive gagal. Menurutmu, apa yang
dilakukan solve() secara internal sehingga tidak mengalami
masalah yang sama? (Tidak perlu jawaban teknis lengkap — cukup intuisi
awal; akan dibahas tuntas sebagai Partial Pivoting minggu
depan.)
Fungsi solve() menggunakan metode numerik yang lebih
aman dibandingkan gauss_naive sederhana.
Secara intuisi, solve() dapat mengatur atau menukar
baris sehingga tidak menggunakan pivot yang bernilai 0.
Dengan begitu, masalah pembagian dengan pivot 0 dapat dihindari dan solusi yang benar dapat diperoleh.
Hal ini berkaitan dengan konsep Partial Pivoting yang akan dibahas lebih lanjut pada pertemuan berikutnya.
Praktikum ini menunjukkan bahwa representasi floating-point memiliki keterbatasan ketelitian. Machine epsilon dan gap antar bilangan floating-point bergantung pada skala bilangan.
Selain itu, eliminasi Gauss tanpa pivoting dapat mengalami masalah ketika pivot bernilai 0. Strategi pivoting dapat digunakan untuk menghindari masalah tersebut dan memperoleh solusi yang lebih stabil.