Numerical Stability

Dua algoritma yang setara secara matematis bisa memberikan hasil yang sangat berbeda di komputer. Penyebabnya bukan bug, melainkan cara algoritma memperlakukan galat pembulatan floating point: algoritma stabil menjaga galat tetap kecil, algoritma tidak stabil membesarkannya — kadang ribuan kali lipat. Chapter ini membedah empat kasus klasik: penjumlahan, sistem persamaan linear, iterasi rekursif, dan eliminasi Gauss, lalu menutup dengan teknik mitigasi yang bisa langsung dipraktikkan.

Algoritma Stabil vs Tidak Stabil #

Definisi kerjanya: algoritma stabil jika galat hasil akhir sebanding dengan galat pembulatan input — tidak ada pembesaran. Contoh termudah adalah penjumlahan. Loop penjumlahan biasa (naive summation) menambahkan elemen satu per satu; ketika partial sum sudah jauh lebih besar dari elemen berikutnya, elemen kecil itu “tenggelam” karena pembulatan. Compensated summation (Kahan) menyimpan galat pembulatan yang hilang pada variabel koreksi c, lalu mengembalikannya pada penjumlahan berikutnya.

Data ujinya ekstrem tapi realistis: satu nilai 1e9 dan 1000 nilai 1e-9. Jumlah eksaknya 1000000000.000001:

w = [1e9, repmat(1e-9, 1, 1000)];
sNaive = 0;
for k = 1:numel(w)
    sNaive = sNaive + w(k);              % nilai kecil tenggelam
end
sKahan = 0;  c = 0;
for k = 1:numel(w)
    y = w(k) - c;
    t = sKahan + y;
    c = (t - sKahan) - y;                % catat galat pembulatan yang hilang
    sKahan = t;
end
fprintf('jumlah eksak : %.6f\n', 1e9 + 1000*1e-9);
fprintf('naive         : %.6f\n', sNaive);
fprintf('kahan         : %.6f\n', sKahan);

Output yang Dihasilkan:

jumlah eksak : 1000000000.000001
naive         : 1000000000.000000
kahan         : 1000000000.000001

Loop naif kehilangan seluruh kontribusi 1000 nilai kecil (hasil 1000000000.000000), sementara Kahan memulihkannya utuh. Kedua algoritma menjumlahkan bilangan yang sama — hanya urutan dan pencatatan galatnya yang berbeda.

Kondisi Matriks: Masalah yang Sensitif terhadap Gangguan #

Kadang masalahnya bukan algoritmanya, melainkan soalnya sendiri. Angka kondisi cond(A) mengukur sensitivitas solusi sistem A*x = b terhadap gangguan: semakin dekat A ke singular, semakin besar pembesaran galat. Aturan kasarnya: galat relatif x bisa mencapai cond(A) kali galat relatif b.

Matriks hampir singular A = [1 1; 1 1.00001] punya cond(A) = 400002. Gangguan 1e-4 pada b cukup untuk membalikkan solusi:

A  = [1 1; 1 1.00001];                   % hampir singular
b  = [3; 3.00001];
bp = [3.0001; 3.00002];                  % gangguan sekitar 1e-4
x1 = A\b;   x2 = A\bp;
B  = [2 1; 1 3];                         % well-conditioned
y1 = B\b;   y2 = B\bp;
fprintf('kond buruk: x awal = %.1f, %.1f\n', x1(1), x1(2));
fprintf('           x baru = %.4f, %.1f\n', x2(1), x2(2));
fprintf('           cond(A) = %g, residu baru = %.2e\n', ...
        cond(A), norm(A*x2 - bp));
fprintf('kond baik : x awal = %.6f, %.6f\n', y1(1), y1(2));
fprintf('           x baru = %.6f, %.6f\n', y2(1), y2(2));
fprintf('           cond(B) = %g\n', cond(B));

Output yang Dihasilkan:

kond buruk: x awal = 2.0, 1.0
           x baru = 11.0001, -8.0
           cond(A) = 400002, residu baru = 6.28e-16
kond baik : x awal = 1.199998, 0.600004
           x baru = 1.200056, 0.599988
           cond(B) = 2.61803

Gangguan mini menggeser solusi sistem buruk secara dramatis, sementara sistem sehat hanya bergeser halus. Perhatikan juga residu sistem buruk: 6.28e-16, nyaris nol. Solusi yang salah total tetap menghasilkan residu kecil — residu kecil bukan jaminan solusi akurat; hanya kondisi matriks yang menjelaskan selisihnya.

Pertumbuhan Galat pada Iterasi: Fibonacci #

Rekursi menghadapi bahaya tambahan: galat satu langkah menjadi input langkah berikutnya. Fibonacci punya dua rumus: iteratif (f_n = f_n-1 + f_n-2) dan rumus Binet (f_n = phi^n/sqrt(5) dengan phi = (1+sqrt(5))/2). Keduanya tepat secara matematis; secara numerik hanya satu yang bisa dipercaya.

Phi tersimpan sebagai bilangan double dengan galat pembulatan kecil. Pada rumus Binet, pangkat n membesarkan galat relatif phi kira-kira n kali — makin besar n, makin liar hasilnya. Iterasi penjumlahan hanya menambahkan galat pembulatan lokal di tiap langkah dan tetap terkendali:

fibs = zeros(1, 80);
fibs(1) = 1;  fibs(2) = 1;
for k = 3:80
    fibs(k) = fibs(k-1) + fibs(k-2);     % iteratif: pembulatan terkendali
end
phi = (1 + sqrt(5))/2;
for n = 71:2:77
    fibBinet = round(phi^n/sqrt(5));     % Binet: galat dibesarkan pangkat
    fprintf('n=%d: iteratif=%.0f  binet=%.0f  galat=%.0f\n', ...
            n, fibs(n), fibBinet, fibBinet - fibs(n));
end

Output yang Dihasilkan:

n=71: iteratif=308061521170129  binet=308061521170130  galat=1
n=73: iteratif=806515533049393  binet=806515533049395  galat=2
n=75: iteratif=2111485077978050  binet=2111485077978055  galat=5
n=77: iteratif=5527939700884757  binet=5527939700884771  galat=14

Binet masih akurat hingga n = 70, lalu mulai melenceng: galat tumbuh dari 1 menjadi 14 hanya dalam enam langkah — laju pertumbuhannya sebanding dengan angka Fibonacci itu sendiri. Iterasi penjumlahan tetap eksak di rentang ini. Pilih rumus yang galatnya tidak diperbesar struktur perhitungannya.

Pivoting pada Eliminasi Gauss #

Eliminasi Gauss membagi setiap baris dengan elemen pivot. Ketika pivot sangat kecil, pembagian itu menghasilkan multiplier raksasa yang menelan digit presisi baris di bawahnya. Solusinya pivoting: tukar baris supaya elemen terbesar menjadi pivot — langkah murah yang menyelamatkan akurasi. Operator backslash MATLAB menjalankan eliminasi LU lengkap dengan partial pivoting secara otomatis.

Perbandingannya terlihat pada sistem dengan pivot mungil 1e-20:

A2 = [1e-20 1; 1 1];   b2 = [1; 0];
An = A2;   bn = b2;
m = An(2,1)/An(1,1);                     % multiplier raksasa 1e20
An(2,:) = An(2,:) - m*An(1,:);
bn(2)   = bn(2) - m*bn(1);
x2n = bn(2)/An(2,2);
x1n = (bn(1) - An(1,2)*x2n)/An(1,1);     % bagi pivot mungil: galat meledak
xb = A2\b2;                              % backslash: LU dengan pivoting
fprintf('tanpa pivoting : x = %.6f, %.6f, residu = %.2e\n', ...
        x1n, x2n, norm(A2*[x1n; x2n] - b2));
fprintf('backslash      : x = %.6f, %.6f, residu = %.2e\n', ...
        xb(1), xb(2), norm(A2*xb - b2));

Output yang Dihasilkan:

tanpa pivoting : x = 0.000000, 1.000000, residu = 1.00e+00
backslash      : x = -1.000000, 1.000000, residu = 0.00e+00

Solusi sebenarnya adalah (-1, 1). Eliminasi tanpa pivoting melaporkan (0, 1) dengan residu besar 1.0 — hasilnya jelas rusak; pembagian dengan pivot 1e-20 menghancurkan komponen pertama. Backslash bertukar baris lebih dulu lalu menghitung dengan residu nol persis.

Teknik Mitigasi #

Empat jurus yang berlaku umum:

1. Gunakan builtin: backslash, sum, integral, dan polyfit dirancang dengan strategi numerik yang sudah diuji puluhan tahun — jangan menulis ulang kecuali ada alasan kuat.
2. Scaling: samakan skala satuan (kilometer vs meter, 1e9 vs 1e-9) sebelum berhitung; pembulatan diperlakukan adil antar besaran.
3. Formulasi alternatif: Kahan untuk penjumlahan, iteratif untuk Fibonacci, log-domain untuk perkalian rantai panjang — cari rumus yang tidak membesarkan galat.
4. Presisi lebih tinggi: vpa (Symbolic Math Toolbox, MATLAB-only) menghitung dengan digit arbitrer untuk verifikasi; gunakan untuk mengecek, bukan untuk produksi, karena jauh lebih lambat.

Berkas contoh.m pada bagian bawah halaman menjalankan keempat eksperimen chapter ini berurutan: Kahan, kondisi dua matriks, tabel galat Fibonacci, dan duet eliminasi tanpa pivoting versus backslash.

Ringkasan #

Stabilitas numerik menentukan apakah galat pembulatan diperbesar atau dikendalikan. Penjumlahan Kahan menyelamatkan nilai kecil dari partial sum raksasa; kondisi matriks cond(A) memperingatkan sensitivitas solusi terhadap gangguan b — residu kecil saja tidak cukup sebagai bukti kebenaran. Iterasi yang salah formulasi (Binet) menumbuhkan galat secara eksponensial sementara iterasi stabil tetap eksak. Backslash menyelesaikan sistem dengan pivoting otomatis; mitigasi umum: percayakan komputasi pada builtin, samakan skala, pilih formulasi stabil, dan verifikasi dengan presisi tinggi.