LAMPIRAN
21
Lampiran 1 Pembuktian Teorema 2.4.15 (Rumus Proyeksi Ortogonal) Misalkan suatu hiperbidang di ℝ𝑁 mempunyai persamaan 𝐚𝑇 𝐱 = 𝑏 dan misalkan 𝐱 ∗ adalah sebarang titik di ℝ𝑁 , maka proyeksi ortogonal, 𝐱 𝑝 , dari 𝐱 ∗ terhadap hiperbidang tersebut dinyatakan dengan 𝐱𝑝 = 𝐱∗ +
𝑏 − 𝐚𝑇 𝐱 ∗ 𝐚. 𝐚 𝟐
Bukti: Teorema ini dapat dibuktikan dengan membuktikan bahwa (i) titik 𝐱 𝑝 terletak pada hiperbidang 𝐚𝑇 𝐱 = 𝑏 (dalam hal ini, 𝐚𝑇 𝐱 𝑝 = 𝑏) dan (ii) vektor 𝐱 𝑝 − 𝐱 ∗ ortogonal terhadap hiperbidang 𝐚𝑇 𝐱 = 𝑏 (dalam hal ini, vektor 𝐱 𝑝 − 𝐱 ∗ sejajar dengan vektor 𝐚). Pertama, karena 𝐱𝑝 = 𝐱∗ +
𝑏 − 𝐚𝑇 𝐱 ∗ 𝐚, 𝐚 𝟐
maka 𝑏 − 𝐚𝑇 𝐱 ∗ 𝐚 𝐚 𝟐 𝑏 − 𝐚𝑇 𝐱 ∗ 𝑇 = 𝐚𝑇 𝐱 ∗ + 𝐚 𝐚 𝐚 𝟐 𝑇 ∗ 𝑏−𝐚 𝐱 = 𝐚𝑇 𝐱 ∗ + 𝐚 𝟐 𝐚 𝟐 = 𝐚𝑇 𝐱 ∗ + 𝑏 − 𝐚𝑇 𝐱 ∗ = 𝑏.
𝐚𝑇 𝐱 𝑝 = 𝐚𝑇 𝐱 ∗ +
Jadi, terbukti bahwa titik 𝐱 𝑝 terletak pada hiperbidang 𝐚𝑇 𝐱 = 𝑏. Kedua, karena 𝐱𝑝 = 𝐱∗ +
𝑏 − 𝐚𝑇 𝐱 ∗ 𝐚, 𝐚 𝟐
maka 𝑏 − 𝐚𝑇 𝐱 ∗ 𝐚. 𝐚 𝟐 ∗ Karena vektor 𝐱 𝑝 − 𝐱 merupakan merupakan vektor 𝐚 dikalikan dengan suatu skalar, maka vektor 𝐱 𝑝 − 𝐱 ∗ sejajar dengan vektor 𝐚. Jadi, terbukti bahwa vektor 𝐱 𝑝 − 𝐱 ∗ ortogonal terhadap hiperbidang 𝐚𝑇 𝐱 = 𝑏. Bukti lengkap. ∎ 𝐱𝑝 − 𝐱∗ =
22
Lampiran 2 Pembuktian Persamaan (6) Diketahui Persamaan (3), yaitu 𝑓𝑖 𝐱 = 𝐱 +
𝑏𝑖 − 𝐚𝑇𝑖 𝐱 𝐚𝑖 𝐚𝑖 2
untuk 𝑖 = 1, 2, … , 𝑀.
Akan dibuktikan bahwa Persamaan (3) ekivalen dengan Persamaan (6), yaitu 𝑓𝑖 𝐱 = 𝑃𝑖 𝐱 +
𝑏𝑖 𝐚𝑖
2
untuk 𝑖 = 1, 2, … , 𝑀
𝐚𝑖
dengan 𝑃𝑖 = 𝐼 − Bukti: 𝑏𝑖 − 𝐚𝑇𝑖 𝐱 𝐚𝑖 𝐚𝑖 2 𝐚𝑇𝑖 𝐱 𝑏𝑖 =𝐱− 𝐚 + 𝐚 𝐚𝑖 2 𝑖 𝐚𝑖 2 𝑖 1 𝑏𝑖 =𝐱− 𝐚𝑖 𝐚𝑇𝑖 𝐱 + 𝐚 2 𝐚𝑖 𝐚𝑖 2 𝑖 1 𝑏𝑖 = 𝐼𝐱 − 𝐚𝑖 𝐚𝑇𝑖 𝐱 + 𝐚 2 𝐚𝑖 𝐚𝑖 2 𝑖 1 𝑏𝑖 = 𝐼− 𝐚𝑖 𝐚𝑇𝑖 𝐱 + 𝐚. 2 𝐚𝑖 𝐚𝑖 2 𝑖
𝑓𝑖 𝐱 = 𝐱 + ↔ 𝑓𝑖 𝐱 ↔ 𝑓𝑖 𝐱 ↔ 𝑓𝑖 𝐱 ↔ 𝑓𝑖 𝐱
Bukti lengkap. ∎
1 𝐚𝑖
2
𝐚𝑖 𝐚𝑇𝑖 .
23
Lampiran 3 Pembuktian Persamaan (9) Diketahui Persamaan (4), yaitu 𝐹 𝐱; 𝐛 = 𝑓1 ∘ 𝑓2 ∘ ⋯ ∘ 𝑓𝑀 𝐱 = 𝑓1 𝑓2 ⋯ 𝑓𝑀 𝐱
⋯
,
Persamaan (6), yaitu 𝑓𝑖 𝐱 = 𝑃𝑖 𝐱 +
𝑏𝑖 𝐚𝑖
2
untuk 𝑖 = 1, 2, … , 𝑀,
𝐚𝑖
Persamaan (8), yaitu 𝑀
𝑅𝐛 = 𝑖=1
𝑏𝑖 𝐚𝑖
2
𝑄𝑖−1 𝐚𝑖 ,
dan 𝑄𝑖 = 𝑃1 𝑃2 … 𝑃𝑖 (𝑖 = 1, 2, … , 𝑀) dengan 𝑄0 = 𝐼. Akan dibuktikan bahwa Persamaan (4) ekivalen dengan Persamaan (9), yaitu 𝐹 𝐱; 𝐛 = 𝑄𝐱 + 𝑅𝐛 dengan matriks 𝑄 = 𝑄𝑀 . Bukti: 𝐹 𝐱; 𝐛 = 𝑄𝐱 + 𝑅𝐛 𝑀
𝑏𝑖 𝐚𝑖
↔ 𝐹 𝐱; 𝐛 = 𝑃1 𝑃2 … 𝑃𝑀 𝐱 + 𝑖=1
2
𝑄𝑖−1 𝐚𝑖
𝑏1 𝑏2 𝑏𝑀 𝑄0 𝐚1 + 𝑄1 𝐚2 + ⋯ + 𝑄 𝐚 2 2 𝐚1 𝐚2 𝐚𝑀 2 𝑀−1 𝑀 𝑏1 𝑏2 𝑏𝑀 = 𝑃1 𝑃2 … 𝑃𝑀 𝐱 + 𝐚 + 𝑃 𝐚 + ⋯+ 𝑃 𝑃 … 𝑃𝑀−1 𝐚𝑀 𝐚1 2 1 𝐚2 2 1 2 𝐚𝑀 2 1 2 𝑏2 𝑏𝑀 𝑏1 = 𝑃1 𝑃2 … 𝑃𝑀 𝐱 + 𝐚 + ⋯+ 𝑃 … 𝑃𝑀−1 𝐚𝑀 + 𝐚 𝐚2 2 2 𝐚𝑀 2 2 𝐚1 2 1 𝑏2 𝑏𝑀 = 𝑓1 𝑃2 … 𝑃𝑀 𝐱 + 𝐚2 + ⋯ + 𝑃 … 𝑃𝑀−1 𝐚𝑀 2 𝐚2 𝐚𝑀 2 2 𝑏𝑀 𝑏2 = 𝑓1 𝑃2 𝑃3 … 𝑃𝑀 𝐱 + ⋯ + 𝑃3 … 𝑃𝑀−1 𝐚𝑀 + 𝐚 2 𝐚𝑀 𝐚2 2 2 𝑏𝑀 = 𝑓1 𝑓2 𝑃3 … 𝑃𝑀 𝐱 + ⋯ + 𝑃 … 𝑃𝑀−1 𝐚𝑀 𝐚𝑀 2 3
↔ 𝐹 𝐱; 𝐛 = 𝑃1 𝑃2 … 𝑃𝑀 𝐱 + ↔ 𝐹 𝐱; 𝐛 ↔ 𝐹 𝐱; 𝐛 ↔ 𝐹 𝐱; 𝐛 ↔ 𝐹 𝐱; 𝐛 ↔ 𝐹 𝐱; 𝐛 ⋮
↔ 𝐹 𝐱; 𝐛 = 𝑓1 𝑓2 ⋯ 𝑃𝑀 𝐱 + ↔ 𝐹 𝐱; 𝐛 = 𝑓1 𝑓2 ⋯ 𝑓𝑀 𝐱 Bukti lengkap. ∎
𝑏𝑀 𝐚𝑀
⋯
.
2
𝐚𝑀 ⋯
24
Lampiran 4 Bukti Persamaan (10) Berdasarkan Persamaan (5) dan (9), untuk 𝐛 = 𝟎 diperoleh 𝐱
𝑛+1
=𝐹 𝐱
𝑛
; 𝟎 = 𝑄𝐱
𝑛
,
Akan dibuktikan (dengan menggunakan Induksi Matematika) bahwa 𝐱
𝑛
= 𝑄𝑛 𝐱
𝐱
1
0
,
untuk setiap 𝑛 ∈ ℕ. Bukti: Untuk 𝑛 = 1, Ruas kiri: = 𝑄𝐱
0
.
Ruas kanan: 𝑄𝐱 0 . Ruas kiri sama dengan ruas kanan, berarti benar untuk 𝑛 = 1. Anggap benar untuk 𝑛 = 𝑚, 𝑚 ∈ ℕ, yaitu 𝐱
𝑚
= 𝑄𝑚 𝐱
0
.
Akan dibuktikan bahwa benar untuk 𝑛 = 𝑚 + 1, 𝑚 ∈ ℕ, yaitu 𝐱
𝑚 +1
= 𝑄𝑚 +1 𝐱
0
Bukti: 𝐱
𝑚 +1
Terbukti benar untuk 𝑛 = 𝑚 + 1, 𝑚 ∈ ℕ. Bukti lengkap. ∎
= 𝑄𝐱 𝑚 = 𝑄𝑄𝑚 𝐱 0 = 𝑄𝑚 +1 𝐱 0 .
.
25
Lampiran 5 Pembuktian Persamaan (11) Berdasarkan Persamaan (5) dan (9), untuk 𝐛 = 𝟎 diperoleh 𝑛+1
𝐱
𝑛
=𝐹 𝐱
; 𝟎 = 𝑄𝐱
𝑛
,
dan berdasarkan bagian pertama dari Teorema 4.3.6, diperoleh 𝑄 = 𝑃𝒦 + 𝑄 = 𝑃𝒦 + 𝑄𝑃ℐ . Akan dibuktikan (dengan menggunakan Induksi Matematika) bahwa 𝑛
𝐱
0
= 𝑃𝒦 𝐱
+ 𝑄𝑛 𝐱
0
,
∈ Im 𝐴𝑇
∈ Ker 𝐴
untuk setiap 𝑛 ∈ ℕ. Bukti: Untuk 𝑛 = 1, Ruas kiri: 𝐱
1
= 𝑄𝐱
0
.
Ruas kanan: 𝑃𝒦 𝐱
0
+ 𝑄𝐱
0
= 𝑃𝒦 + 𝑄 𝐱
0
= 𝑄𝐱
0
.
∈ Im 𝐴𝑇
∈ Ker 𝐴
Ruas kiri sama dengan ruas kanan, berarti benar untuk 𝑛 = 1. Anggap benar untuk 𝑛 = 𝑚, 𝑚 ∈ ℕ, yaitu 𝐱
𝑚
0
= 𝑃𝒦 𝐱
+ 𝑄𝑚 𝐱
0
.
∈ Im 𝐴𝑇
∈ Ker 𝐴
Akan dibuktikan bahwa benar untuk 𝑛 = 𝑚 + 1, 𝑚 ∈ ℕ, yaitu 𝑚 +1
𝐱
0
= 𝑃𝒦 𝐱
+ 𝑄 𝑚 +1 𝐱
0
.
+ 𝑄𝑃ℐ 𝑃𝒦 𝐱
0
∈ Im 𝐴𝑇
∈ Ker 𝐴
Bukti: 𝐱
𝑚 +1
= 𝑄𝐱
𝑚
= 𝑃𝒦 + 𝑄
0
𝑃𝒦 𝐱
+ 𝑄𝑚 𝐱
∈ Ker 𝐴
= 𝑃𝒦 𝑃𝒦 𝐱
0
𝑚
+ 𝑃𝒦 𝑄 𝐱
= 𝑃𝒦 𝐱 =
∈ Ker 𝐴 𝑃𝒦 𝐱 0 ∈ Ker 𝐴
∈ Im 𝐴𝑇 0
∈ Im 𝐴𝑇 𝑚 +1
∈ Ker 𝐴 0
0
+ 𝟎 + 𝑄𝟎 + 𝑄
∈ Ker 𝐴
𝐱
∈ Im 𝐴𝑇
+ 𝑄 𝑚 +1 𝐱 ∈ Im 𝐴𝑇
Terbukti benar untuk 𝑛 = 𝑚 + 1, 𝑚 ∈ ℕ. Bukti lengkap. ∎
0
.
0
+ 𝑄 𝑄𝑚 𝐱
0
∈ Im 𝐴𝑇
26
Lampiran 6 Pembuktian Persamaan (12) Berdasarkan Persamaan (5) dan (9), diperoleh 𝑛+1
𝐱
=𝐹 𝐱
𝑛
; 𝐛 = 𝑄𝐱
𝑛
+ 𝑅𝐛.
Akan dibuktikan (dengan menggunakan Induksi Matematika) bahwa 𝑛−1
𝐱
𝑛
= 𝑄𝑛 𝐱
0
𝑄𝑘 𝑅
+
𝐛,
𝑘=0
untuk setiap 𝑛 ∈ ℕ. Bukti: Untuk 𝑛 = 1, Ruas kiri: 𝐱
1
0
= 𝑄𝐱
+ 𝑅𝐛.
Ruas kanan: 0
𝑄𝐱
0
𝑄𝑘 𝑅
+
𝐛 = 𝑄𝐱
0
+ 𝑄0 𝑅𝐛 = 𝑄𝐱
0
+ 𝑅𝐛.
𝑘=0
Ruas kiri sama dengan ruas kanan, berarti benar untuk 𝑛 = 1. Anggap benar untuk 𝑛 = 𝑚, 𝑚 ∈ ℕ, yaitu 𝑚 −1
𝐱
𝑚
= 𝑄𝑚 𝐱
0
𝑄𝑘 𝑅
+
𝐛.
𝑘=0
Akan dibuktikan bahwa benar untuk 𝑛 = 𝑚 + 1, 𝑚 ∈ ℕ, yaitu 𝑚
𝐱
𝑚 +1
=𝑄
𝑚 +1
𝐱
0
𝑄𝑘 𝑅
+
𝐛.
𝑘=0
Bukti: 𝐱
𝑚 +1
= 𝑄𝐱
𝑚
+ 𝑅𝐛 𝑚 −1 𝑚
=𝑄 𝑄 𝐱
0
𝑄𝑘 𝑅
+
𝐛 + 𝑅𝐛
𝑘=0 𝑚 −1
= 𝑄𝑚 +1 𝐱
0
𝑄𝑘 𝑅
𝐛 + 𝑅𝐛
𝑄𝑘 +1 𝑅
𝐛 + 𝑅𝐛
+𝑄 𝑘=0 𝑚 −1
= 𝑄𝑚 +1 𝐱
0
+ 𝑘=0 𝑚
= 𝑄𝑚 +1 𝐱
0
+
𝑄𝑘 𝑅
𝐛 + 𝑄0 𝑅𝐛
𝑄𝑘 𝑅
𝐛.
𝑘=1 𝑚
= 𝑄𝑚 +1 𝐱
0
+ 𝑘=0
Terbukti benar untuk 𝑛 = 𝑚 + 1, 𝑚 ∈ ℕ. Bukti lengkap. ∎
27
Lampiran 7 Program MATLAB dari algoritme untuk metode Kaczmarz Fungsi kaczmarz function x = kaczmarz(A,b,x0) %--------------------------------------------------------------------------------% Fungsi kaczmarz menghasilkan penyelesaian hampiran atas sistem persamaan linear % (SPL) Ax=b : % a11*x1 + a12*x2 + ... + a1N*xN = b1 % a21*x1 + a22*x2 + ... + a2N*xN = b2 % : % aM1*x1 + aM2*x2 + ... + aMN*xN = b3 % menggunakan metode Kaczmarz. %--------------------------------------------------------------------------------% Input: % A = matriks berorde M x N (matriks koefisien dari SPL) % a11 a12 ... a1N % a21 a22 ... a2N % : % aM1 aM2 ... aMN % b = vektor kolom berorde M x 1 (vektor konstanta dari SPL) % b1 % b2 % : % bM % x0 = vektor kolom berorde N x 1 (penyelesaian hampiran awal) % % Input tambahan: % iter = banyaknya iterasi % harus berupa bilangan asli % default: 100 % atau % tol = batas toleransi % harus berupa bilangan real positif % default: 10^(-6) % % Output: % x = vektor kolom berorde Nx1 (penyelesaian hampiran) % x1 % x2 % : % xN %--------------------------------------------------------------------------------% Cek banyaknya input if (nargin < 3) error('Input terlalu sedikit.'); elseif (nargin > 3) error('Input terlalu banyak.'); end [M,N] = size(A); % Orde A, b, dan x0 harus cocok if (size(b,1) ~= M) error('Orde matriks A dan vektor b tidak cocok.'); elseif (size(b,2) ~= 1) error('b harus berbentuk vektor kolom.'); elseif (size(x0,1) ~= N || size(x0,2) ~= 1) error('Orde x0 tidak cocok dengan masalah.') end x = zeros(N,1); % Penentuan kriteria pemberhentian disp('Apakah kriteria pemberhentian akan ditentukan sendiri?'); disp('Pilih : 1 atau 0'); disp('1 : Ya'); disp('0 : Tidak'); jawab = input('Jawab : '); if (isempty(jawab)) error('Tidak ada jawaban.');
28
Lampiran 7 (lanjutan) elseif (jawab ~= 1 && jawab ~= 0) error('Jawaban harus berupa bilangan 1 atau 0.'); elseif (jawab == 0) disp('Kriteria pemberhentian:') disp('Banyaknya iterasi = 100 atau batas toleransi = 10^(-6).'); maksiter = 100; % Nilai default tol = 10^(-6); % Nilai default iter = 1; for k = 1:M x = x0+(b(k)-A(k,:)*x0)*A(k,:)'/(norm(A(k,:))^2); x0 = x; end d = norm(A*x-b); while (iter <= maksiter && d > tol) iter = iter+1; for i = 1:iter for k = 1:M x = x0+(b(k)-A(k,:)*x0)*A(k,:)'/(norm(A(k,:))^2); x0 = x; end end d = norm(A*x-b); end disp(['Pada iterasi ke-', num2str(iter)]); disp(['diperoleh penyelesaian hampiran dengan norm sisaan sebesar ', num2str(d)]); else disp('Pilih kriteria pemberhentian (1, 2, atau 3)'); disp('1 : Banyaknya iterasi'); disp('2 : Batas toleransi'); disp('3 : Keduanya'); itertol = input('Pilihan : '); if (isempty(itertol)) error('Tidak ada kriteria pemberhentian yang dipilih.'); elseif (itertol ~= 1 && itertol ~= 2 && itertol ~=3) error('Input harus berupa bilangan 1, 2 atau 3.'); elseif (itertol == 1) disp('Masukkan banyaknya iterasi!'); disp('Banyaknya iterasi harus berupa bilangan asli.'); iter = input('Banyaknya iterasi = '); if (isempty(iter)) error('Tidak ada input sebagai banyaknya iterasi.'); elseif (iter < 1 || (fix(iter)-iter) ~= 0) error('Banyaknya iterasi harus berupa bilangan asli.'); else for i = 1:iter for k = 1:M x = x0+(b(k)-A(k,:)*x0)*A(k,:)'/(norm(A(k,:))^2); x0 = x; end end d = norm(A*x-b); disp(['Pada iterasi ke-', num2str(iter)]); disp(['diperoleh penyelesaian hampiran dengan norm sisaan sebesar ', num2str(d)]); end elseif (itertol == 2) disp('Masukkan batas toleransi!'); disp('Batas toleransi harus berupa bilangan real positif.'); tol = input('Batas toleransi = '); if (isempty(tol)) error('Tidak ada input sebagai batas toleransi.'); elseif (tol <= 0) error('Batas toleransi harus berupa bilangan real positif.'); else iter = 1;
29
Lampiran 7 (lanjutan) for k = 1:M x = x0+(b(k)-A(k,:)*x0)*A(k,:)'/(norm(A(k,:))^2); x0 = x; end d = norm(A*x-b); while d > tol iter = iter+1; for k = 1:M x = x0+(b(k)-A(k,:)*x0)*A(k,:)'/(norm(A(k,:))^2); x0 = x; end d = norm(A*x-b); end disp(['Pada iterasi ke-', num2str(iter)]); disp(['diperoleh penyelesaian hampiran dengan norm sisaan sebesar ', num2str(d)]); end else disp('Pertama, masukkan banyaknya iterasi!'); disp('Banyaknya iterasi harus berupa bilangan asli.'); iter = input('Banyaknya iterasi = '); if (isempty(iter)) error('Tidak ada input sebagai banyaknya iterasi.'); elseif (iter < 1 || (fix(iter)-iter) ~= 0) error('Banyaknya iterasi harus berupa bilangan asli.'); else maksiter = iter; end disp('Kedua, masukkan batas toleransi!'); disp('Batas toleransi harus berupa bilangan real positif.'); tol = input('Batas toleransi = '); if (isempty(tol)) error('Tidak ada input sebagai batas toleransi.'); elseif (tol <= 0) error('Batas toleransi harus berupa bilangan real positif.'); end iter = 1; for k = 1:M x = x0+(b(k)-A(k,:)*x0)*A(k,:)'/(norm(A(k,:))^2); x0 = x; end d = norm(A*x-b); while (iter <= maksiter && d > tol) iter = iter+1; for i = 1:iter for k = 1:M x = x0+(b(k)-A(k,:)*x0)*A(k,:)'/(norm(A(k,:))^2); x0 = x; end end d = norm(A*x-b); end disp(['Pada iterasi ke-', num2str(iter)]); disp(['diperoleh penyelesaian hampiran dengan norm sisaan sebesar ', num2str(d)]); end end end
30
Lampiran 8 Program MATLAB untuk membangkitkan matriks koefisien dan vektor konstanta dari SPL ke-1 Fungsi blur function [A,b] = blur(N,band,sigma) %--------------------------------------------------------------------------------% Masalah uji BLUR: Memburamkan gambar digital. % % Matriks A adalah matriks simetrik berukuran (N^2)x(N^2), matriks Toeplitz blok % ganda yang memodelkan pemburaman suatu gambar berukuran NxN dengan suatu fungsi % pemancaran titik Gauss. Ini disimpan dalam format matriks sparse. % % Pada setiap blok Toeplitz, hanya elemen-elemen matriks dalam suatu jarak band-1 % dari diagonal yang merupakan taknol (band adalah setengah bandwidth). Jika band % tidak dirinci, digunakan band = 3. % % Parameter sigma mengontrol lebar dari fungsi pemancaran titik Gauss dan % akibatnya jumlah pemulusan. Jika sigma tidak dirinci, digunakan sigma = 0.7. % % Vektor x adalah suatu versi kolom yang ditumpuk dari suatu gambar uji sederhana, % sedangkan vektor b adalah suatu versi kolom yang ditumpuk dari gambar yang % diburamkan, yaitu b = A*x. %--------------------------------------------------------------------------------% Inisialisasi if (nargin < 2) band = 3; end band = min(band,N); if (nargin < 3) sigma = 0.7; end % z A A A
Mengonstruksi matriks koefisien A = [exp(-([0:band-1].^2)/(2*sigma^2)),zeros(1,N-band)]; = toeplitz(z); = sparse(A); = (1/(2*pi*sigma^2))*kron(A,A);
% Mengonstruksi vektor konstanta b x = zeros(N,N); N2 = round(N/2); N3 = round(N/3); N6 = round(N/6); N12 = round(N/12); T = zeros(N6,N3); for i = 1:N6 for j = 1:N3 if ((i/N6)^2 + (j/N3)^2 < 1) T(i,j) = 1; end end end T = [fliplr(T),T]; T = [flipud(T);T]; x(2+[1:2*N6],N3-1+[1:2*N3]) = T; T = zeros(N6,N3); for i = 1:N6 for j = 1:N3 if ((i/N6)^2 + (j/N3)^2 < 0.6) T(i,j) = 1; end end end T = [fliplr(T),T]; T = [flipud(T);T]; x(N6+[1:2*N6],N3-1+[1:2*N3]) = x(N6+[1:2*N6],N3-1+[1:2*N3]) + 2*T;
31
Lampiran 8 (lanjutan) f = find(x==3); x(f) = 2*ones(size(f)); T = triu(ones(N3,N3)); [mT,nT] = size(T); x(N3+N12+[1:nT],1+[1:mT]) = 3*T; T = zeros(2*N6+1,2*N6+1); [mT,nT] = size(T); T(N6+1,1:nT) = ones(1,nT); T(1:mT,N6+1) = ones(mT,1); x(N2+N12+[1:mT],N2+[1:nT]) = 4*T; x = reshape(x(1:N,1:N),N^2,1); b = A*x; end
32
Lampiran 9 Program MATLAB untuk membangkitkan matriks koefisien dan vektor konstanta dari SPL ke-2 Fungsi phillips function [A,b] = phillips(n) %--------------------------------------------------------------------------------% Diskretisasi persamaan Fredholm jenis pertama ditemukan oleh D. L. Phillips. % Didefinisikan fungsi % phi(x) = | 1 + cos(x*pi/3) , |x| < 3 . % | 0 , |x| >= 3 % maka kernel K, penyelesaian f, dan ruas kanan g diberikan oleh: % K(s,t) = phi(s-t) , % f(t) = phi(t) , % g(s) = (6-|s|)*(1+.5*cos(s*pi/3)) + 9/(2*pi)*sin(|s|*pi/3) . % Interval integrasi adalah [-6,6]. % % Didiskretisasi dengan metode Galerkin. % % Input : n harus kelipatan 4. %--------------------------------------------------------------------------------% Cek input if (rem(n,4)~=0) error('n harus kelipatan 4'); end % Mengonstruksi matriks koefisien A h = 12/n; n4 = n/4; r1 = zeros(1,n); c = cos((-1:n4)*4*pi/n); r1(1:n4) = h + 9/(h*pi^2)*(2*c(2:n4+1) - c(1:n4) - c(3:n4+2)); r1(n4+1) = h/2 + 9/(h*pi^2)*(cos(4*pi/n)-1); A = toeplitz(r1); % Mengonstruksi vektor konstanta b b = zeros(n,1); c = pi/3; for i = n/2+1:n t1 = -6 + i*h; t2 = t1 - h; b(i) = t1*(6-abs(t1)/2) ... + ((3-abs(t1)/2)*sin(c*t1) - 2/c*(cos(c*t1) - 1))/c ... - t2*(6-abs(t2)/2) ... - ((3-abs(t2)/2)*sin(c*t2) - 2/c*(cos(c*t2) - 1))/c; b(n-i+1) = b(i); end b = b/sqrt(h); end
33
Lampiran 10 Program MATLAB untuk membangkitkan matriks koefisien dan vektor konstanta dari SPL ke-3 Fungsi tomo function [A,b,x] = tomo(N,f) %--------------------------------------------------------------------------------% TOMO: Membuat suatu masalah uji tomografi 2D % % Fungsi ini membuat suatu masalah uji tomografi dua dimensi sederhana. Suatu % domain 2D [0,N]x[0,N] dibagi ke dalam N^2 sel dengan ukuran satuan dan sejumlah % round(f*N^2) sinar pada arah acak yang menembus domain ini. Nilai default % adalah f=1. % % Setiap sel bersesuaian dengan suatu nilai dalam vektor x dan setiap sinar % bersesuaian dengan elemen dalam vektor konstanta b, yaitu % sum_{cells in ray} x_{cell j} * length_{cell j} % dengan length_{cell j} adalah panjang sinar pada sel ke-j. % % Matriks A adalah sparse, dan setiap baris (bersesuaian dengan suatu sinar) % memegang nilai length_{cell j} pada posisi ke-j. Oleh karena itu: % b = A*x . % Apabila suatu penyelesaian x_reg telah dihitung, ini dapat divisualisasikan % dengan mengartikan imagesc(reshape(x_reg,N,N)). % % Penyelesaian eksak, reshape(x,N,N), identik dengan gambar eksak pada fungsi % blur. %--------------------------------------------------------------------------------% Nilai default if nargin==1 f = 1; end % Mengonstruksi matriks koefisien A A = spalloc(N,N,2*N); x = (0:N)'; y = x; for i = 1:round(f*N^2) x0 = N*rand; x1 = N*rand; y0 = N*rand; y1 = N*rand; a = (y1-y0)/(x1-x0); b = y0 - a*x0; yp = a*x + b; xp = (y - b)/a; xp = [x;xp]; yp = [yp;y]; [xp,I] = sort(xp); yp = yp(I); I = find(xp >= 0 & xp <= N & yp >= 0 & yp <= N); xp = xp(I); yp = yp(I); I = find(diff(xp)==0); xp(I) = []; yp(I) = []; d = sqrt( diff(xp).^2 + diff(yp).^2 ); xm = 0.5*(xp(1:end-1)+xp(2:end)); ym = 0.5*(yp(1:end-1)+yp(2:end)); j = ( floor(xm) )*N + floor(ym) + 1; A(i,j) = d'; end
34
Lampiran 10 (lanjutan) x = zeros(N,N); N2 = round(N/2); N3= round(N/3); N6 = round(N/6); N12 = round(N/12); T = zeros(N6,N3); for i = 1:N6 for j = 1:N3 if ((i/N6)^2 + (j/N3)^2 < 1) T(i,j) = 1; end end end T = [fliplr(T),T]; T = [flipud(T);T]; x(2+(1:2*N6),N3-1+(1:2*N3)) = T; T = zeros(N6,N3); for i = 1:N6 for j = 1:N3 if ((i/N6)^2 + (j/N3)^2 < 0.6) T(i,j) = 1; end end end T = [fliplr(T),T]; T = [flipud(T);T]; x(N6+(1:2*N6),N3-1+(1:2*N3)) = x(N6+(1:2*N6),N3-1+(1:2*N3)) + 2*T; f = find(x==3); x(f) = 2*ones(size(f)); T = triu(ones(N3,N3)); [mT,nT] = size(T); x(N3+N12+(1:nT),1+(1:mT)) = 3*T; T = zeros(2*N6+1,2*N6+1); [mT,nT] = size(T); T(N6+1,1:nT) = ones(1,nT); T(1:mT,N6+1) = ones(mT,1); x(N2+N12+(1:mT),N2+(1:nT)) = 4*T; x = reshape(x(1:N,1:N),N^2,1); % Mengonstruksi vektor konstanta b b = A*x; end
35
Lampiran 11 Norm sisaan dari penyelesaian hampiran atas ketiga SPL yang digunakan yang diperoleh pada iterasi ke-0 sampai iterasi ke-30 Iterasi ke0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
SPL ke-1 15.4392 10.6021 7.2495 4.7998 3.0725 2.0196 1.3670 0.9012 0.5614 0.3748 0.2860 0.2406 0.2104 0.1871 0.1685 0.1538 0.1420 0.1322 0.1237 0.1163 0.1097 0.1039 0.0986 0.0938 0.0894 0.0853 0.0816 0.0782 0.0750 0.0720 0.0693
Norm sisaan dari penyelesaian hampiran SPL ke-2 SPL ke-3 15.2906 207.3483 14.9492 31.1586 13.0457 14.0032 3.9800 8.8950 1.9818 6.6116 0.9924 5.5046 0.5148 4.8414 0.2770 4.3556 0.1597 3.9637 0.1062 3.6390 0.0859 3.3698 0.0798 3.1479 0.0782 2.9659 0.0775 2.8169 0.0766 2.6945 0.0754 2.5934 0.0737 2.5089 0.0717 2.4373 0.0695 2.3756 0.0671 2.3217 0.0646 2.2739 0.0621 2.2308 0.0595 2.1917 0.0570 2.1559 0.0545 2.1228 0.0520 2.0923 0.0497 2.0639 0.0474 2.0375 0.0452 2.0129 0.0430 1.9901 0.0410 1.9688
