Kuliah #7 Pemodelan TK Lanjut S 2 (Tambahan) CONTOH RINGKAS: Solusi SPANL (Sistem Persamaan Aljabar Non Linear)

Kuliah #7 – Pemodelan TK Lanjut S2 (Tambahan)  CONTOH RINGKAS: Solusi SPANL (Sistem Persamaan Aljabar Non Linear)  Prof. Dr. Ir. Setijo Bismo, DEA.  ‐  Departemen Teknik Kimia FTUI, Oktober 2015 

A. Sistem Persamaan Aljabar Non Linear dengan 2 Variabel: x dan y Diberikan SPANL dengan 2 buah variabel berikut ini:

f1 ( x, y ) ≡ x 2 − y − 1 = 0, f 2 ( x, y ) ≡ y 2 − x = 0. Carilah solusi SPANL di atas menggunakan Metode Newton-Raphson dengan harga-harga awal sebagai berikut: a. x0 = 0, y0 = − 1 b.

( (1, 1)

)

A.1. Solusi dengan menggunakan MS-Excel: ¨ Langkah #1: susunlah fungsi SPANL tersebut dalam bentuk “vektorial”, sebagai berikut:

⎡ x 2 − y − 1⎤ Fˆ ( x, y ) = ⎢ ⎥ 2 y − x ⎥⎦ ⎣⎢ ¨ Langkah #2: susunlah matriks Jacobi dari fungsi SPANL di atas sebagai berikut:

⎡ ∂f1 ⎢ ∂x  ⎢ J ( x, y ) = ⎢ ∂f 2 ⎢ ∂x ⎣

∂f1 ⎤ ∂y ⎥ ⎡ 2 x −1 ⎤ ⎥ = ⎢ ⎥ ∂f 2 ⎥ ⎣ −1 2 y ⎦ ∂y ⎥⎦

¨ Langkah #3: jika kita susun-ulang dan gunakan formula rekursif Newton-Raphson menjadi:

⎡ x⎤ ⎡ x⎤ ⎡Δx⎤ = + ⎢ y⎥ ⎢ y⎥ ⎢ Δ y ⎥ ; dengan n adalah “iterasi yang ke..” ⎣ ⎦ n +1 ⎣ ⎦ n ⎣ ⎦n atau menjadi

Xˆ n +1 = Xˆ n + Δ Xˆ n dengan

 −1 Δ Xˆ n = − J ( xn , yn ) ⋅ Fˆ ( xn , yn ) Jika kita gunakan harga awal

(x

0

= 0, y0 = − 1) , maka tahap demi tahap menjalankan

iterasinya, adalah sebagai berikut: #Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [1] 

[Iterasi yang ke-0]: ƒ Ini adalah “awal dari iterasi” atau n = 0 , sehingga dapat dihitung ƒ Nilai dari matriks

 J ( x0 , y0 ) pada

⎡ E28 = 2*D28 F28 = -1 ⎤ = ⎡ 0 −1⎤ , ⎢ E29 = -1 F29 =2*D29 ⎥⎦ ⎢⎣ −1 −2 ⎥⎦ ⎣

kemudian ƒ Nilai dari vektor

Fˆ ( x0 , y0 ) pada

ƒ Sehingga vektor

Δ Xˆ 0 dapat dihitung pada sel-sel H28 dan H29 atau

⎡ G28 = D28^2-D29-1⎤ = ⎡0, 000000E + 00 ⎤ ⎢1, 000000E + 00 ⎥ ⎢ G29 = D29^2-D28 ⎥ ⎣ ⎦ ⎣ ⎦

⎡ H28 ⎤ ⎢ H29 ⎥ ⎣ ⎦

dengan “formula MS-Excel {=MMULT(MINVERSE(E28:F29);G28:G29)}”, yang hasilnya adalah ⎡−1,000000E + 00⎤ ⎢ 1,000000E + 00 ⎥ ⎣ ⎦

ƒ Hasilnya, untuk harga awal ( x0 = 0, y0 = − 1) , dapat dilihat pada tabel berikut ini:

ƒ Kemudian, untuk pengisian sel-sel berikutnya, dilakukan pada “iterasi yang ke-1” seperti di bawah ini. [Iterasi yang ke-1]: ƒ Pada tahap ini, harga n adalah 1 (satu), dapat dihitung Xˆ 1 = Xˆ 0 + Δ Xˆ 0 pada sel-sel ⎡ D30 ⎤ atau lebih lengkapnya adalah ⎡ x1 = D30 = D28-H28⎤ dan hasilnya ⎢ y = D31 = D29-H29 ⎥ ⎣ 1 ⎦

⎢ D31⎥ ⎣ ⎦ adalah ⎡ 1,0000 ⎤ ⎢ −1,0000 ⎥ ⎣ ⎦

x ƒ Selanjutnya, begitu harga-harga ⎡⎢ 1 ⎤⎥ telah diketahui, dapat dihitung nilai dari ⎣ y1 ⎦

unsur-unsur matriks Jacobi

 J ( x1 , y1 ) pada

⎡ E30 = 2*D30 ⎢ E31 = -1 ⎣

F30 = -1 ⎤ = F31 = 2*D31 ⎥⎦

⎡ 2,0000 −1,0000 ⎤ , kemudian dihitung juga... ⎢ −1,0000 −2,0000 ⎥ ⎣ ⎦

ƒ Nilai dari vektor

Fˆ ( x1 , y1 ) pada

⎡ G30 = D30^2-D31-1⎤ = ⎡1,000000E + 00 ⎤ ⎢ G31 = D31^2-D30 ⎥ ⎢0,000000E + 00 ⎥ ⎦ ⎣ ⎦ ⎣

⎡ H30 ⎤ ⎢ H31⎥ ⎣ ⎦ dengan “formula MS-Excel {=MMULT(MINVERSE(E30:F31);G30:G31)}”,

ƒ Sehingga vektor

Δ Xˆ 1 dapat dihitung pada sel-sel H30 dan H31 atau

yang hasilnya adalah

⎡ 4,000000E − 01 ⎤ ⎢ −2,000000E − 01⎥ ⎣ ⎦

#Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [2] 

ƒ Hasilnya, secara lengkap, dapat dilihat pada tabel berikut ini:

ƒ Kemudian, untuk selanjutnya, sebenarnya pekerjaan kita sudah lebih mudah. Karena kita akan dapat “memanfaatkan” protokol dan atau prosedur-prosedur rutin dari MS-EXCEL yang akan dapat memudahkan pekerjaan kita, seperti di bawah ini. [Iterasi yang ke-2]: ƒ Pada tahap ini, yaitu pada harga n adalah 2 (dua), lakukanlah pembiruan sel-sel (cell block):“ D30 sampai H31 ” seperti ilustrasi di bawah ini:

KLIK tetikus “kiri” dan tahan DI SINI

ƒ Perhatikan tanda panah merah pada ilustrasi di atas: KLIK tombol tetikus (mouse) kiri pada titik tersebut, tetap tahan terus sambil diseret ke bawah..., sampai sel H33 ! Hasilnya adalah sebagai berikut:

ƒ Dapat dilihat, bahwa semua “sel-sel sudah terisi” semuanya (yang diblok BIRU seperti di atas). Perlu diketahui, pengisian kolom n (iterasi) adalah manual. ƒ Dalam aksi seperti di atas, penting untuk diperhatikan adalah isi dari “pasangan sel-sel” G32 & G33 (= harga vektor fungsi SPANL) dan atau H32 & H33 (harga selisih x lama dan baru, atau Δ Xˆ ): jika harga-harga pasangan sel tersebut lebih kecil dari pasangan “homolog” di atasnya, maka berarti prosedur yang kita lakukan sudah benar! Artinya, prosedur komputasi yang kita lakukan sudah menuju jalur konvergen! #Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [3] 

[Iterasi yang ke-3]...dan seterusnya... ƒ Dalam hal ini, jika prosedur komputasi yang kita lakukan sudah dalam jalur konvergen, maka langkah-langkah pembiruan sel-sel (cell block) dapat diteruskan kembali, yaitu dari:“ D32 sampai H33 ” kemudian “seretlah KLIK tetikus” kita sampai didapatkan pasangan-pasangan sel G34 & G35 dan atau H34 & H35 pada harga yang terkecil (sesuai kriteria penghentian iterasi), seperti ilustrasi di bawah ini:

ƒ Atau hasil akhirnya adalah sebagai berikut:

ƒ Dari ilustrasi di atas, dapat kita lihat, bahwa: sesungguhnya “penghentian iterasi” sudah dapat dikatakan “memenuhi syarat” pada saat sampai pada iterasi yang ke4 ( n = 4), yaitu pada saat:

⎡1,311636E-06 ⎤ dan atau Fˆ ( x4 , y4 ) = ⎢ ⎣1,820782E-05 ⎥⎦

⎡ -6,468290E-06 ⎤ Δ Xˆ 4 = ⎢ ⎣ -8,101816E-06 ⎥⎦ ƒ Secara konkret, dapat dikatakan: bila kita hentikan iterasi (perhitungan) pada:

→ n = 4: bila kriteria atau ε =

1,0 × 10

−4

#Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [4] 

→ n = 5: bila kriteria atau ε = → n = 6: bila kriteria atau ε =

1,0 × 10

−10

mendekati eps - machine

ƒ Artinya, kriteria berpengaruh pada saat mana kita harus menghentikan komputasi seperti di atas.

Untuk harga awal ( x0 = 1, y0 = 1) , hasil akhirnya dapat dilihat pada tabel di bawah ini:

A.2. Solusi dengan menggunakan Bahasa Pemrogram: ¨ Misalkan, kita hendak menggunakan Program EZY Pascal, maka “Listing” (coding) programnya dapat ditulis sebagai berikut: Program SPANL2Pers; Const NEQ = 4; Type Real = Extended; SubsEl = 1..NEQ; IVektor = Array[SubsEl] of Integer; RVektor = Array[SubsEl] of Real; RMatriks = Array[SubsEl] of RVektor;

{

Procedure PivotGauss(A : RMatriks; Var X : RVektor; B : RVektor; N : SubsEl); ----------------------------------------------------------PROGRAM : SOLUSI PERSAMAAN LINIER DENGAN STRATEGI PIVOTING DALAM ELEMINASI GAUSSIAN (L.Atkinson, p100-101) VERSI : Presisi ganda A [I] : Matrik-segi dari sistem linier B [I] : Vector parameter/harga dari persamaan linier X [O] : Vector variabel/jawab yang diinginkan N [I] : Jumlah persamaan ======> A x X = B ----------------------------------------------------------- }

#Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [5] 

Var P : IVektor; I,J,K,PI,PIVI : Integer; ABSA,FMUL,GMAX,PIVOT : Real; Begin If (N > NEQ) then Begin Write('(PivotGauss : Problem dimensi (MAKS. ',NEQ); Halt End; For I:=1 to N do P[I] := I; For J:=1 to N-1 do Begin { --- Mencari Pivot sepanjang baris J --- } K := J; GMAX := ABS(A[P[K],J]); For I:=J+1 to N do Begin ABSA := ABS(A[P[I],J]); If (ABSA > GMAX) then Begin GMAX := ABSA; K := I End End; PIVI := P[K]; If (K <> J) then Begin { --- Menukar P(J) dengan P(K) --- } I := P[J]; P[J] := P[K]; P[K] := I End; { --- Akhir pencarian Pivot --- } PIVOT := A[PIVI,J]; For I:=J+1 to N do Begin PI := P[I]; FMUL := A[PI,J]/PIVOT; A[PI,J] := 0.0; If (ABS(FMUL) > 0.0) then Begin { - Mengurangi baris A(PI) dengan A(PIVI) - } For K:=J+1 to N do A[PI,K] := A[PI,K] - FMUL*A[PIVI,K]; B[PI] := B[PI] - FMUL*B[PIVI] End; { --- Akhir dari pengurangan --- } End; End; PI := P[N]; X[N] := B[PI]/A[PI,N]; For I:=N-1 downto 1 do Begin PI := P[I]; { -- Calculate the DOT PRODUCT -- } For K:=I+1 to N do B[PI] := B[PI] - A[PI,K]*X[K]; #Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [6] 

X[I] := B[PI]/A[PI,I] End; End; Procedure FSPANL(x : RVektor; var f : RVektor; n : SubsEl); {Fungsi vektorial dari SPANL} Var i : SubsEl; fx : RVektor; Begin fx[1] := x[1]*x[1] - x[2] - 1; fx[2] := -x[1] + x[2]*x[2]; For i := 1 to n do f[i] := -fx[i]; End; Procedure JACSPANL(x : RVektor; var A : RMatriks; n : SubsEl); {Matriks JACOBI dari Fungsi SPANL} Var i,j,k : SubsEl; fk,fk1 : RVektor; dx : Real; Begin A[1,1] := 2*x[1]; A[1,2] := -1.0; A[2,1] := -1.0; A[2,2] := 2*x[2]; End; Function NormVect(x : RVektor; n : SubsEl) : Real; Var Sum : Real; i : SubsEl; Begin Sum := 0.0; For i := 1 to n do Sum := Sum + Sqr(x[i]); NormVect := Sum End; Var ITER,ITMAX : Integer; I,N : SubsEl; Aij : RMatriks; FB,DX,X : RVektor; Solved,Converged,ItMaxReached : Boolean; FTOL,TOL,NV : Real; Begin Write('Jumlah persamaan/variabel: '); Readln(N); Writeln('Harga Awal untuk X (X0):'); For I := 1 to N do Begin Write('X0[',I,']= '); Readln(X[I]); End; Write('ITMAX: '); Readln(ITMAX); Write(' FTOL: '); Readln(FTOL); Write(' TOL: '); Readln(TOL); {Menghitung harga norma fungsi SPANL} FSPANL(X,FB,N); #Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [7] 

NV := NormVect(FB,N); ITER := 0; If (NV > FTOL) then {Memeriksa harga norma fungsi SPANL} Repeat Inc(ITER); JACSPANL(X,Aij,N); PivotGauss(Aij,DX,FB,N); For I := 1 to N do X[I] := X[I] + DX[I]; FSPANL(X,FB,N); Solved := NormVect(FB,N) < FTOL; If not Solved then Begin FSPANL(X,FB,N); End; Converged := NormVect(DX,N) < TOL; ItMaxReached := ITER >= ITMAX; Until Solved or Converged or ItMaxReached; If not ItMaxReached then Begin Writeln('Hasil Komputasi dari Vektor Jawab [X] adalah:'); For I := 1 to N do writeln('X[',I,'] = ',X[I]); Writeln('Jumlah Iterasi: ',ITER); End else Writeln('Jumlah Iterasi terlampaui (>=',ITMAX,')'); Readln End.

¨ Hasilnya, setelah eksekusi program di atas maka hasilnya adalah sebagai berikut:

¨ Perhatikan: bahwa program EZY Pascal seperti di atas sangat lah panjang sehingga tidak praktis untuk perhitungan-perhitungan (komputasi) yang bersifat segera dan ringkas.

#Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [8] 

B. Sistem Persamaan Aljabar Non Linear dengan 2 Variabel:

x , y , dan z

Di bawah ini diberikan sebuah SPANL dengan 3 buah variabel berikut ini:

f1 ( x, y, z ) ≡ 3 x y + 2 yz f 2 ( x, y, z ) ≡ e x + yz f3 ( x, y, z ) ≡ xz + y cos( z ) Sekarng, agar saudara dapat berlatih sendiri menyelesaikan SPANL 3 persamaan seperti di atas dengan menggunakan Metode Newton-Raphson, maka carilah solusi SPANL di atas dalam spreadsheet MS-EXCEL dengan harga-harga awal ( x0 = 0, y0 = 1, z0 = 2 ) .

#Solusi dengan menggunakan MS-Excel: ¨ Langkah #1: susunlah fungsi SPANL tersebut dalam bentuk “vektorial”, sebagai berikut:

⎡3 x y + 2 yz ⎤ ⎢ ⎥ Fˆ ( x, y ) = ⎢ e x + yz ⎥ ⎢ xz + y cos( z ) ⎥ ⎣ ⎦ ¨ Langkah #2: susunlah matriks Jacobi dari fungsi SPANL di atas sebagai berikut:

⎡ ∂f1 ⎢ ⎢ ∂x ⎢ ∂f 2 J ( x, y , z ) = ⎢ ⎢ ∂x ⎢ ∂f ⎢ 3 ⎣⎢ ∂y

∂f1 ∂y ∂f 2 ∂y ∂f 3 ∂y

∂f1 ⎤ ⎥ ∂z ⎥ ⎡3 y ∂f 2 ⎥ ⎢ x ⎥ = ⎢e ∂z ⎥ ⎢⎣ z ∂f3 ⎥ ⎥ ∂z ⎦⎥

(3 x

+ 2 z) z

cos( z )

⎤ ⎥ ⎥ x − y sin( z ) ⎥⎦ 2y y

¨ Langkah #3: jika kita susun-ulang dan gunakan formula rekursif Newton-Raphson menjadi:

Jika kita gunakan harga awal ( x0 = 0, y0 = 1, z0 = 2 ) , maka hasil spreadsheet dalam MS-EXCEL yang dimaksudkan adalah sebagai berikut: #Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015  [9] 

Tugas buat saudara:

Coba kerjakan sendiri sampai dapat hasilnya seperti pada tabel di atas!

----- Selamat belajar dan sampai jumpa lagi di lain waktu -----

#Property of  Prof. Dr. Ir. Setijo Bismo, DEA.                                Departemen Teknik Kimia FTUI, Oktober 2015    [10]