Anum Week 1

Prerequisite

Bab ini akan mencakup beberapa hal di aljabar linear, jadi cobalah untuk mengingat kembali apa itu Operasi Baris Elementer (OBE).

Persamaan Linear yang Mudah

Persamaan linear dalam aljabar linear biasa dituliskan sebagai berikut:

\[A\vec{x}=\vec{b}\]

Solusi yang mudah ditemukan adalah jika bentuknya $U\vec{x}=\vec{b}$ atau $L\vec{x}=\vec{b}$.

Saya harap anda sudah mengetahui itu di aljabar linear jadi saya tidak akan membahasnya terlalu lama.

Untuk mengubah $A\vec{x}=\vec{b}$ menjadi eselon baris, kita biasanya melakukan OBE. Perlu diperhatikan bahwa jika kita melanjutkanya, biasanya kita bisa langsung mendapatkan $\vec{x}$ dari bentuk eselon baris tereduksi. Tetapi, kita ingin komputer untuk melakukan hal ini. Oleh karena itu, dibuatlah algoritma yang dapat mengubah matriks segitiga atas ($U$) atau matriks segitiga bawah ($L$) yang kemudian bisa kita buat menjadi bentuk eselon baris tereduksi.

Berikut adalah algoritma untuk mengubah $U\vec{x}=\vec{b}$. Yaitu mendapatkan SPL yang sudah berbentuk matriks segitiga atas.

function x = Atas(U,b)
  n = length(b);
  x = zeros(n,1);
  for i=n:-1:1
    x(i) = (b(i) - U(i,i+1:n) * x(i+1:n)) / U(i,i);
  end
end

Dan berikut adalah algoritma yang serupa untuk matriks segitiga bawah

function x = Bawah(L,b)
  n = length(b);
  x = zeros(n,1);
  for i=1:n
    x(i) = (b(i)-L(i,1:i-1) * x(1:i-1)) / L(i,i);
  endfor
endfunction

Berikut adalah demo untuk matriks segitiga atas:

A = rand(5)
b = rand(5,1)

L = tril(A)

x = Bawah(L, b)
L*x-b

Harusnya, $L*x-b$ hasinya adalah 0 atau bilangan yang sangat dekat dengan 0 (ingat materi soal mesin epsilon).

Cobalah pikirkan apa itu $L*x-b$. Hint: hasinya mendekati 0.

Selanjutnya, berikut ini adalah demo untuk matriks segitiga bawah:

A = rand(5)
b = rand(5,1)

U = triu(A)

x = Atas(U, b)
U*x-b

Jika kamu mengerti $Lx-b$, seharusnya kamu juga mengerti $Ux-b$.

Oh iya, coba analisis kompleksitas dari kedua algoritma diatas.

Persamaan Linear yang Sebenarnya

Setelah mengetahui bagaimana cara menyelesaikan persamaan linear matsriks segitiga atas maupun matriks segitiga bawah, kita kemudian perlu mengetahui bagaiman acara untuk menyelesaikan matriks yang sebenarnya bukan? Yaitu bentuk matriks

\[A\vec{x}=\vec{b}\]

Untuk menyelesaikan persoalan tersebut, mari kita ingat kembali materi gauss jordan untuk menyelesaikan persamaan matriks. Ingat bahwa pada gauss jordal, kita membuat $A\vec{x}=\vec{b}$ menjadi matriks segitiga atas atau bahwah yang kemudian kita solve. Ingat juga bahwa kita menggunakan OBE untuk melakukannnya. Berikut adalah algoritma gauss jordan yang mengubah A dan b (matriks augmented) menjadi matriks segitiga atas atau matriks segitiga bawah.

function [U, b1] = GaussOnly(A, b)
  n = length(b);
  C = [A b];
  %kolom (paling kiri sampai sebelum habis)
  for j=1:n-1
    %baris (kedua paling atas sampai bawah)
    for i=j+1:n
      %ratio (kiri (dia menurun) / paling atas kiri)
      %notice bahwa dia dibagi elemen diagonal ( C(j,j) )
      m = C(i, j)/C(j,j);
      %update
      C(i,:) = C(i,:)-m*C(j,:);
    end
  end
  U = C(:,1:n)
  b1 = C(:,n+1)
end

Perhatikan bahwa fungsi diatas mengembalikan U dan b1, dimana U adalah matriks segitiga atas dan b1 adalah b yang telah terkena treatment OBE.

Selanjutnya, karena kita sedang berada di kelas analisis numerik, kita akan coba untuk lihat edge case. Seperti yang telah dijelaskan pada blog sebelumnya, bahwa karena mesin tempat kita bekerja memiliki jumlah bit yang finite, kita perlu memerhatikan jika selisih 2 pengurangan dapat lebih kecil dari mesin epsilon atau pembagian dengan 0. Hal ini dapat terjadi jika kita memiliki persamaan berikut:

\[\begin{cases} 0x+y=1 \\ 3x+2y=5 \end{cases}\]

Cobalah kamu melakukan eliminasi gauss secara manual dengan tidak melakukan operasi OBE tukar baris. Hal ini jelas mustahil bukan? Hal tersebut mustahil untuk dilakukan karena pembagian dengan 0.

Sayangnya, algoritma yang kita buat diatas (yaitu GaussOnly) hanya menggunakan OBE dimana kita mengalikan sebuah skalar dan menjumlahkan 2 persamaan (baris).

Karena itu, kita ingin membuat sebuah algoritma yang dapat melakukan OBE tukar baris. Berikut adalah algoritmanya:

function [U, bt] = GaussPivot(A, b)
  % A non singular juga bisa
  % Complexity sama (tapi agak lama dikit)

    % Augment b in A
    A = [A, b];
    n = size(A, 1);

    for k = 1:n-1
        % Find the row index where the maximum absolute value resides
        [~, t] = max(abs(A(k:n, k)));
        t = t + k - 1; % Index adjustment as the column decreases

        % Swap rows k and t
        A([k, t], :) = A([t, k], :);

        for i = k+1:n
            m = A(i, k) / A(k, k);
            % Elimination of the i-th row
            A(i, :) = A(i, :) - m * A(k, :);
        end
    end

    % Extract U (upper triangular matrix) and bt (augmented b)
    U = A(1:n, 1:n);
    bt = A(:, n+1);
end

Oh ya, jika kamu perhatikan, sebenarnya kita melakukan beberapa kali pembagian. Hal menarik dari pembagian ini adalah hasilnya sebenarnya bisa dibuat menjadi matriks dimana matriks itu adalah salah satu hasil dekomposisi dari matriks awal yang kita miliki. Hal ini dinamakan LU decomposition.

Secara sederhana, dekomposisi ini bisa membuat matriks A yang kita miliki menjadi matriks L dan matriks U. Saya harap Anda masih mengingat matriks tersebut, yaitu matriks segitiga atas dan segitiga bawah. Jika kita memiliki matriks tersebut, kita bisa menyelesaikan persamaan linear seperti ini:

\[\begin{align*} A\vec{x}&=\vec{b} \\ LU\vec{x}&=\vec{b} \end{align*}\]

Perhatikan bahwa kita bisa membuat permisalan $U\vec{x}=\vec{y}$.

Hal tersebut kemudian bisa kita subtitusi pada persamaan diatas yang akan menyebabkan:

\[L\vec{y}=\vec{b}\]

yang kemudian bisa kita cari x nya dari:

\[U\vec{x}=\vec{y}\]

Berikut adalah algoritma dekomposisi LU.

function [L, U] = LU_Fac(A)

  [n, n] = size(A);

  L = eye(n);

  for j=1:n-1;
    for i=j+1:n;
      L(i,j) = A(i,j)/A(j,j);
      A(i,:) = A(i,:) - L(i,j)*A(j,:);
    endfor
  endfor
  U = A;
endfunction

Dengan begitu, kamu bisa menyelesaikan sistem persamaan linear dengan cara berikut (dengan asumsi Anda sudah menyimpan algoritma Atas, Bawah, dan LU_Fac):

A = rand(5)
b = rand(5,1)

[L, U] = LU_Fac(A)

%selesaikan Ly=b
y = Bawah(L, b)

%selesaikan Ux=y
x = Atas(U, y)

%Verify answer
A*x-b

Selamat, Anda sudah mengetahui beberapa algoritma untuk menyelesaikan SPL. Pertanyaannya adalah bagaimana kompleksitas dan error propagation pada setiap algoritma yang ada?