Numerical Differentiation

Turunan mengukur laju perubahan — kecepatan dari posisi, kemiringan kurva, medan gaya dari energi. Sayangnya, dua situasi umum membuat turunan simbolik tidak tersedia: data eksperimen hanya berupa tabel titik-titik, dan fungsi yang kita hadapi kadang berupa black box tanpa bentuk analitik. Di situ turunan numerik mengambil alih: hampirkan turunan dari nilai fungsi di beberapa titik berdekatan. Chapter ini membangun tiga rumus beda hingga, menyelidiki jarak h yang optimal (terlalu besar dan terlalu kecil sama-sama berbahaya), lalu memakai fungsi gradient() untuk turunan data pada grid satu dan dua dimensi.

Tiga Rumus Beda Hingga #

Aturan dasarnya: turunan f'(x) adalah kemiringan garis singgung, dan kemiringan bisa dihampiri dengan kemiringan garis potong (secant) lewat dua titik berdekatan. Tiga pilihan titik yang lazim:

Mekanisme di balik layar: garis potong maju dan mundur condong mengikuti kelengkungan f, sedangkan central difference menarik garis yang “mengapit” titik x sehingga kelengkungan kiri-kanan saling meniadakan — itulah sebabnya galatnya turun satu orde. Notasi O(h) dan O(h^2) dibaca: bila h diperkecil 10 kali, galat O(h) mengecil 10 kali, galat O(h^2) mengecil 100 kali.

f = @(x) sin(x);
x0 = 1;
h = 0.01;
eksak = cos(x0);
fwd = (f(x0 + h) - f(x0)) / h;
bwd = (f(x0) - f(x0 - h)) / h;
ctr = (f(x0 + h) - f(x0 - h)) / (2*h);
fprintf('eksak    = %.8f\n', eksak);
fprintf('forward  = %.8f, galat = %.2e\n', fwd, abs(fwd - eksak));
fprintf('backward = %.8f, galat = %.2e\n', bwd, abs(bwd - eksak));
fprintf('central  = %.8f, galat = %.2e\n', ctr, abs(ctr - eksak));

Output yang Dihasilkan:

eksak    = 0.54030231
forward  = 0.53608598, galat = 4.22e-03
backward = 0.54450062, galat = 4.20e-03
central  = 0.54029330, galat = 9.00e-06

Mari kita bedah baris-baris kode di atas:

1. Baris 1-3 (f, x0, h): fungsi uji sin, titik evaluasi x = 1, dan jarak h = 0.01 — satuan-satuan yang akan kita variasikan sepanjang chapter.
2. Baris 4 (eksak = cos(x0)): turunan analitik sin'(x) = cos(x) dijadikan pembanding; di titik ini nilainya 0.54030231.
3. Baris 5-6 (fwd, bwd): dua rumus satu sisi memiringkan garis potong ke depan dan ke belakang; galatnya sama besar (4.22e-03 vs 4.20e-03) karena keduanya ditumpangi galat O(h) yang sama.
4. Baris 7 (ctr): central difference mengapit titik x secara simetris — galatnya 9.00e-06, hampir 500 kali lebih kecil, bukti langsung beda orde O(h) vs O(h^2).

Galat O(h^2): Setiap h Dibelah Dua #

Sifat O(h^2) paling jelas terlihat bila h dibelah dua berulang kali — galat central harus ikut mengecil menjadi seperempatnya:

for h = [0.1 0.05 0.025]
    d = (sin(x0 + h) - sin(x0 - h)) / (2*h);
    fprintf('h = %.3f, galat = %.3e\n', h, abs(d - eksak));
end

Output yang Dihasilkan:

h = 0.100, galat = 9.001e-04
h = 0.050, galat = 2.251e-04
h = 0.025, galat = 5.628e-05

Setiap pembelahan h (0.1 menjadi 0.05 menjadi 0.025) mengecilkan galat tepat menjadi seperempatnya: 9.001e-04 -> 2.251e-04 -> 5.628e-05. Pola ini adalah sidik jari O(h^2) — perbaiki h dengan faktor k, galat membaik dengan faktor k^2.

Trade-off h: Truncation vs Round-off #

Godaan berikutnya: kalau h makin kecil makin akurat, kenapa tidak langsung pakai h = 1e-10? Karena dua galat bertarung di dalam rumus central. Galat truncation sebesar O(h^2) mengecil saat h mengecil, tetapi pengurangan f(x+h) - f(x-h) membuang digit signifikan saat h terlalu kecil (pola cancellation dari chapter Numerical Errors), sehingga galat round-off sekitar eps/h justru MEMBESAR:

fprintf('%8s %12s\n', 'h', 'galat');
for h = [1e-1 1e-2 1e-3 1e-4 1e-5 1e-6 1e-7 1e-8 1e-9 1e-10]
    d = (sin(x0 + h) - sin(x0 - h)) / (2*h);
    fprintf('%8.0e %12.2e\n', h, abs(d - eksak));
end

Output yang Dihasilkan:

       h        galat
   1e-01     9.00e-04
   1e-02     9.00e-06
   1e-03     9.01e-08
   1e-04     9.00e-10
   1e-05     1.11e-11
   1e-06     2.77e-11
   1e-07     1.94e-10
   1e-08     2.58e-09
   1e-09     2.97e-09
   1e-10     5.85e-08

Tabel di atas membentuk kurva berbentuk mangkuk. Dari h = 1e-1 sampai h = 1e-5, galat menurun rapi mengikuti O(h^2) — wilayah truncation. Di bawah h = 1e-5, galat berbalik MEMBESAR — wilayah round-off yang dipicu pembatalan pada pengurangan. Titik terbaiknya sekitar h = 1e-5 dengan galat 1.11e-11, dan teori memang meramalkan h optimal di sekitar eps^(1/3) ≈ 6e-6 untuk central difference.

gradient: Turunan Data pada Grid #

Untuk data yang sudah tersimpan sebagai vektor nilai pada grid berjarak tetap, MATLAB menyediakan gradient(). Mekanisme di balik layar: fungsi ini memakai central difference di titik interior dan rumus satu sisi (forward/backward) di dua titik tepi:

x = 0:0.1:pi;
y = sin(x);
g = gradient(y, 0.1);
fprintf('g(11)  (x = 1.0, interior) = %.6f, cos(1) = %.6f\n', g(11), cos(1));
fprintf('g(1)   (x = 0.0, tepi)     = %.6f, cos(0) = %.6f\n', g(1), cos(0));

Output yang Dihasilkan:

g(11)  (x = 1.0, interior) = 0.539402, cos(1) = 0.540302
g(1)   (x = 0.0, tepi)     = 0.998334, cos(0) = 1.000000

Mari kita bedah baris-baris kode di atas:

1. Baris 1-2 (x = 0:0.1:pi; y = sin(x)): grid berjarak tetap 0.1 beserta nilainya — bentuk data khas hasil simulasi atau sensor.
2. Baris 3 (g = gradient(y, 0.1)): argumen kedua menyatakan jarak grid h; keluaran g berukuran sama dengan y.
3. Baris 4 (g(11)): titik interior di x = 1.0 memakai central difference dengan h = 0.1 — hasilnya 0.539402, persis seperti rumus central kita di tabel sebelumnya.
4. Baris 5 (g(1)): titik tepi x = 0 tidak punya tetangga kiri, jadi gradient jatuh ke forward difference — galatnya 1.7e-03, lebih besar dari interior.

Turunan pada Grid 2D #

Data dua dimensi — citra, medan suhu, permukaan ketinggian — punya dua turunan arah. Fungsi gradient() menerima matriks dan mengembalikan dua komponen sekaligus: Gx (turunan arah kolom/x) dan Gy (turunan arah baris/y):

[X, Y] = meshgrid(0:0.5:2, 0:0.5:2);
Z = X.^2 + Y.^2;
[Gx, Gy] = gradient(Z, 0.5, 0.5);
fprintf('Gx(4,4) = %.4f (dZ/dX eksak = 2*x = 3)\n', Gx(4,4));
fprintf('Gy(4,4) = %.4f (dZ/dY eksak = 2*y = 3)\n', Gy(4,4));

Output yang Dihasilkan:

Gx(4,4) = 3.0000 (dZ/dX eksak = 2*x = 3)
Gy(4,4) = 3.0000 (dZ/dY eksak = 2*y = 3)

Mari kita bedah baris-baris kode di atas:

1. Baris 1-2 (meshgrid, Z): meshgrid membangun grid koordinat; permukaan Z = X.^2 + Y.^2 dipilih karena turunannya eksak diketahui (2x dan 2y) sehingga galat bisa diperiksa.
2. Baris 3 ([Gx, Gy] = gradient(Z, 0.5, 0.5)): dua argumen jarak — 0.5 untuk arah x (antar kolom) dan 0.5 untuk arah y (antar baris); dua matriks keluaran berukuran sama dengan Z.
3. Baris 4-5 (Gx(4,4), Gy(4,4)): pada titik grid keempat-baris-keempat-kolom (x = 1.5, y = 1.5), kedua komponen menghasilkan 3.0000 — tepat sama dengan turunan analitik karena fungsi kuadratik polos diparameterisasi grid berjarak tetap.

Berkas contoh.m pada bagian bawah halaman menjalankan seluruh contoh di atas — tiga rumus beda, tabel trade-off h, dan turunan grid 1D/2D.

Ringkasan #

Turunan numerik dihampiri lewat beda hingga: forward/backward dengan galat O(h) dan central difference dengan galat O(h^2) — pilihan default kita. Pemilihan h adalah tarik ulur dua galat: truncation O(h^2) vs round-off eps/h, dengan titik manis sekitar h = 1e-5 (kurva galat berbentuk mangkuk). Untuk data grid, gradient(y, h) menerapkan central difference di interior dan rumus satu sisi di tepi, sementara [Gx, Gy] = gradient(Z, hx, hy) memperluasnya ke permukaan dua dimensi. Analisis trade-off yang sama akan kita pakai lagi saat menghampiri integral di chapter berikutnya.