KOMPUTASI PROSES TEKNIK KIMIA MENGGUNAKAN MATLAB

KOMPUTASI PROSES TEKNIK KIMIA MENGGUNAKAN MATLAB Teguh Kurniawan, ST.

5

F(X,Y )

0

-5

-10 2

3 2 0

1 0 -1

-2 Y

-2 -3

X

Jurusan Teknik Kimia Fakultas Teknik Universitas Sultan Ageng Tirtayasa 2006

Kata Pengantar Segala rasa syukur saya panjatkan pada Allah SWT yang telah memberikan segala daya sehingga terwujudlah buku ajar dihadapan para pembaca. Momen bulan Ramadhan 1427 H telah memberikan kekuatan tersendiri bagi saya untuk memberikan karya yang sebaik-baiknya. Meskipun demikian, saya sadar tiada karya manusia yang sempurna. Oleh karena itu, diawal halaman pengantar saya ungkapkan harapan yang sangat besar kepada para pembaca seluruhnya untuk memberikan saran demi perbaikan buku ajar. Perkembangan perangkat keras (hardware-microprocessor) yang amat cepat, membuat teknologi perangkat lunak (software) menjadi cepat pula berkembang. Berbagai macam perangkat lunak untuk mendukung penyelesaian persoalan keteknikan (engineering) bermunculan. Di mulai dengan Basic, Pascal, FORTRAN sebagai program pendahulu. Bahkan FORTRAN sempat menjadi primadona “senjata ampuh” bagi orang teknik. Belakangan muncul POLYMATH, MATHCAD dan MATLAB sebagai generasi baru perangkat lunak bidang teknik dan MIPA. Perangkat lunak yang terakhir disebut, saat ini sudah sangat banyak penggunanya baik dari kalangan akademisi maupun praktisi karena berbagai keunggulannya dalam menyelesaikan persoalan numerik dan simulasi proses. Salah satu bidang keteknikan yang sarat dengan perhitungan numerik, pemodelan serta simulasi fenomena proses adalah bidang teknik kimia. Penggunaan paling banyak di jurusan teknik kimia adalah pada bagian teknik reaktor, pengendalian proses, kinetika katalisis, termodinamika, operasi perpindahan kalor, peristiwa perpindahan dan tentunya komputasi kroses teknik kimia. Seiring dengan makin luasnya penggunaan MATLAB di kalangan industri, maka sebagai langkah untuk mempersiapkan para mahasiswa agar siap terjun di dunia kerja sudah saatnya bagi setiap kampus yang berorientasi pada kebutuhan pasar untuk mengajarkan MATLAB. Materi buku ajar disusun dengan pendekatan penyelesaian kasus-kasus persoalan numerik teknik kimia dengan menggunakan subrutin MATLAB versi 7.0 seperti roots, fzero, fsolve, quad, ode23, fminsearch dll. Namun demikian dalam buku ajar ini juga tetap dimasukkan berbagai materi asas numerik untuk mencegah terjadinya pengetahuan semu (black box) seperti pemrograman Gauss, Newton-Raphson, bisection dalam bahasa MATLAB. Harapan penulis dengan hadirnya buku ajar ini dapat membantu memperkenalkan MATLAB sekaligus menjadi acuan praktis untuk memahami MATLAB sebagai alat untuk menyelesaikan persoalan-persoalan dalam dunia teknik, khususnya teknik kimia. Akhir kata saya ucapkan terima kasih kepada Universitas Sultan Ageng Tirtayasa (UNTIRTA) atas segala dukungannya. Cilegon, 30 Oktober 2006

Teguh Kurniawan, ST Staf Pengajar Komputasi Proses Teknik Kimia Jurusan Teknik Kimia-Fakultas Teknik Universitas Sultan Ageng Tirtayasa

i

DAFTAR ISI

Halaman Kata Pengantar..........................................................................................................i Daftar Isi..................................................................................................................ii Daftar Tabel............................................................................................................iii Daftar Gambar........................................................................................................iv Bab 1 Pengenalan MATLAB...................................................................................1 Bab 2 Sistem Persamaan Linier.............................................................................27 Bab 3 Persamaan Tak Linier..................................................................................40 Bab 4 Optimisasi....................................................................................................62 Bab 5 Regresi Linier dan Non Linier.....................................................................65 Bab 6 Integrasi.......................................................................................................78 Bab 7 Persamaan Diferensial Biasa.......................................................................83 Bab 8 Persamaan Diferensial Parsial...................................................................109 Daftar Pustaka.........................................................................................................v

ii

Bab 1 Pengenalan MATLAB &Pengantar Pemrograman

Bab 1 Pengenalan Matlab & Pengantar Pemrograman 1.1 Perangkat Lunak MATLAB MATLAB merupakan perangkat lunak produk dari The MathWorks,Inc yang memadukan kemampuan perhitungan, pencitraan, dan permograman dalam satu paket. MATLAB merupakan bahasa komputasi teknik yang lebih mudah dan lebih canggih dalam penggunaannya dibandingkan dengan bahasa teknik pendahulunya seperti FORTRAN, BASIC, PASCAL. Sebetulnya MATLAB tidaklah berbeda dengan kalkulator scientific yang sehari-hari kita (orang teknik) kenal. Bedanya MATLAB adalah kalkulator super canggih, karena MATLAB memiliki keunggulan sbb: 1. Menghitung sampai dengan ketelitian 16 angka dibelakang koma, sehingga perhitungan lebih akurat. 2. Menyediakan fasilitas untuk membuat program sesuai dengan kebutuhan kita. 3. Mampu menampilkan data-data dalam grafik 2-D hingga 3-D dengan pewarnaan yang akan memudahkan interpretasi data yang kita miliki. 4. Menyediakan perintah-perintah praktis untuk menyelesaikan berbagai macam

persoalan

matematis

seperti

persamaan

pangkat

tinggi

(polinomial), persamaan linier, persamaan tak linier, optimasi fungsi, persamaan diferensial biasa, persamaan diferensial parsial, fungsi integral, interpolasi data, operasi aljabar, operasi matrik, korelasi data-data dan masih banyak lagi. 5. Memiliki kemudahan dalam mengelola data-data yang sangat banyak dalam bentuk vektor/matrik. 6. Memiliki fasilitas toolbox yang berisi subrutin untuk menyelesaikan persoalan tertentu dan dapat dengan mudah dimodifikasi serta ditambah untuk pengembangan lebih lanjut.

Halaman 1 dari 101


Secara garis besar lingkungan kerja MATLAB terdiri atas beberapa unsur, yaitu: 1. Command window (layar kendali) 2. Workspace (rak data) 3. Command history (layar pengingat) 4. M-file (editor )Æ akan dibahas pada bagian khusus.

Workspace berfungsi sbg tempat menyimpan secara ototmatis segala variabel masukan dan hasil

Command window merupakan jendela utama MATLAB. Tempat untuk mengeksekusi perintah menampilkan masukan dan hasil

Command history adalah tempat menyimpan secara otomatis segala perintah yang telah dituliskan pada command windows.

Gambar 1.1 Lingkungan kerja MATLAB 7.0 Untuk lebih jelas mengenai lingkungan kerja MATLAB perhatikan contoh berikut ini (lihat gambar 1.2). Pada command window ketikkan a = 2 dan b = 4, maka secara otomatis MATLAB akan menyimpan variabel a dengan harga 2 dan variabel b dengan harga 4 pada workspace. Variabel a dan b dapat dipanggil setiap saat dibutuhkan. Misalkan kita ingin menghitung perkalian a dan b, kemudian menyimpannya dengan nama variabel c. Pada command window ketikkan c = a*b, maka MATLAB akan memanggil harga a dan b kemudian melakukan operasi perkalian dan menyimpan hasilnya dengan nama variabel c. Segala sesuatu yang telah diketikkan pada command window disimpan dalam command history dan dapat dipanggil kembali dengan menggunakan key arrow atas dan bawah (↑↓).

Halaman 2 dari 101


Perintah memasukan data variabel a

Menyimpan secara otomatis harga variabel a, b , dan c

Perintah memasukan data variabel b

Perintah menghitung harga variabel c

Menyimpan secara otomatis perintah-perintah yang telah diketikkan di command window

Gambar 1.2 Sistem kerja MATLAB Sekali kita mendefinisikan sebuah variabel, MATLAB akan menyimpan dalam workspace untuk “selamanya”. Untuk menghapus seluruh variabel yang telah dibuat dapat menggunakan perintah clear. Variabel a, b, dan c yang telah tersimpan akan hilang. Jika ingin membersihkan layar command window tanpa menghapus variabel-variabelnya kita dapat menggunakan perintah clc. Beberapa nama variabel yang telah didefinisikan oleh MATLAB sebagai berikut: Tabel 1.1 Variabel terdefinisi dalam MATLAB Var

Keterangan

eps

Bilangan yang jika ditambahkan dengan suatu bilangan lain tidak mengubah besar bilangan lain itu. Epsilon berharga 2.2204e-016

pi

3.1416....

inf

Tak berhingga (Infinity). Simbol matematika = ~

NaN Bilangan tak tentu (Not a Number) contoh 0/0, ~ - ~ i,j

Bilangan imajiner (akar dari -1)

Halaman 3 dari 101


Untuk menamakan sebuah variabel sebaiknya tidak memakai nama-nama yang telah didefinisikan oleh MATLAB.

1.2 Matrik, Vektor dan MATLAB MATLAB adalah singkatan dari matrix laboratory. Oleh karena itu pemahaman terhadap konsep matrik harus memadai agar dapat memanfaatkan MATLAB sebagai bahasa komputasi dengan maksimal. Vektor merupakan matrik yang hanya terdiri atas satu kolom atau satu baris saja. Penulisan matrik di MATLAB Tanda pisah antar elemen matrik Tanda koma (,) atau spasi digunakan untuk memisahkan elemen-elemen satu baris. Tanda titik koma(;) digunakan untuk memisahkan elemen-elemen satu kolom. >> a=[1,2,3] a = 1

2

3

>> b=[1;2;3] b = 1 2 3

>> A=[1 2 3;4 5 6;7 8 9] A = 1

2

3

4

5

6

7

8

9

Matrik transposisi >> A' ans =

1

4

7

2

5

8

3

6

9

Halaman 4 dari 101


Menentukan ukuran matrik >> size(A) ans = 3

3

Menentukan determinan matrik >> det(A) ans = 0

Menentukan invers matrik >> inv(A) Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.541976e-018. ans = 1.0e+016 * -0.4504

0.9007

-0.4504

0.9007

-1.8014

0.9007

-0.4504

0.9007

-0.4504

Perhitungan

invers

matrik

A

menggunakan

MATLAB

ternyata

memunculkan peringatan yang menyatakan bahwa matrik A adalah singular (tak wajar). Hal ini bisa diketahui lebih awal dengan melihat harga determinan A. Apabila determinan A berharga nol dapat dipastikan matrik A adalah matrik singular. Vektor baris adalah matrik yang terdiri atas satu baris saja. >> B=[2:6] B = 2

3

4

5

6

Penulisan seperti di atas akan menghasilkan vektor baris dengan selisih 1 >> C=[2:2:6] C = 2

4

6

Penulisan seperti di atas akan menghasilkan vektor baris dengan selisih 2

Halaman 5 dari 101


Vektor kolom adalah matrik yang terdiri atas satu kolom saja >> V=[2:0.5:4]' V = 2.0000 2.5000 3.0000 3.5000 4.0000

Penulisan seperti di atas akan menghasilkan vektor kolom dengan selisih 0.5

Menentukan ukuran vektor >> length(V) ans = 5

Matrik kerancang Matrik kerancang adalah matrik berdimensi besar yang sebagian besar elemennya adalah nol. Misalnya: ⎡1 ⎢0 ⎢ ⎢0 ⎢ ⎢0 ⎢1 ⎢ ⎢7 ⎢0 ⎣

0 0 0 0 0 7 0

0 8 0 0 0 0 0

0 0 7 0 0 0 0

3 0 9 0 0 0 0

0 0 0 0 0 0 9

6⎤ 0 ⎥⎥ 0⎥ ⎥ 0⎥ 0⎥ ⎥ 0⎥ 0 ⎥⎦

Matrik segitiga atas Matriks segitiga atas (disimbolkan U atau R) Adalah matriks bujur sangkar yang semua elemen di bawah diagonalnya nol ⎡1 ⎢0 ⎢ ⎢0 ⎢ ⎣0

4⎤ 5 6 7 ⎥⎥ = U4 0 8 9⎥ ⎥ 0 0 10 ⎦ 2 3

Halaman 6 dari 101


Matrik segitiga bawah Matrik segitiga bawah (disimbolkan L) Adalah matriks bujur sangkar yang semua elemen di atas diagonalnya nol. ⎡1 ⎢2 ⎢ ⎢3 ⎢ ⎣4

0 5 6 7

0 0⎤ 0 0 ⎥⎥ = L4 8 0⎥ ⎥ 9 10 ⎦

Matrik identitas Matrik identitas adalah matrik yang elemen diagonalnya bernilai 1 dan elemen lainnya bernilai nol. ⎡1 ⎢0 ⎢ ⎢0 ⎢ ⎣0

0 1 0 0

0 0 1 0

0⎤ 0 ⎥⎥ = I4 0⎥ ⎥ 1⎦

Perhatikan cara membuat matrik identitas tanpa harus mengetik elemen per elemen anggota matrik sbb: >> diag(ones(4,1)) ans = 1

0

0

0

0

1

0

0

0

0

1

0

0

0

0

1

1

0

0

0

0

1

0

0

0

0

1

0

0

0

0

1

atau >> eye(4) ans =

Matriks diagonal Matrik diagonal adalah matrik yang elemen selain diagonalnya bernilai nol.

Halaman 7 dari 101


⎡1 ⎢0 ⎢ ⎢0 ⎢ ⎣0

0 2 0 0

0 0 3 0

0⎤ 0 ⎥⎥ = D4 0⎥ ⎥ 4⎦

Perhatikan cara membuat matrik diagonal tanpa harus mengetik elemen per elemen sbb: >> diag(1:4) ans = 1

0

0

0

0

2

0

0

0

0

3

0

0

0

0

4

Matrik Tridiagonal ⎡ −2 1 0 0 0 0 0 0 0 ⎤ ⎢ 1 −2 1 0 0 0 0 0 0 ⎥ ⎢ ⎥ ⎢ 0 1 −2 1 0 0 0 0 0 ⎥ ⎢ ⎥ ⎢ 0 0 1 −2 1 0 0 0 0 ⎥ ⎢ 0 0 0 1 −2 1 0 0 0 ⎥ ⎢ ⎥ ⎢ 0 0 0 0 1 −2 1 0 0 ⎥ ⎢ 0 0 0 0 0 1 −2 1 0 ⎥ ⎢ ⎥ ⎢ 0 0 0 0 0 0 1 −2 1 ⎥ ⎢⎢ 0 0 0 0 0 0 0 1 −2 ⎥⎥ ⎣ ⎦

Perhatikan cara membuat matrik diagonal tanpa harus mengetik elemen per elemen sbb: >> diag(-2*ones(9,1))+diag(ones(8,1),1)+diag(ones(8,1),-1) ans = -2

1

0

0

0

0

0

0

0

1

-2

1

0

0

0

0

0

0

0

1

-2

1

0

0

0

0

0

0

0

1

-2

1

0

0

0

0

0

0

0

1

-2

1

0

0

0

0

0

0

0

1

-2

1

0

0

0

0

0

0

0

1

-2

1

0

0

0

0

0

0

0

1

-2

1

0

0

0

0

0

0

0

1

-2

Halaman 8 dari 101


Aljabar Matrik Operasi aljabar matrik maupun skalar menggunakan simbol yang tidak jauh berbeda. Berikut ini hirarki operasi aljabar dalam MATLAB. Pertama ^ kedua * ketiga / atau \ dan terakhir + dan -. Keterangan: ^

Pangkat

*

Perkalian

/

Pembagian matrik kanan (mis: B/A = B*inv(A))

\

Pembagian matrik kiri (mis: A\B = inv(A)*B)

+

Penambahan

-

Pengurangan

Penjumlahan dan pengurangan Hanya dapat dilakukan jika matrik-matrik yang akan dijumlahkan dan dikurangkan memiliki orde sama. ⎡ 2 3 1 6 ⎤ ⎡ 2 3 1 6 ⎤ ⎡ 4 6 2 12 ⎤ ⎢ 1 4 5 2 ⎥ + ⎢1 4 5 2 ⎥ = ⎢ 2 8 10 4 ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎡ 2 3 1 6 ⎤ ⎡ 2 3 1 6 ⎤ ⎡0 0 0 0⎤ ⎢1 4 5 2⎥ − ⎢1 4 5 2⎥ = ⎢0 0 0 0⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦

>> A =[2 3 1 6;1 4 5 2] A =

2

3

1

6

1

4

5

2

4

6

2

12

2

8

10

4

0

0

0

0

0

0

0

0

>> A+A ans =

>> A-A ans =

Halaman 9 dari 101


Perkalian matrik Syarat Æ jumlah kolom A = jumlah kolom baris B

AB

AB ≠ BA Misal ⎡1 ⎤ B = ⎢⎢ 2 ⎥⎥ ⎢⎣ 3 ⎥⎦

A = [1 2 3]

⎡1 ⎤ AB = [1 2 3] ⎢⎢ 2 ⎥⎥ = 1 + 4 + 9 = 14 ⎣⎢ 3 ⎦⎥

(1 x 3) (3 x 1)=(1 x 1) ⎡1 ⎤ ⎡1 2 3⎤ ⎢ ⎥ BA = ⎢ 2 ⎥ [1 2 3] = ⎢⎢ 2 4 6 ⎥⎥ ⎢⎣ 3 ⎥⎦ ⎢⎣ 3 6 9 ⎥⎦

Operasi perkalian matrik dalam MATLAB dilakukan dengan simbol * >> A=[1,2,3] A = 1

2

3

>> B=[1;2;3] B = 1 2 3 >> A*B ans = 14 >> B*A ans = 1

2

3

2

4

6

3

6

9

Halaman 10 dari 101


Pembagian matrik kanan

xA = c x = cA−1 x = c/ A Misalkan: ⎡1 2 3⎤ x ⎢⎢ 2 5 4 ⎥⎥ = [ 20 15 −8] ⎢⎣ 4 3 1 ⎥⎦ >> A=[1 2 3;2 5 4;4 3 1] A = 1

2

3

2

5

4

4

3

1

>> c=[20 15 -8] c = 20

15

-8

>> x=c/A x = -8.6667

3.0952

5.6190

Pembagian matrik kiri

Ax = c x = A−1c x = A\c Misalkan: ⎡1 2 3⎤ ⎡ 20 ⎤ ⎢ 2 5 4 ⎥ x = ⎢15 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⎣ 4 3 1 ⎥⎦ ⎢⎣ −8⎥⎦

>> A=[1 2 3;2 5 4;4 3 1] A = 1

2

3

2

5

4

4

3

1

Halaman 11 dari 101

Bab 1 Pengenalan MATLAB &Pengantar Pemrograman >> c=[20;15;-8] c = 20 15 -8

>> x=A\c x = -1.0000 -4.7143 10.1429

Beberapa fungsi built in matrik -perintah yang telah terdefinisi dalam MATLAB- akan sangat berguna untuk mempermudah pekerjaan perhitungan matrik, disajikan sebagai berikut.

expm

eksponensial dari sebuah matrik

logm

logaritma dari sebuah matrik

sqrtm

akar kuadrat dari sebuah matrik

Operasi elemen matrik Seringkali dibutuhkan operasi antar elemen-elemen matrik, oleh karena itu MATLAB telah menyediakan perintah untuk melakukan operasi elemen matrik dengan simbol .* (titik diikuti dengan bintang). Perkalian elemen hanya dapat dilakukan untuk orde matrik yang sama. 1. .* perkalian antar elemen matrik. A.*B adalah perkalian antar elemen per elemen matrik A dengan B. A dan B harus memiliki ukuran yang sama kecuali jika salah satunya adalah skalar (bilangan tunggal). 2. ./ Pembagian elemen kanan. A./B adalah matrik dengan elemen-elemen A(i,j)/B(i,j). A dan B harus memiliki ukuran yang sama, kecuali jika salah satunya adalah skalar.

Halaman 12 dari 101


3. .\ Pembagian elemen kiri. A.\B adalah matrik dengan elemen-elemen B(i,j)/A(i,j). A dan B harus memiliki ukuran yang sama, kecuali jika salah satunya adalah skalar. 4. .^ Pangkat elemen. A.^B is adalah matrik dengan elemen-elemen A(i,j) pangkat B(i,j). A dan B harus memiliki ukuran yang sama, kecuali jika salah satunya adalah skalar. Berikut ini masing-masing contoh operasi elemen matrik. >> A=[1 2;3 4] A = 1

2

3

4

>> A.*A ans = 1

4

9

16

>> A./A' ans =

1.0000

0.6667

1.5000

1.0000

>> A.\A' ans = 1.0000

1.5000

0.6667

1.0000

>> A.^A ans = 1

4

27

256

Halaman 13 dari 101


1.3 Membuat Grafik Grafik 2 Dimensi Perintah menggambar grafik 2D

plot(x,y) Misalkan: x

1

2

3

4

5

y

2.7

7.4

20.1

54.6

148.4

>> x=[1,2,3,4,5] x = 1

2

3

4

5

>> y=[2.7,7.4,20.1,54.6,148.4] y = 2.7000

7.4000

20.1000

54.6000

148.4000

150

y

100

50

>> plot(x,y) >> xlabel('x') >> ylabel('y')

0

1

1.5

2

2.5

3 x

3.5

4

4.5

5

Gambar 1.3 Grafik 2 Dimensi Grafik 3 Dimensi Perintah menggambar grafik 3D

surf(x,y,z) Misalkan: x

y

z(x=1)

z(x=2)

z(x=3)

1

1

2

5

10

2

2

5

8

13

3

3

10

13

18

4

17

20

25

Halaman 14 dari 101

Bab 1 Pengenalan MATLAB &Pengantar Pemrograman >> x=[1 2 3] x =

1

2

3

>> y=[1 2 3 4]

y =

1

2

3

4

>> z=[2 5 10;5 8 13;10 13 18;17 20 25]

z = 25

5

10

20

5

8

13

15

10

13

18

17

20

25

z

2

10 5

>> surf(x,y,z) >> xlabel('x')

0 4 3

3

2.5

>> zlabel('z')

2

2

>> ylabel('y') y

1.5 1

1

x

Gambar 1.4 Grafik 3 Dimensi Untuk mempercantik tampilan dan mempermudah penafsiran grafik dengan menambah legenda warna ketikkan perintah berikut ini. >> shading interp >> colorbar

Halaman 15 dari 101


Gambar 1.5 Grafik 3 Dimensi yang diperhalus

Grafik 3 Dimensi Semu Apabila penafsiran grafik 3D seperti tercetak di muka masih dirasakan sulit, MATLAB telah menyediakan perintah untuk membuat grafik 3D menjadi grafik 2D. >> pcolor(x,y,z) >> xlabel('x') >> ylabel('y') >> zlabel('z') >> shading interp >> colorbar

Gambar 1.6 Grafik 3 Dimensi semu

1.4 Algoritma & Pemrograman Algoritma adalah urutan langkah-langkah logis yang dibutuhkan untuk melakukan suatu tugas spesifik. Algoritma dapat dituliskan dalam bentuk kalimat, namun lebih umum dituliskan dalam bentuk diagram alir (flow chart).

Halaman 16 dari 101


Tabel 1.2 Simbol algoritma pemrograman Simbol

Nama

Fungsi

Garis alir

Menyatakan aliran logika

Terminal

Menyatakan awal atau akhir suatu program.

Proses

Menyatakan perhitungan atau manipulasi data.

Masukan/keluaran

Menyatakan masukan atau keluaran data dan informasi

Kondisi/keputusan

Menyatakan sebuah perbandingan, pertanyaan atau keputusan yang menentukan lintasan mana yang akan diikuti.

Konektor

Menyatakan perpindahan halaman.

M-file Sampai saat ini kita masih menjalankan perintah-perintah serta masukan data dengan mengetikkannya secara langsung dalam command window. Tentunya akan sangat merepotkan jika kita dihadapkan pada persoalan yang menuntut pembuatan program yang sangat panjang berpuluh-puluh bahkan beratus-ratus baris perintah. Untuk kemudahan dalam membuat program, MATLAB menyediakan fasilitas m-file atau editor sebagai tempat untuk mengetikkan perintah dan menyimpan program-program yang dibuat. Penulisan program dalam

Halaman 17 dari 101


m-file dapat dilakukan dengan dua cara, yaitu skrip m-file dan fungsi m-file. Berikutnya akan dibahas satu per satu.

Cara membuka M-File: File/New/M-File

Gambar 1.7 Membuka m-file Aturan penamaan m-file Penamaan m file baik untuk skrip maupun fungsi memiliki aturan tertentu yang harus dipatuhi. Berikut ini aturan penamaan m-file pada MATLAB 7. 1. Penamaan harus dimulai dengan huruf latin (a-z atau A-Z) baru kemudian boleh diikuti dengan angka. Huruf kapital dengan huruf kecil tidaklah sama (FILE ≠ file) 2. Tidak boleh ada spasi, titik, koma, titik koma, dan segala macam tanda baca lainnya kecuali underscore ( _ ). 3. Nama sebuah fungsi m-file sebaiknya disamakan dengan nama fungsinya. 4. Sebaiknya tidak menggunakan nama yang telah didefinisikan sebagai fungsi MATLAB tertentu, contoh roots, fzero,zeros dll.

Halaman 18 dari 101


Skrip m-file Skrip adalah file sederhana yang tidak memiliki input argumen dan output argumen. Definisi lain yang mudah diingat, skrip adalah penulisan program MATLAB dalam m-file dengan bentuk bukan fungsi. Sebagai contoh berikut ini adalah penulisan perintah-perintah dalam m-file untuk membuat grafik 3D yang telah dituliskan sebelumnya secara langsung pada command window. File ini disimpan dengan nama coba_m_file.m

Gambar 1.7 Skrip coba_m_file.m Untuk memberikan komentar dalam m-file dapat dilakukan dengan menambahkan % sebelum mengetikkan komentar atau keterangan yang diperlukan seperti terlihat pada coba_m_file di atas. Eksekusi atau menjalankan skrip tersebut dapat dilakukan dengan berbagai cara yang berbeda sbb: Halaman 19 dari 101


1. Tekan tombol F5 pada keyboard, atau 2. Klik debug kemudian run, atau 3. Aktifkan command window. Ketikkan nama file yang akan dieksekusi. >>coba_m_file.

Fungsi m-file Seperti yang telah dijelaskan sebelumnya selain dengan skrip kita dapat juga melakukan pemrograman dalam bentuk lain yaitu fungsi m-file. Penjelasan mengenai cara membuat fungsi m-file dilakukan dengan pendekatan contoh soal kasus 1. Namun sebelum menginjak pada pembahasan cara membuat m-file, saya akan mengajak untuk melihat beberapa fungsi yang telah ada dalam MATLAB sebagai fungsi built in sebagai berikut. Tabel 1.3 Fungsi built in MATLAB

Fungsi

Keterangan

sin(x) sind(x) cos(x) cosd(x) tan(x) tand(x) log(x) log10(x) log2(x) exp(x) sqrt(x)

harga sinus dari x, radian harga sinus dari x, derajat harga kosinus dari x, radian harga kosinus dari x, derajat harga tangen dari x, radian harga tangen dari x, derajat logaritma dengan basis bilangan natural e dari x logaritma dengan basis bilangan 10 dari x logaritma dengan basis bilangan 2 dari x eksponensial dari x akar kuadrat dari x

Kasus 1 [volume tangki penyimpan] Senyawa kimia yang mudah menguap pada temperatur kamar biasa disimpan dalam fasa cair pada tekanan uapnya. Dalam kasus ini n-butana (C4H10) di simpan pada tekanan 2,581 bar dan temperatur 300 K. Penyimpanan skala besar (bulk>50 m3) n-butana seringkali dilakukan dalam tangki yang berbentuk bola (spherical). Sebuah tangki penyimpan n-butana berbentuk bola. Hitunglah volume tangki jika bola memiliki jari-jari 2,3,……9,10 m!.

Halaman 20 dari 101


Jawaban: 4 Vbola = π r 3 3

Algoritma pemrograman Mulai

Masukan harga Jari-jari, r (m)

Hitung harga Volume bola V = 4/3*π*r^3

Harga Volume Bola (m3)

Selesai Penulisan program untuk kasus 1 kita dilakukan dengan dua cara, yaitu dalam bentuk skrip dan fungsi .

Halaman 21 dari 101


Penulisan program dalam bentuk skrip % kasus_1.m clc clear r = 2:10 V =4/3*pi*r.^3 % Membuat grafik V terhadap r plot(r,V) xlabel('jari-jari [m]') ylabel('Volume [m^3 ]')

Eksekusi kasus_1.m dalam command window >>kasus_1 r = 2

3

4

5

V = 1.0e+003 * Columns 1 through 5 0.0335 0.1131 Columns 6 through 9 1.4368 2.1447

6

7

8

0.2681

0.5236

3.0536

4.1888

9

10

0.9048

4500 4000 3500

Volume [m3 ]

3000 2500 2000 1500 1000 500 0

2

3

4

5

6 jari-jari [m]

7

8

9

10

Gambar 1.8 Volume vs jari-jari tangki penyimpan [skrip]

Halaman 22 dari 101


Penulisan program dalam bentuk fungsi %kasus1.m function V = kasus1(r) V = 4/3*pi*r.^3 % Membuat gambar plot(r,V) xlabel('jari-jari [m]') ylabel('volume [m^3 ]')

Eksekusi fungsi kasus1.m di command window >> kasus1(2:10) ans = 1.0e+003 * Columns 1 through 5 0.0335 0.1131 Columns 6 through 9 1.4368 2.1447

0.2681

0.5236

3.0536

4.1888

0.9048

Eksekusi sebuah fungsi dapat pula dilakukan dengan perintah berikut ini.

feval(‘fungsi’,x1,...,xn) x1,....,xn adalah varibel bebas yang akan dievaluasi. >> feval('kasus1',[2:10]) ans = 1.0e+003 * Columns 1 through 5 0.0335

0.1131

0.2681

0.5236

3.0536

4.1888

0.9048

Columns 6 through 9 1.4368

2.1447 4500 4000 3500

Volume [m3 ]

3000 2500 2000 1500 1000 500 0

2

3

4

5

6 jari-jari [m]

7

8

9

10

Gambar 1.8 Volume vs jari-jari tangki penyimpan [fungsi]

Halaman 23 dari 101


Kontrol aliran MATLAB memiliki kontrol aliran yang berguna dalam menentukan berbagai keputusan selanjutnya sebuah program, diantaranya adalah for, if, while, dan switch. Pernyataan relasi yang sering digunakan dalam kontrol aliran adalah sebagai berikut: ==

sama dengan

>

lebih besar dari

<=

kurang dari sama dengan

~

operator logika tidak

>=

lebih dari sama dengan

&

operator logika dan

~=

tidak sama dengan

|

operator logika atau

<

kurang dari

for % ideal.m clear clc R = 8.314;

%J/mol.K

T =[310:10:400];

%K

P =1e5;

%Pa

for i = 1:10 V(i)=R*T(i)/P;

%m3/mol

end V

if %nilai.m x = input('masukkan nilai ujian='); if x >= 80 disp('Nilai A') elseif x >= 65 disp('Nilai B') elseif x >= 55 disp('Nilai C') elseif x >= 45 disp('Nilai D') else disp('Nilai E') end

Halaman 24 dari 101


While %diff.m dif=1; x2=7 while dif > 0.0005 x1=x2-cos(x2)/(1+x2); dif=abs(x2-x1); x2=x1 end

Switch %pilih.m method = 'Bilinear';

switch lower(method) case {'linear','bilinear'} disp('Method is linear') case 'cubic' disp('Method is cubic') case 'nearest' disp('Method is nearest') otherwise disp('Unknown method.') end

1.5 Manfaatkan fasilitas help! Masih sangat banyak sekali bahasan MATLAB yang berlum tercakup dalam buku ajar ini. Semua hal yang berkaitan dengan operasional MATLAB sudah ada dalam help MATLAB. Kita tinggal membuka dan mempelajarinya sendiri.

Apabila

menemukan

kesulitan

dalam

melakukan

pemrograman

menggunakan MATLAB, kita dapat memanfaatkan fasilitas help. Caranya dengan mengetikkan help kemudian ketikan topik yang kita cari.

Halaman 25 dari 101


Tugas 1: Pengenalan MATLAB dan Membuat Program Sederhana Nomor 1 (Tutorial MATLAB) Baca tutorial “Cepat Mahir MATLAB”, Bab 1 Memulai Menggunakan MATLAB dan Bab 5 Fungsi M-File. Nomor 2 (Persamaan Antoine) Buat sebuah algoritma dan program dalam M-file untuk menghitung tekanan uap murni n-heksana dalam rentang temperatur 25 - 100 oC, dengan menggunakan persamaan Antoine sbb:

ln P = A − B /(T + C ) Dengan A = 14.0568 T = Temperatur (K) B = 2825.42

P = Tekanan uap murni (kPa)

C = -42.7089 Buat pula grafik P terhadap T-nya menggunakan rutin plot dalam MATLAB. Nomor 3 (Equimolar Counterdiffusion) Gas amoniak (A) berdifusi melalui pipa sepanjang 0,10 m yang berisi gas N2 (B) pada tekanan 1,0132 x 105 Pa dan temperatur 298 K. Tekanan pada titik 1 PA,1 = 1,013 x 104 Pa dan pada titik 2 PA,2 = 0,507 x 104 Pa. Diffusivitas DAB = 0,230 x 10-4 m2/s. Laju diffusi gas amoniak (A) dapat dievaluasi menggunakan Hukum Fick’s berikut ini: J A* =

DAB ( PA1 − PA 2 ) [kmol. A /( s.m 2 )] RT Δz

R = 8314 J/(kmol.K)

Buat sebuah algoritma dan program MATLAB berupa suatu fungsi dalam M-file untuk menghitung laju diffusi gas amoniak. Petunjuk : program terdiri atas 2 buah m-file. 1 buah untuk menulis fungsi, 1 buah untuk mengeksekusi fungsi.

_____________________________________o0o_____________________________________

Halaman 26 dari 101

Bab 2 Sistem Persamaan Linier

Bab 2 Sistem Persamaan Linier 2.1 Definisi Vektor-Vektor TSL (Terhubung Secara Linier) Sehimpunan vektor x1, x2, ......, xn yang masing-masing berukuran m disebut terhubung secara linier (TSL) jika dapat ditemukan sehimpunan konstanta k1 , k2 ,...., kn dengan harga tidak semuanya nol, sehingga:

k1 x1 + k2 x2 + .... + kn xn = 0 Catatan: 0 adalah vektor yang semua elemennya nol berukuran m. Contoh himpunan vektor TSL ⎡1 ⎤ ⎡2⎤ ⎡3⎤ ⎢ ⎥ ⎢ ⎥ x1 = ⎢ 2 ⎥ , x2 = ⎢ 4 ⎥ , x3 = ⎢⎢6 ⎥⎥ ⎢⎣ 3 ⎥⎦ ⎢⎣ 6 ⎥⎦ ⎢⎣9 ⎥⎦

⎡1 ⎤ ⎡ 2⎤ ⎡ 3⎤ ⎡ 0 ⎤ ⎢ ⎥ ⎢ ⎥ (−2) ⎢ 2⎥ + (1) ⎢ 4 ⎥ +(0) ⎢⎢6⎥⎥ = ⎢⎢0⎥⎥ ⎢⎣ 3 ⎥⎦ ⎢⎣ 6 ⎥⎦ ⎢⎣9 ⎥⎦ ⎢⎣0⎥⎦ k1

k2

k3

2.2 Definisi Vektor-Vektor TTSL Sehimpunan vektor x1, x2, ......, xn yang masing-masing berukuran m disebut tak terhubung secara linier (TTSL) jika TIDAK dapat ditemukan sehimpunan konstanta seperti pada vektor TSL. Penentuan keterhubungan linier dari vektor-vektor dengan cara yang telah didefinisikan diatas tidaklah praktis karena hanya mengandalkan prinsip cobacoba. Metode praktis yang dapat diterapkan untuk menentukan keterhubungan linier dari vektor-vektor adalah ortogonalisasi Gram Schmidt.

Halaman 27 dari 101


2.3 Konsep Keortogonalan Sehimpunan dari n buah vektor berdimensi m (x1, x2, ......, xn) disebut himpunan ortogonal/saling tegak lurus jika kedua syarat di bawah ini terpenuhi.

( xi ) x j T ( xi ) x j T

= 0 untuk i ≠ j dan ≠ 0 untuk i = j dengan i, j = 1, 2,..., n

Sehimpunan vektor-vektor yang ortogonal pasti TTSL, namun TTSL belum tentu ortogonal. Diagram vennya sbb:

S

TTSL Ortogonal

Contoh: ⎡3⎤ ⎡0.6 ⎤ (1) x1 = ⎢ ⎥ , x2 = ⎢ ⎥ ⎣ −1⎦ ⎣1.8 ⎦

⎡ 0.6 ⎤ x1T x2 = [3 −1] ⎢ ⎥ = 1.8 − 1.8 = 0 ⎣1.8 ⎦ ⎡3⎤ x1T x1 = [3 −1] ⎢ ⎥ = 9 + 1 = 10 ⎣ −1⎦

ortogonal

Pada praktek komputasi seringkali ditemukan

( xi )

T

x j ≈ 0 hampir ortogonal.

Untuk ( xi ) x j ≤ 10−5 secara praktis bisa dikatakan ortogonal. T

Ortogonalisasi Gram-Schmidt

Menegakkan sehimpunan vektor-vektor ortogonal y1,y2,....,yn dari suatu himpunan vektor-vektor x1, x2,....,xn.

Halaman 28 dari 101


1. y1 = x1 ⎡ ( y1 )T x2 ⎤ 2. y2 = x2 − ⎢ ⎥ y1 T ⎣⎢ ( y1 ) y1 ⎥⎦ ⎡ ( y1 )T x3 ⎤ ⎡ ( y2 )T x3 ⎤ 3. y3 = x3 − ⎢ ⎥ y1 − ⎢ ⎥ y2 T T ⎣⎢ ( y1 ) y1 ⎦⎥ ⎣⎢ ( y2 ) y2 ⎥⎦ ⎡ ( y1 )T xi ⎤ ⎡ ( yi −1 )T xi ⎤ i. yi = xi − ⎢ ⎥ y1 − KK − ⎢ ⎥ yi −1 T T ⎣⎢ ( y1 ) y1 ⎦⎥ ⎣⎢ ( yi −1 ) yi −1 ⎦⎥ Vektor-vektor y1 , y2 ,K , yn saling tegak lurus satu sama lain dan y1 juga tegak lurus pada x1 , x2 ,K , x j −1 .

Pengujian Keterhubungan Linier Vektor-Vektor

Konstruksi vektor-vektor yang saling tegak lurus (ortogonal) dengan metode Gram-Schmidt dari vektor-vektor yang akan diuji. Jika (V j ) V j = 0 (atau T

mendekati nol), maka himpunan vektor yang diuji TSL.

Contoh: Ortogonalisasikan vektor-vektor kolom dari matrik A berikut ini. ⎡1 2 3⎤ A = ⎢⎢ 4 5 6 ⎥⎥ ⎢⎣7 8 9 ⎥⎦

Jawaban: x1 x2 x3 ⎡1 2 3⎤ A = ⎢⎢ 4 5 6 ⎥⎥ ⎢⎣7 8 9 ⎥⎦

Menghitung y1 y1 = x1 = [1 4 7 ]

T

Halaman 29 dari 101


Menghitung y2

( y1 )

⎡1 ⎤ y1 = [1 4 7 ] ⎢⎢ 4 ⎥⎥ = 1 + 16 + 49 = 66 ⎢⎣ 7 ⎥⎦

( y1 )

⎡ 2⎤ x2 = [1 4 7 ] ⎢⎢ 5 ⎥⎥ = 2 + 20 + 56 = 78 ⎢⎣ 8 ⎥⎦

T

T

⎡ ( y )T x ⎤ 78 T T y2 = x2 − ⎢ 1 T 2 ⎥ y1 = [ 2 5 8] − [1 4 7 ] = [ 0.8182 0.2727 -0.2727 ] 66 ⎢⎣ ( y1 ) y1 ⎥⎦

Menghitung y3

( y1 )

T

⎡3⎤ x3 = [1 4 7 ] ⎢⎢ 6 ⎥⎥ = 3 + 24 + 63 = 90 ⎢⎣9 ⎥⎦

( y2 )

⎡3⎤ x3 = [ 0.8182 0.2727 -0.2727 ] ⎢⎢6 ⎥⎥ = 1.6364 ⎢⎣9 ⎥⎦

( y2 )

⎡ 0.8182 ⎤ y2 = [ 0.8182 0.2727 -0.2727 ] ⎢⎢ 0.2727 ⎥⎥ = 0.8182 ⎢⎣-0.2727 ⎥⎦

T

T

⎡ ( y1 )T x3 ⎤ ⎡ ( y2 )T x3 ⎤ y3 = x3 − ⎢ ⎥ y1 − ⎢ ⎥ y2 T T ⎢⎣ ( y1 ) y1 ⎥⎦ ⎢⎣ ( y2 ) y2 ⎥⎦ 90 1.6364 T y3 = [3 6 9] − [1 4 7 ] − [ 0.8182 0.2727 −0.2727] = [ 0 0 0] 66 0.8182 Dapat disimpulkan bahwa banyaknya vektor-vektor kolom A yang TTSL adalah 2.

Halaman 30 dari 101


2.4 Norma vektor dan matrik

Misalkan x adalah vektor kolom/baris berdimensi n A adalah matrik bujursangkar berdimensi n x n

Norma ke-1 n

x 1 = ∑ xi = x1 + x2 + ... + xn i =1

n

A 1 = maks ∑ aij j

(norma kolom)

i =1

Norma ke-2 1

⎡ n 2⎤ x 2 = ⎢ ∑ ( xi ) ⎥ ⎣ i =1 ⎦ ⇒ Panjang Vektor

2

A 2 ≡ Norma spektral

Memiliki sifat yang lebih baik daripada A 1 dan

A ∞ , sayangnya tidak mempunyai ungkapan matematis yang sederhana sehingga

jarang digunakan dalam komputasi numeris sebagai gantinya digunakan norma frobenius.

Norma Frobenius

A

F

⎡ n n 2⎤ = ⎢ ∑∑ ( aij ) ⎥ ⎣ i =1 j =1 ⎦

1

2

Norma ke- ∞

x

∞

= maks xi i

n

A ∞ = maks ∑ aij i

j =1

⇒ norma baris

Halaman 31 dari 101


Contoh: x = [1 −2 3] 3

x 1 = ∑ xi = x1 + x2 + x3 = 1 + −2 + 3 = 6 i =1

1

2 ⎡ 3 2⎤ x 2 = ⎢ ∑ ( xi ) ⎥ = x12 + x22 + x32 = (1) 2 + (−2) 2 + (3) 2 = 1 + 4 + 9 = 14 ⎣ i =1 ⎦

x

∞

= maks xi = 3 = 3 i

⎡ 1 −2 3⎤ A = ⎢⎢ −4 5 6 ⎥⎥ ⎢⎣ 7 −8 9 ⎥⎦ 1 + −4 + 7 = 12

1 + −2 + 3 = 6 −4 + 5 + 6 = 15

A

∞

= 24

7 + −8 + 9 = 24

A F = 1+ 4 + 9 +16 + 25 + 36 + 49 + 64 + 81 = 285

−2 + 5 + −8 = 15 3 + 6 + 9 = 18 A 1 = 18

2.5 Martabat matrik

Martabat sebuah matrik adalah jumlah maksimum kolom-kolom TTSL dari matrik yang bersangkutan. Banyak maksimum kolom TTSL sama dengan banyak maksimum baris TTSL. Perintah dalam MATLAB untuk menentukan martabat suatu matrik adalah rank. Contoh: >> A=[1 2 3;4 5 6;7 8 9] A = 1

2

3

4

5

6

7

8

9

Halaman 32 dari 101

Bab 2 Sistem Persamaan Linier >> rank(A) ans = 2

Tugas 2 Ortogonalisasi Gram-Schmidt (pilihlah salah satu soal berikut ini!)

Nomor 1 Konstruksi sekumpulan vektor kolom yang saling tegak lurus dari vektor-vektor kolom matrik A di bawah ini, dengan menggunakan ortogonalisasi Gram-Schmdt. ⎡ −1 0 −1 −1 −1⎤ ⎢ −1 −1 0 −2 2 ⎥ ⎢ ⎥ A = ⎢ 1 −1 2 0 4 ⎥ ⎢ ⎥ ⎢3 1 2 4 0⎥ ⎢⎣ 0 1 −1 1 −3⎥⎦ Nomor 2 Buatlah

sebuah

pemrograman

program

MATLAB.

ortogonalisasi Jangan

lupa

Gram-Schmidt untuk

dalam

menyertakan

bahasa algoritma

pemrogramannya.

2.6 Metode eliminasi Gauss

3x1 + 18 x2 + 9 x3 = 18 2 x1 + 3x2 + 3x3 = 117 4 x1 +

x2 + 2 x3 = 283

Matrik perbesarannya adalah sebagai berikut:

⎡ 3 18 9 | 18 ⎤ A = ⎢⎢ 2 3 3 |117 ⎥⎥ ⎢⎣ 4 1 2 |283⎥⎦ Langkah 1 Buat elemen di bawah a11 menjadi nol Hitung: l21 =

a21 2 = a11 3

Halaman 33 dari 101


l31 =

a31 4 = a11 3

Hitung: Baris ke-2 baru = baris ke-2 lama – l21 x baris ke-1 = [ 2 3 3 |117 ] − 2 3 x [3 18 9|18] = [ 0 −9 −3 |105] Baris ke-3 baru = baris ke-3 lama – l31 x baris ke-1

= [ 4 1 2 | 283] − 4 3 x [3 18 9 |18] = [ 0 −23 −10 | 259]

9 | 18 ⎤ ⎡ 3 18 ⎢ A = ⎢0 −9 −3 | 105 ⎥⎥ ⎢⎣0 −23 −10 | 259 ⎥⎦ Langkah 2 Hitung:

l32 =

a32 23 = a22 9

Hitung: Baris ke-3 baru = baris ke-3 lama – l32 x baris ke-2 = [ 0 −23 −10 |105] − 23 9 x [ 0 −9 −3 | 259] = [ 0 0 − 7 3 | − 28 3 ]

9 | 18 ⎤ ⎡ 3 18 ⎢0 −9 −3 | 105 ⎥ ⎢ ⎥ ⎢⎣0 0 −7 / 3 | −28 / 3⎥⎦ 3 x1 + 18 x2 + 9 x3 = 18 − 9 x2 − 3 x3 = 105 − 73 x3 = − 283

x3 = 4

Halaman 34 dari 101


x2 = −13 x1 = 72 Berikut ini permrograman metode Gauss dalam bahasa MATLAB 7 function x = Gauss(A , c) %GAUSS Solves a set of linier algebraic equations by the Gauss % elimination method. % GAUSS(A,C) finds unknowns of a set of linier algebraic % equations. A is the matrix of coefficients and C is the % vector of constants. % % See also JORDAN, JACOBI. %(c) by N. Mostoufi & A. Constantinides %January 1,1999 c = (c(:).’)’

; %Make sure it's a column vector

n = length(c); [nr nc] = size(A); % Check coefficient matrix and vector of constants if nr ~= nc error('Coefficient matrix is not square.') end if nr ~= n error('Coefficient matrix and vector of constants do not have the same length') end % Check if the coefficient matrix is singular if det(A) == 0 fprintf('\n Rank = %7.3g\n',rank(A)) error('The coefficient matrix is singular.') end unit = eye(n); % Unit matrix order = [1 : n]; % Order of unknowns aug = [A c]; % Augmented matrix % Gauss elimination for k = 1 : n-1 pivot = abs(aug(k , k)); prow = k; pcol = k;

Halaman 35 dari 101


% Locating the maximum pivot element for row = k : n for col = k : n if abs(aug(row , col)) > pivot pivot = abs(aug(row , col)); prow = row; pcol = col; end end end % Interchanging the rows pr = unit; tmp = pr(k , :); pr(k , :) = pr(prow , :); pr(prow , :) = tmp; aug = pr * aug; % Interchanging the columns pc = unit; tmp = pc(k , :); pc(k , :) = pc(pcol , :) ; pc(pcol , :) = tmp; aug(1 : n, 1 : n) = aug(1 : n , 1 : n) * pc; order = order * pc; % Keep track of the column interchanges % Reducing the elements below diagonal to zero in the column k lk = unit; for m = k + 1 : n lk(m , k) = - aug(m , k) / aug(k , k); end aug = lk * aug; end x = zeros(n , 1); % Back substitution t(n) = aug(n , n + 1) / aug(n , n); x(order(n)) = t(n); for k = n - 1 : -1 : 1 t(k) = (aug(k,n+1) - sum(aug(k,k+1:n).*t(k+1:n))) / aug(k,k); x(order(k)) = t(k); end

Kasus1 Kukus lewat jenuh bertemperatur 130 oC mengalir dalam sebuah pipa yang memiliki diameter dalam 20 mm (D1), dan diameter luar 25 mm (D2). Pipa diinsulasi setebal 40 mm [(D3 – D2)/2]. Koefisien konveksi kukus (hi) = 1700

Halaman 36 dari 101


W/m2.K, dan koefisien konveksi udara (ho) = 3 W/m2.K. Konduktivitas termal pipa (ks) = 45 W/m.K, dan insulasi (ki) = 0,064 W/m.K. Temperatur udara di luar insulasi = 25 oC. Perkirakan temperatur T1, T2, dan T3.

T2

T3

T1 Kukus, TS

Udara, Ta

Perpindahan panas dari kukus ke pipa.

hiπ D1 (TS − T1 ) =

(T1 − T2 ) ln ( D2 / D1 ) / ( 2π k s )

Perpindahan panas dari pipa ke insulasi

(T1 − T2 ) (T2 − T3 ) = ln ( D2 / D1 ) / ( 2π k s ) ln ( D3 / D2 ) / ( 2π ki ) Perpindahan panas dari insulasi ke udara

(T2 − T3 ) = hOπ D3 (T3 − Ta ) ln ( D3 / D2 ) / ( 2π ki ) Ada tiga persamaan linier yang berhasil dirumuskan dari peneracaan energi tersebut. ⎡ ⎤ ⎡ ⎤ 2k s 2k s + hi D1 ⎥ T1 − ⎢ ⎢ ⎥ T2 = hi D1TS ln D / D ln D / D ( ) ( ) 2 1 2 1 ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ks ks ki ki + ⎢ ⎥ T1 − ⎢ ⎥ T2 + ⎢ ⎥ T3 = 0 ⎣ ln ( D2 / D1 ) ⎦ ⎣ ln ( D2 / D1 ) ln ( D3 / D2 ) ⎦ ⎣ ln ( D3 / D2 ) ⎦ ⎡ ⎤ ⎡ ⎤ 2 ki 2 ki + hO D3 ⎥ T3 = − hO D3Ta ⎢ ⎥ T2 − ⎢ ⎣ ln ( D3 / D2 ) ⎦ ⎣ ln ( D3 / D2 ) ⎦

Ubah sistem persamaan linier menjadi bentuk matrik Ax = c, menjadi sbb:

Halaman 37 dari 101


⎡⎡ ⎤ ⎤ ⎡ ⎤ 2k s 2k s + hi D1 ⎥ −⎢ 0 ⎢⎢ ⎥ ⎥ ⎢ ⎣ ln ( D2 / D1 ) ⎥ ⎦ ⎣ ln ( D2 / D1 ) ⎦ ⎢ ⎥ ⎡ T1 ⎤ ⎡ hi D1TS ⎤ ⎤ ⎡ ⎤ ⎡ ⎤ ks ks ki ki ⎢ ⎡ ⎥ ⎢T ⎥ = ⎢ 0 ⎥ −⎢ + ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥⎢ 2⎥ ⎢ D D D D D D D D ln / ln / ln( / ) ln / ( ( ( 2 1)⎦ 2 1) 3 2 ⎦ 3 2 )⎦ ⎣ ⎣ ⎣ ⎢ ⎥ ⎢⎣T3 ⎥⎦ ⎢⎣ −hO D3Ta ⎥⎦ ⎢ ⎡ ⎤ ⎡ ⎤⎥ 2 ki 2 ki h D 0 − + ⎢ ⎢ ⎥ ⎢ O 3 ⎥⎥ D D D D ln / ln / ( ) ( ) ⎢⎣ 3 2 3 2 ⎣ ⎦ ⎣ ⎦ ⎥⎦ Berikut ini pemrograman MATLAB-nya. %kasus2.m clc clear % Input data Ts = 130; % oC Ta = 25; % oC D1 = 20e-3; % Diameter dalam pipa, m D2 = 25e-3; % Diameter luar pipa, m Ith = 40e-3; % Tebal insulasi, m D3 = (D2 + 2*Ith); % Diameter pipa + insulasi hi = 1700; % Koefisien transfer panas bagian dalam (W/m2.K) ho = 3 ; % koefisien transfer panas bagian luar (W/m2.K) ks = 45; % Konduktivitas panas baja (W/m.K) ki = 0.064; % Konduktivitas panas insulasi (W/m.K) % Matriks koefisien variabel A = [2*ks/log(D2/D1)+hi*D1 , -2*ks/log(D2/D1) , 0 ks/log(D2/D1) , -(ks/log(D2/D1)+ki/log(D3/D2)) , ki/log(D3/D2) 0 , 2*ki/log(D3/D2) , -(2*ki/log(D3/D2)+ho*D3)]; % Matriks konstanta c = [hi*D1*Ts ; 0 ; -ho*D3*Ta]; % Menyelesaikan sis pers. linier dengan fungsi invers MATLAB T = inv(A)*c

Eksekusi persamaan di command window >>kasus2 T= 129.7858 129.7678 48.1191

Halaman 38 dari 101


Tugas 3: Menyelesaikan Sistem Persamaan Aljabar Linier Secara Simultan Suatu Sistem Distribusi Uap Dalam Sebuah Pabrik Kimia. Nomor 1

Sebuah sistem persamaan linier dirumuskan dari Neraca massa & energi distribusi uap pabrik (ditampilkan di bawah). Sistem tersebut terdiri dari 14 buah variabel xi dengan i = 3,...,16 belum diketahui, dan yi adalah parameter yang telah diketahui. xi dan yi dalam 1000 lb/h. Dengan menggunakan MATLAB hitunglah 14 variabel (xi, i=3,…,16) yang belum diketahui itu. 181.60 − x3 − 132.57 − x4 − x5 = − y1 − y2 + y5 + y4 = 5.1 1.17 x3 − x6 = 0 132.57 − 0.745 x7 = 61.2 x5 + x7 − x8 − x9 − x10 + x15 = y7 + y8 − y3 = 99.1 x8 + x9 + x10 + x11 − x12 − x13 = − y7 = −8.4 x6 − x15 = y6 − y5 = 24.2 −1.15(181.60) + x3 − x6 + x12 + x16 = 1.15 y1 − y9 + 0.4 = −19.7 181.60 − 4.594 x12 − 0.11x16 = − y1 + 1.0235 y9 + 2.45 = 35.05 −0.0423(181.60) + x11 = 0.0423 y1 = 2.88 −0.016(181.60) + x4 = 0 x8 − 0.147 x16 = 0 x5 − 0.07 x14 = 0 −0.0805(181.60) + x9 = 0 x12 − x14 + x16 = 0.4 − y9 = −97.9

Nomor 2

Ketik ulang program Gauss.m pada m-file MATLAB. Gunakan fungsi gauss itu untuk menyelesaikan sistem persamaan linier pada nomor 1.

________________________________o0o_______________________________

Halaman 39 dari 101

Bab 3 Persamaan Tak Linier


Persamaan Linier

y = mx + c

y

y=x

LINIER

x

Gambar 3.1 Kurva linier Persamaan Tak Linier Contoh:

y y = exp( x)

NON-LINIER Gambar 3.2 Kurva tak linier

x

Halaman 40 dari 101


Berikut ini beberapa contoh Persamaan Tak Linier Tabel 3.1 Contoh Persamaan Tak linier Jenis Pers.

Contoh

Tak Linier

x2 − 4 x + 3 = 0

Persamaan Kuadrat Persamaan Polinomial

x 4 + 6 x3 + 7 x 2 − 6 x − 8 = 0

Persamaan Transenden

sin x − 2 exp(− x 2 ) = 0

Persamaan Logaritmik

ln(1 + x 2 ) − 2 exp(− x 2 ) = 0

Dalam aplikasinya di bidang teknik kimia, persamaan tak linier memiliki peranan yang sangat penting. Tabel 3.2 Aplikasi Persamaan Tak Linier dalam bidang teknik kimia Aplikasi Pers. Tak Linier

Neraca Massa dan Energi,

Contoh

εΔH + N o 0

Tout out

∫C

out P ,i

To

Termodinamika Persamaan gas nyata/kubik, Kesetimbangan reaksi kimia,

Operasi Teknik Kimia, dll.

P=

RT a − 2 V −b V

dT − N

Tin in

∫C

in P ,i

=0

To

(1

ΔC p dT ΔG0o − ΔH 0o ΔH 0o 1 ΔC p ln K + + + ∫ dT − ∫ =0 RT0 RT T T0 R R T T0 T

⎛ n α j z jF F ⎞ ⎜⎜ ∑ ⎟⎟ − F (1 − q ) = 0 − α φ j 1 = j ⎝ ⎠

o

T

o

(2

1) Persamaan kubik tersebut diusulkan oleh Johannes Diderik van der Waals (1873), Fisikawan Belanda, peraih nobel Fisika pada tahun 1910. 2) Persamaan Underwood untuk distilasi multikomponen.

Halaman 41 dari 101


Klasifikasi persamaan tak Linier Berdasarkan jumlah banyaknya persamaan tak linier dibagi menjadi dua 1. Persamaan Tunggal Persamaan tak linier hanya satu buah. Contoh: x 2 + 2 x − 3 = 0

f ( x) = 0

2. Persamaan Serentak/Sistem Persamaan Persamaan tak linier terdiri atas minimal dua buah.

f 1 ( x1 , x 2 ,..., x N ) = 0

f 2 ( x1 , x 2 ,..., x N ) = 0 ... f N ( x1 , x 2 ,..., x N ) = 0

Contoh :

2 x − sin{( x + y ) / 2} = 0 2 x − cos{( x − y ) / 2} = 0

Solusi Persamaan Tak Linier Tunggal

Mencari akar-akar x yang membuat harga y atau f(x) menjadi nol.

f ( x ) = 0 → x = .... ? Ada berbagai metode numerik yang dapat digunakan untuk menyelesaikan persaman tak linier tunggal, diantaranya: z Metode Penyetengahan Interval z Metode Substitusi Berurut z Metode Wegstein z Metode Interpolasi Linear z Metode Newton-Raphson,dll.

Dalam diktat ini hanya akan diterangkan metode penyetengahan interval dan metode Newton-Raphson..

Halaman 42 dari 101


Metode penyetengahan interval

Tabel 3.3 Karakteristik metode penyetengahan interval No 1.

Keunggulan Sederhana

Kelemahan Tebakan awal terdiri atas dua buah [a,b] dan harus memenuhi f(a)*f(b)<0

2.

Pasti konvergen

Laju konvergensi relatif lebih lambat daripada metode Newton-Raphson

Algoritma metode penyetengahan interval

1. Mulai 2. Definisikan persamaan tak linier yang akan dicari akarnya dan tetapkan toleransinya. Misalkan: f(x) = log(x) . Tol = 1e-5 3. Tetapkan dua buah tebakan awal a dan b. Misalkan: a = 0.1 dan b = 10. 4. Hitung harga f(a) dan f(b) Hitung: f(0.1) = log(0.1) = -1 f(10) = log(10) = 1 5. Periksa apakah syarat f(a) * f(b) < 0 terpenuhi. Jika syarat terpenuhi proses dilanjutkan ke langkah berikutnya. Jika tidak kembali ke langkah 2. Periksa: f(0.1) * f(10) = -1*1 = -1 < 0. Syarat terpenuhi proses berlanjut ke langkah berikutnya. 6. Hitung harga m = (a + b) / 2 Hitung: m = (0.1 + 10)/2 = 5.05 7. Hitung harga f(m) Hitung: f(5.05) = log(5.05) = 0.7033 8. Periksa apakah f(a) * f(m) > 0 terpenuhi. Jika terpenuhi, harga a diganti dengan m. Jika tidak terpenuhi harga b yang diganti dengan m. Periksa: f(0.1) * f(5.05) = -1*0.7033 = -0.7033 < 0. Syarat tidak terpenuhi harga b diganti dengan m, sehingga b yang baru = 5.05.

Halaman 43 dari 101


9. Periksa kriteria iterasi, |(a - b)/a| > Tol. Jika kriteria iterasi terpenuhi proses kembali ke langkah 5. Jika tidak dilanjutkan ke langkah berikutnya.. Periksa: |(0.1 – 5.05)/0.1| = 49.5 > 1e-5. Kriteria iterasi terpenuhi, maka kembali ke langkah 5. 10. Hitung harga akar pembuat nol, x = (a + b) / 2 11. Selesai Algoritma yang telah disusun dimuka, dapat juga dituliskan dalam bentuk diagram Alir (flow chart) sebagai berikut: 1

mulai

m=(a+b)/2

Definisikan f(x) dan toleransi

Hitung harga: f(m)

Tetapkan harga a dan b

f(a)*f(m)

y Hitung harga: f(a), f(b)

tid

a=m f(a)=f(m)

b=m f(b)=f(m)

ya

tidak

|(a-b)/a|>tol f(a)*f(b)< 0

ya

tidak x*=(a+b)/2

1 Selesai

Gambar 3.3 Algoritma metode penyengahan interval

Halaman 44 dari 101


Pemrograman metode penyetengahan interval dengan menggunakan bahasa MATLAB. biseksi.m function x = biseksi(fungsi,a,b,tol) % a = tebakan awal pertama, b = tebakan awal kedua % tol = toleransi while abs((a - b)/a) > tol fa = feval(fungsi,a); fb = feval(fungsi,b); if fa*fb > 0 error('masukan tebakan a dan b yang berbeda') end m = (a + b)/2; fm = feval(fungsi,m); if fm*fa > 0; a = m; else b = m; end end x=(a+b)/2;

Contoh: Carilah akar-akar persamaan kuadrat x 2 + 4 x + 3 = 0 dengan menggunakan metode penyetengahan interval!. %kuadrat.m function y = kuadrat(x) y = x^2+4*x+3;

Perintah pada command window sbb: >>biseksi(‘kuadrat’,-2,1,1e-6) ans = -1.0000 >> biseksi('kuadrat',-2,-4,1e-6) ans = -3.0000

Dari perhitungan menggunakan metode bisection diperoleh akar-akar dari persamaan kuadrat adalah [-1,-3].

Halaman 45 dari 101


Metode Newton-Raphson

Tabel 4.3 Karakteristik metode penyetengahan interval No 1. 2.

Keunggulan

Kelemahan

Hanya butuh satu tebakan

Kekonvergenan adakalanya gagal

awal.

dicapai.

Laju konvergensi cepat

y f ( xn ) 0

xn

xn +1

x

Gambar 3.4 Metode newton-Raphson

f '( xn ) = gradien = xn +1 = xn −

f ( xn +1 ) − f ( xn ) 0 − f ( xn ) = xn +1 − xn xn +1 − xn

f ( xn ) f '( xn )

Metode Newton-Raphson.

f ( xn ) xn +1 = xn − f '( xn )

Algoritma Metode Newton-Raphson

1. Mulai 2. Definisikan persamaan tak linier dan turunannya. Misalkan: f(x) = log(x).

f ' ( x) = 1 /( x ln 10) 3. Tetapkan harga tebakan awal ( x0 ) dan besar toleransinya Misalkan: x0 = 0.1 . Tol = 1e-5 4. Nyatakan x = x 0 dan x0 = x + 1 .

Halaman 46 dari 101


x = 0.1 x0=0.1+1=1.1 5. Periksa kriteria iterasi |(x – x0)/x| > Tol. Jika kriteria iterasi terpenuhi proses dilanjutkan. Jika kriteria iterasi tidak terpenuhi proses dihentikan. Akar pembuat nol diperoleh. 6. Nyatakan x0 = x.

x0 = 0.1 7. Hitung harga f(x0) dan f’(x0).

f(0.1) = -1 f’(0.1) = 1/(0.1*ln10) = 4.343 8. Hitung harga x = x0 −

x = 0.1 −

f ( x0 ) f ' ( x0 )

(−1) = 0.3303 4.343

9. Kembali ke langkah 5 10. Selesai.

Halaman 47 dari 101


Mulai

Definisikan f(x) dan f’(x), x0, tol Nyatakan: x = x0 x0 = x + 1

tidak ya Nyatakan: x0 = x

Tampilkan x* = x

Hitung harga: f(x0) dan f’(x0)

Selesai

Hitung harga: x=x0-f(x0)/f’(x0)

Gambar 3.5 Algoritma metode Newton-Raphson

Halaman 48 dari 101


%NewRap.m function [x iter] = NewRap(fungsi,dfungsi,x0,tol) % fungsi = fungsi yang akan dicari akar-akarnya % dfungsi = turunan pertama fungsi % x0 = tebakan awal % tol = toleransi itermax = 100; iter = 0; x = x0; x0 = x + 1; % loop iterasi while abs((x - x0)/x) > tol & iter <= itermax iter = iter + 1; x0 = x; fx= feval(fungsi,x); df= feval(dfungsi,x); % Rumus Newton-Raphson x = x0 - fx/df; end

%kuadrat.m function y = kuadrat(x) y = x^2+4*x+3;

%dkuadrat.m function dy = dkuadrat(x) dy = 2*x+4;

Kita gunakan contoh kasus yang sama dengan contoh pada metode bisection. y = x2 + 4x + 3 = 0 dy = 2x + 4 dx >> [x iter]=NewRap('kuadrat','dkuadrat',2,1e-6) x = -1.0000

Halaman 49 dari 101

Bab 3 Persamaan Tak Linier iter = 6 >> [x iter]=NewRap('kuadrat','dkuadrat',-4,1e-6) x = -3.0000 iter = 5

Dari perhitungan menggunakan metode Newton Raphson diperoleh akar-akar dari persamaan kuadrat adalah [-1,-3].

Subrutin dalam MATLAB untuk persamaan tak linier tunggal

MATLAB telah menyediakan program untuk menyelesaikan persamaan linier tunggal yang telah menyatu dengan program MATLAB itu sendiri. Ada dua subrutin yang umum digunakan, yaitu roots dan fzero. Tabel 4.4 Perbandingan subrutin roots terhadap fzero Rutin

Keunggulan

roots.m

1. Seluruh akar dapat diketahui

Kelemahan 1. Hanya untuk pers. kuadrat

dengan hanya sekali

dan polinomial.

menjalankan rutin. 2. Tidak membutuhkan tebakan mula. fzero.m

1. Solusi bagi segala jenis pers

1. Hanya satu buah akar

tak linier.

yang dapat diketahui sekali menjalankan rutin. 2. Membutuhkan tebakan mula.

Penggunaan roots: Penulisan perintah roots di Command window MATLAB C(1)*X^N + ... + C(N)*X + C(N+1)

Halaman 50 dari 101


C = [C(1) C(2)........C(N) C(N+1) roots(C) Contoh persamaan kuadrat x 2 + 4 x − 5 = 0 maka C(1)=1, C(2)=4, C(3)= -5. Carilah akar-akar persamaan kuadrat di bawah ini.

x2 + 4 x − 5 = 0 MATLAB Command window >> C=[1 4 -5] C = 1

4

-5

>> roots(C) ans = -5 1

Penggunaan fzero: Penulisan fzero di MATLAB Command window

x =fzero(‘fungsi’,x0) Contoh penggunaan fzero: x2 + 4 x + 3 = 0 Penulisan contoh di MATLAB Command window >> fzero('x^2+4*x+3',0) ans = -1

Halaman 51 dari 101


Untuk keteraturan dan kemudahan pemanggilan akan lebih baik mendefinisikan fungsi pada m-file. %kuadrat.m function y = kuadrat(x) y = x^2+4*x+3

Baru kemudian kita panggil fungsi dari MATLAB Command window >> x = fzero('kuadrat',0) x = -1

Untuk mencari akar lainnya, ubah tebakan awalnya. >> x = fzero('kuadrat',-4) x = -3.0000

Kasus 3 Aplikasi subrutin roots Kasus 3 Tekanan uap n-butana pada temperatur 350 K adalah 9.4573 bar. Hitunglah volume molar uap jenuh dan cair jenuh n-butana pada Kondisi tersebut dengan menggunakan persamaan gas Van der Waals. (R=8.314j/mol.K ;Tc=425.1 K; Pc=37.96 bar)

Jawaban: Persamaan Van der Waals P=

RT a − 2 V −b V

a=

27 R 2Tc 2 1 RTc dan b = 64 Pc 8 Pc

Transformasi ke dalam bentuk umum pers.polinomial

Halaman 52 dari 101


P (V − b)(V 2 ) = RTV 2 − a(V − b) P (V 3 − bV 2 ) = RTV 2 − aV + ab PV 3 − ( Pb + RT )V 2 + aV − ab = 0 % kasus3.m clear clc % Masukan kondisi operasi P = input('masukan tekanan, Pa = '); T = input('masukan temperatur, K = '); R = 8314 ; %J/(kmol.K) Pc = 37.96e5; %Pa Tc = 425.1; %K % Hitung konstanta a & b a = (27/64)*R^2*Tc^2/Pc; b = (1/8)*R*Tc/Pc; % Definisikan koefisien polinomial VdW=[P, -(P*b + R*T), a, -a*b]; vol = roots(VdW) %liter/mol

Eksekusi program kasus3.m. Masukan dan hasil di Command Window >>kasus3 masukan tekanan, Pa = 9.4573e5 masukan temperatur, K = 350 vol = 2.6669 0.3354 0.1910

Kasus 4 Aplikasi subrutin fzero Diketahui sebuah persamaan kapasitas panas sbb.

Cp = 0.716 − 4.257 E −6T +

15.04 T

⎡ kJ ⎤ ⎢ kg.K ⎥ ⎣ ⎦

Tentukan temperatur pada saat Cp = 1 kJ/kg.K !

Halaman 53 dari 101


Langkah 1 Membuat program fungsi yang akan dinolkan. %KapPns.m function f = KapPns(T,cp) %Persamaan tak linier yang akan dinolkan f = cp - 0.716 + 4257e-6*T - 15.04/T^0.5;

Langkah 2 Membuat program pengeksekusi % kasus4.m clear clc cp = input('masukan kapasitas panas,kJ/kg.K = '); T = fzero(@(T) KapPns(T,cp),100)

Langkah 3 Eksekusi program kasus4.m Masukan dan hasil di Command Window >> kasus4 masukan harga kapasitas panas,kJ/kg.K = 1 T= 189.7597

Tugas 4 Menyelesaikan persamaan tak linier tunggal dengan menggunakan subrutin MATLAB

Tekanan uap n-butana pada temperatur 350 K adalah 9.4573 bar. Volume molar uap jenuh dan cair jenuh n-butana pada kondisi tersebut dapat dihitung dengan menggunakan persamaan kubik Redlich-Kwong-Soave sebagai berikut: P=

RT aα − V − b V (V + b)

Dalam bentuk persamaan polinomial menjadi sebagai berikut::

Halaman 54 dari 101


Z 3 − Z 2 + ( A − B − B 2 ) Z − AB = 0 Dengan: B=

bP α aP PV A= 2 2 Z = RT RT RT

a=

0.4278R 2TC2 0.0867 RTC ;b = PC PC 2

⎡ ⎛ T ⎞⎤ 2 α = ⎢1 + S ⎜⎜1 − ⎟⎟ ⎥ ; S = 0.48508 + 1.55171ω − 0.15613ω T ⎢⎣ C ⎠⎥ ⎝ ⎦ (R=8.314j/mol.K ;Tc=425.1 K; Pc=37.96 bar; ω = 0.1931). Hitunglah volume molar uap jenuh dan cair jenuh n-butana pada kondisi itu !!.

Sistem Persamaan Tak Linier

Sistem persamaan tak linier merupkan persamaan tak linier yang terdiri atas lebih dari satu buah persamaan tak linier. Solusi Sistem Persamaan Tak Linier

Metode Newton f ( x1 , x2 ) = 0 f ( x1 , x2 ) = 0

⎡ ∂f1 (1) ⎢ ∂x | x ⎢ 1 ⎢ ∂f 2 (1) ⎢ ∂x | x ⎣ 1

∂f1 (1) ⎤ | x ⎥ (1) ∂x2 ⎡δ ⎤ ⎡ f (1) ⎤ ⎥⎢ 1 ⎥ = −⎢ 1 ⎥ (1) ∂f 2 (1) ⎥ ⎣ ∂ (1) ⎣ f2 ⎦ |x ⎥ 2 ⎦ ∂x2 ⎦

Jδ = − f

x ( n +1) = x ( n ) + ρδ

Halaman 55 dari 101


function [xnew , iter] = Newton(fnctn,x0,rho,tol,varargin) %(c) by N. Mostoufi & A. Constantinides,January 1, 1999 %Initialization if nargin < 4 | isempty(tol) tol = 1e-6; end if nargin < 3 | isempty(rho) rho = 1; end x0 = (x0(:).')'; % Make sure it's a column vector nx = length(x0); x = x0*1.1; xnew = x0; iter = 0; maxiter = 100; while max(abs(x-xnew)) > tol & iter < maxiter iter = iter + 1; x = xnew; fnk = feval(fnctn,x,varargin{:}); % Set dx for derivation for k = 1:nx if x(k) ~= 0 dx(k) = x(k) / 100; else dx(k) = 1/100; end end % Calculation of the Jacobian matrix a = x; b = x; for k = 1 : nx a(k) = a(k) - dx(k); fa = feval(fnctn,a,varargin{:}); b(k) = b(k) + dx(k); fb = feval(fnctn,b,varargin{:}); jacob(:,k) = (fb - fa) / (b(k) - a(k)); a(k) = a(k) + dx(k); b(k) = b(k) - dx(k); end % Next approximation of the roots if det(jacob) == 0 xnew = x + max([abs(dx), 1.1*tol]); else xnew = x - rho * inv(jacob) * fnk; end end if iter >= maxiter disp('Warning : Maximum iterations reached.') end

Halaman 56 dari 101


Contoh sistem persamaan tak linier.

f1 ( x, y ) = x 3 − 3 xy 2 − 1/ 2 = 0 f 2 ( x, y ) = 3 x 2 y − y 3 − 3 / 2 = 0 Langkah 1 Buat terlebih dahulu fungsi sistem persamaan taklinier dalam m-file. %sistem.m function f = sistem(x) f=[x(1)^3-3*x(1)*x(2)^2-0.5 3*x(1)^2*x(2)-x(2)^3-sqrt(3)/2]

Langkah 2 Buat program pengeksekusi menggunakan Newton.m pada m-file yang berbeda atau dapat juga langsung di command window. %run_sistem.m [x iter] = Newton('sistem',[1 2])

Langkah 3 Jalankan program pengeksekusi. Klik debug/run Langkah 2 dapat di loncat dengan menuliskan langsung perintah eksekusi pada Command window >> [x iter] = Newton('sistem',[1 2]) x = 2.5198 1.5874

iter = 7

Subrutin dalam MATLAB untuk sistem persamaan taklinier

Solusi sistem persamaan taklinier dapat menggunakan fsolve pada MATLAB. Contoh: x3 − 3 xy 2 = 1/ 2 3x 2 y − y 3 = 3 / 2

Langkah 1 Buat terlebih dahulu fungsi sistem persamaan taklinier dalam m-file.

Halaman 57 dari 101


Langkah 2 Buat program pengeksekusi menggunakan fsolve pada m-file yang berbeda atau dapat juga langsung di command window.

Langkah 3 Jalankan program pengeksekusi. function f = sistem(x) f=[x(1)^3-3*x(1)*x(2)^2-0.5 3*x(1)^2*x(2)-x(2)^3-sqrt(3)/2]

>>[X,FVAL] = fsolve('sistem',[1 2]) Optimization terminated: first-order optimality is less than options.TolFun. X= 2.5198

1.5874

FVAL = 1.0e-010 * 0.1930 0.0966

Kasus 5 Reaksi reformasi kukus berlangsung menurut rangkaian reaksi kesetimbangan berikut:

CH 4( g ) + H 2O( g ) CO( g ) + 3H 2( g )

R-1

CO( g ) + H 2O( g ) CO2( g ) + H 2

R-2

Pada suhu 2000 K harga konstanta kesetimbangan untuk masing-masing reaksi adalah 1,930x10-4 dan 5,528. Tentukan komposisi kesetimbangan komponenkomponen apabila Gas umpan berkomposisi 20% CH4(g) dan 80% H2O(g) berada pada kondisi suhu 2000 K dan tekanan 1 atm.

Halaman 58 dari 101


Jawaban Misal ditetapkan basis perhitungan 10 mol gas umpan e1 = derajat reaksi dari reaksi pertama e2 = derajat reaksi dari reaksi kedua Fraksi mol kesetimbangan setiap komponen dapat dinyatakan sebagai berikut: YCO =

e1 − e2 10 + 2e1

YH 2 =

e2 10 + 2e1

YCH 4 =

YCO2 =

3e1 + e2 10 + 2e1

YH 2O =

8 − e1 − e2 10 + 2e1

2 − e1 10 + 2e1

Persamaan konstanta kesetimbangan dinyatakan sebagai berikut: K1 =

YCOYH32 P 2 YCH 4 YH 2O

K2 =

YCO2 YH 2 YCOYH 2O

( e1 − e2 )( 3e1 − e2 ) 2 ( 2 − e1 )(8 − e1 − e2 )(10 + 2e1 ) 3

= K1

e2 ( 3e1 + e2 ) =K ( e1 − e2 )( 8 − e1 − e2 ) 2

Berikut ini pemrograman MATLAB-nya. function y = KsT(e,K1,K2) %Sistem Pers.tak linier yang akan dinolkan y = [(e(1)-e(2))*(3*e(1)-e(2))^3 /((2-e(1))*(8-e(1)… - e(2))*(10+2*e(1))^2) - K1 e(2)*(3*e(1)+e(2)) / ((e(1)-e(2))*(8-e(1)-e(2))) - K2];

clear clc K1 = input(‘Masukan konstanta kst. reaksi 1 = '); K2 = input(‘Masukan konstanta kst. reaksi 2 = '); %Pencari nol fungsi KsT.m e = fsolve(@(e) KsT(e,K1,K2),[1 0.5])

Halaman 59 dari 101


Eksekusi di MATLAB command window >>kasus5 Masukan harga konstanta kst. reaksi 1 = 1.93e-4 Masukan harga konstanta kst. reaksi 2 = 5.528 Optimization terminated: first-order optimality is less than options.TolFun. e= 0.7480

0.6920

Tugas 5 Menyelesaikan sistem persamaan tak linier dengan menggunakan subrutin MATLAB

Suatu reaksi elementer A →

B + C berlangsung dalam sebuah reaktor tangki

berpengaduk kontinu. Laju umpan murni A, 12 mol/s pada temperatur 25 oC. Reaksi bersifat eksotermik, untuk itu digunakan air pendingin bertemperatur 50 oC untuk menyerap kalor yang dibebaskan reaksi. Asumsi konstanta kapasitas panas sama baik di sisi reaktan maupun produk, neraca energi untuk sistem ini dirumuskan sebagai berikut: − FAo X ΔH R = FAoC P , A (T − T0 ) + UA(T − Ta )

FA0 = laju molar umpan, mol/s. X

= konversi

∆HR = Kalor reaksi, J/(mol.K)

CP,A = kapasitas panas A, J/(mol.K) T

= temperatur reaktor, oC

T0 = temperatur referensi, 25 oC Ta = temperatur air pendingin, oC U

= koefisien pindah panas total, W/(m2.K)

Halaman 60 dari 101


A

= luas pindah panas, m2

Untuk reaksi orde pertama konversi dirumuskan sebagai berikut: X=

τk 1+τ k

Dengan τadalah waktu tinggal dalam sekon, dan k adalah laju reaksi spesifik dalam s-1 dihitung dengan menggunakan persamaan Arrhenius: k = 650 exp[ −3800 /(T + 273)]

Hitunglah harga temperatur reaktor dan konversinya!. (∆HR=-1500 kJ/mol; τ=10 s; CP,A = 4500 J/(mol.K); UA/FA0 =700 W.s/(mol.K).

_______________________________o0o________________________________

Halaman 61 dari 101

Bab 4 Optimisasi

Bab 4 Optimisasi

Titik maksimum 5

F(X,Y )

0

-5

-10 2

3

Titik 0 minimum

2 1 0 -1

-2 Y

-2 -3

X

Optimisasi Menggunakan MATLAB Untuk mencari harga minimum dan maksimum kita dapat menggunakan perintah fminsearch. Berikut ini cara penulisannya.

[x,fval,exitflag] = fminsearch(fun,x0) keterangan: fun = Fungsi yang akan diminimumkan atau dimaksimumkan x0 = Tebakan awal x = Harga x yang menyebabkan fungsi minimum atau maksimum fval = Nilai maksimum atau minimum. exitflag = Kriteria penghentian proses iterasi. Harga x mencapai kekonvergenan jika exitflag bernilai 1.

Halaman 62 dari 101

Bab 4 Optimisasi

Optimisasi Variabel Tunggal Carilah titik minimum x 2 + 4 x + 3 = 0 dengan menggunakan subrutin fiminsearch

dalam MATLAB. >> [x fval exitflag]=fminsearch('kuadrat',2) x = -2.0000 fval = -1 exitflag = 1

Optimisasi Variabel Jamak

Kasus12 Carilah titik minimum dari persamaan multivariabel berikut ini. y = ( x1 − 3) 2 + 0.5( x2 − 4) 2 + 3

Jawaban %multivaribel.m function y = multivariabel(x); y = (x(1)-3)^2 + 0.5*(x(2)-4)^2 + 3;

%kasus12.m [x,fval,exitflag] = fminsearch('multivariabel',[1,16])

>> kasus12 x = 3.0000

4.0000

fval = 3.0000 exitflag = 1

Halaman 63 dari 101

Bab 4 Optimisasi

Tugas 11 Optimisasi

Nomor 1: Temperatur Optimal Dalam Reaktor kf ⎯⎯ → 8 -5000/T detikA←⎯ ⎯ P Diselenggarakan dalam reactor batch. Diketahui kf = 10 e kr

1 dan k = 1016e-10000/T detik-1, dimana T dalam K Neraca massa P dalam r reaktor partaian: dCP = k f C A − kr CP = k f ( C Ao − CP ) − kr CP dt dCP + ( k f + kr ) CP = k f C Ao dt

(

− ( k f + k r )t CP k f 1 − e = Hasil integrasi persamaan ini adalah: C Ao k f + kr

)

CA0 = konsentrasi A mula-mula Hitung temperatur optimal yang menyebabkan perolehan maksimal produk P pada waktu reaksi 1 detik.

Nomor 2: Rosenbrock

Carilah titik minimum dari fungsi multidimensional Rosenbrock berikut. f ( x) = 100( x2 − x12 ) 2 + (1 − x1 ) 2

________________________________o0o_______________________________

Halaman 64 dari 101

Bab 5 Regresi Linier dan Non Linier

Bab 5 Regresi Linier dan Non Linier 5.1 Regresi Linier Misalkan kita memiliki sekumpulan data-data sbb ( x1 , y1 ), ( x2 , y2 ),....., ( xn , yn ) Kemudian diperoleh suatu persamaan matematik untuk merepresentasikan datadata eksperimental tersebut berupa persamaan linier. y = a0 x + a1 Jumlah kuadrat terkecil dari selisih antara model dengan data sbb n

n

S = ∑ [ ymodel ( xi ) − yi ] = ∑ [ a0 xi + a1 − yi ] 2

i =1

2

i =1

Turunan terhadap parameter a0 n ∂S = 2∑ xi ( a0 xi + a1 − yi ) = 0 ∂a0 i =1 n

∑ x (a x + a − y ) = 0 i =1

0 i

i

n

1

i

n

n

∑x a +∑xa =∑x y i =1

2 i 0

i =1

i 1

i =1

i

i

............................................(i

Turunan terhadap parameter a1 n ∂S = 2∑ ( xi a0 + a1 − yi ) = 0 ∂a1 i =1 n

∑(a x + a − y ) = 0 i =1

0 i

1

i

n

n

i =1

i =1

∑ xi a0 + na1 = ∑ yi .............................................................(ii Persamaan (i dan (ii dapat dijelmakan dalam bentuk matrik sbb: ⎡ n 2 ⎢ ∑ xi ⎢ i =1 ⎢ n ⎢ ∑ xi ⎣ i =1

⎤ ⎡ n ⎤ x xi yi ⎥ ∑ ∑ i⎥ ⎢ ⎡a ⎤ i =1 ⎥ ⎢ 0 ⎥ = ⎢ i =1 ⎥ ⎥ ⎣ a1 ⎦ ⎢ n ⎥ n ⎥ ⎢ ∑ yi ⎥ ⎦ ⎣ i =1 ⎦ n

Halaman 65 dari 101


Harga paramter a0 dan parameter a1 dapat diperoleh dengan menyelesaikan sistem persamaan linier di atas. Kasus 1 Harga konduktivitas alumunium pada berbagai temperatur sbb T (K)

300

400

500

600

800

k (Btu/(h×ft2)(°F/ft) 273 240 237 232 220 Model matematik dapat diwakili dengan menggunakan persamaan linier k = a0T + a1

Untuk mencari harga a0 dan a1 dapat menggunakan metode jumlah selisih kuadrat terkecil seperti yang telah dijelaskan sebelumnya. %konduktivitas.m clear clc T=[4,5,6,8]*100; k=[240,237,232,220]; n=length(k); A = [sum(T.^2),sum(T) sum(T), n]; c = [sum(k.*T) sum(k)]; a = A\c kmod = a(1)*T +a(2); S = sum((k-kmod).^2)

% Absis % Ordinat

Hasil pada command window >>konduktivitas a = -0.0511 261.6571 S = 3.8857

Subrutin MATLAB untuk regresi persamaan linear dan polinomial dapat menggunakan perintah sbb:

[P,S] = polyfit(x,y,n) data-data

n=orde polinom

Kasus 1 merupakan persamaan linier, maka n yang digunakan dalam subrutin polyfit adalah 1. Berikut ini pemrograman MATLAB-nya.

Halaman 66 dari 101

Bab 5 Regresi Linier dan Non Linier %konduktivitas2 T = [4,5,6,8]*100;

% Absis

K = [240,237,232,220];

% Ordinat

[P,S] = polyfit(T,k,1)

>>konduktivitas2 P = -0.0511

261.6571

S = R: [2x2 double] df: 2 normr: 1.9712

5.2 Linierisasi

Seringkali ditemukan persamaan tak linier dalam permasalah real teknk kimia. Tentunya kita tak dapat begitu daja mengalurkan data-data dengan menggunakan pemodelan linier. Agar dapat dimodelkan dengan pemodelan linier, maka persamaan tak linier itu harus dilinierisasi terlebih dahulu. Berikut ini pemaparannya. LINIERISASI

y = ae

ln y = ln(aebx ) ln y = bx + ln a

bx

Tabel 5.1 Hasil linearisasi persamaan-persamaan tak linier Tipe persamaan

absis

ordinat

slope

intersep

y = ax + b

x

y

a

b

y = aebx

x

ln(y)

b

ln(a)

x

x/y

a

b

1/x

y

a

b

y=

x (ax + b)

y = a/ x+b

Halaman 67 dari 101


y = ax b + c

ln(x)

ln(y-c)

b

ln(a)

Kasus Suatu reaksi berorde n memiliki laju reaksi sbb:

r = kC An Apabila volume reaktor partaian (batch) konstan, persamaan laju reaksi menjadi sbb:

dC A = − kC An dt Tentukan orde laju reaksi tersebut jika diketahui data-data eksperimen sbb: Tabel 5.2 Konsentrasi, waktu, dan laju perubahannya Waktu (s)

CA (mol/liter) dCA/dt (mol/liter.s)

0

1.0

-0.10000

10

0.50

-0.02500

20

0.33

-0.01110

30

0.25

-0.00625

40

0.20

-0.00400

Jawaban: dC A = − kC An dt

⎛ dC ⎞ ln ⎜ − A ⎟ = n ln C A + ln k ⎝ dt ⎠

%linierisasi t = [0:10:40];

% waktu

CA = [1,0.50,0.33,0.25,0.20];

%konsentrasi A

dCAdt = -[0.1,0.025,0.0111,0.00625,0.00400]; %Laju y = log(-dCAdt);

Halaman 68 dari 101

Bab 5 Regresi Linier dan Non Linier x = log(CA);

[P S] = polyfit(x,y,1); n = P(1)

%orde reaksi

k = exp(P(2))

%konstanta laju reaksi

Hasil pada command window n = 1.9982 k = 0.1002

5.3 Regresi Linier Peubah Banyak

y ( x1 , x2 ) = a0 x1 + a1 x2 + a2 n

S = ∑ ( yi − a0 x1,i − a1 x2,i − a2 ) 2 i =1

n ∂S = −2∑ x1,i ( yi − a0 x1,i − a1 x2,i − a2 ) = 0 ∂a0 i =1 n

∑ (x i =1

1,i

n

∑x

yi − x1,2i a0 − x1,i x2,i a1 − x1,i a2 ) = 0

2

1,i

i =1

n

n

n

i =1

i =1

i =1

a0 + ∑ x1,i x2,i a1 + ∑ x1,i a2 = ∑ x1,i yi

n ∂S = −2∑ x2,i ( yi − a0 x1,i − a1 x2,i − a2 ) = 0 ∂a1 i =1 n

∑ (x i =1

2,i

n

∑x i =1

yi − x1,i x2,i a0 − x 2,2 i a1 − x2,i a2 ) = 0 n

n

n

i =1

i =1

x a + ∑ x 2 ,i a1 + ∑ x2,i a2 = ∑ x2,i yi

1,i 2,i 0

i =1

2

n ∂S = −2∑ ( yi − a0 x1,i − a1 x2,i − a2 ) = 0 ∂a2 i =1

n

∑x i =1

n

n

i =1

i =1

a + ∑ x2,i a1 + na2 = ∑ yi

1,i 0

Halaman 69 dari 101


Ketiga buah persamaan linier tersebut dapat dijelmakan dalam matrik sbb: ⎡ n 2 ⎢ ∑ x1,i ⎢ i =1 ⎢ n ⎢ ∑ x1,i x2,i ⎢ i =1 ⎢ n ⎢ ∑ x1,i ⎣ i =1

n

∑ x1,i x2,i i =1

n

∑x i =1

2 2,i

n

∑x i =1

2,i

⎤ ⎡ n ⎤ x ∑ 1,i ⎥ ⎢ ∑ x1,i yi ⎥ i =1 ⎥ ⎡ a0 ⎤ ⎢ i =1 ⎥ n ⎥⎢ ⎥ ⎢ n ⎥ x2,i ⎥ ⎢ a1 ⎥ = ⎢ ∑ x2,i yi ⎥ ∑ i =1 ⎥ ⎣⎢ a2 ⎦⎥ ⎢ i =1 ⎥ ⎥ ⎢ n ⎥ n ⎥ ⎢ ∑ yi ⎥ ⎦ ⎣ i =1 ⎦ n

Kasus Perhatikan data-data reaksi non-isotermal suatu reaksi reversibel berikut ini: r A ⎯⎯ →P

Tabel 5.3 Data laju reaksi CA(mol/liter)

Temperatur (K)

Laju reaksi (mol/liter.s)

1

373

1.508

0.023

395

2.936

1.15

365

1.293

0.87

400

3.242

1.05

405

4.566

0.75

388

1.899

0.55

410

2.780

0.65

380

1.255

Anggap reaksi tersebut memenuhi model persamaan laju sbb: ⎛ −E ⎞ r = k0C n exp ⎜ ⎟ ⎝ RT ⎠

Perkirakan harga k0, E dan n dari data-data yang tersedia. Jawaban. ⎛ −1 ⎞ ln r = n ln C + E ⎜ ⎟ + ln k0 ⎝ RT ⎠ yi = ln ri

x1,i = ln Ci

x2,i = −

1 RTi

Halaman 70 dari 101

Bab 5 Regresi Linier dan Non Linier % multivariabel clear clc CA = [1,0.023,1.15,0.87,1.05,0.75,0.55,0.65];

%mol/liter

T = [373,395,365,400,405,388,410,380];

%K

r = [1.508,2.936,1.293,3.242,4.566,1.899,2.780,1.255];%mol/liter.s

y = log(r); x1 = log(CA); x2 =-1./(0.082*T);

A = [sum(x1.^2),sum(x1.*x2), sum(x1) sum(x1.*x2),sum(x2.^2),sum(x2) sum(x1),sum(x2),length(r)]; c = [sum(x1.*y) sum(x2.*y) sum(y)]; a = A\c

Hasil pada Command window a = -0.0108 320.9052 10.8467

5.4 Regresi Non Linier

Pada bagian sebelumnya kita telah mempelajari regresi persamaan tak linier dengan terlebih dahulu melakukan linierisasi. Namun tidak semua persamaan tak linier dapat memberikan parameter yang akurat dengan linierisasi. Pada bagian ini kita akan mempelajari regresi persamaan tak linier sehingga kita tidak lagi harus melinierisasikan persamaan tak linier. Perhatikan fungsi tak linier sbb:

⎛ a1 ⎞ y = exp ⎜ a0 + ⎟ ( x + a2 ) ⎠ ⎝

Persamaan Antoine

a0, a1, dan a2 merupakan parameter.

Halaman 71 dari 101


⎛ ⎛ a1 ⎞ ⎞ S = ∑ ⎜⎜ yi − exp ⎜ a0 + ⎟⎟ ( xi + a2 ) ⎠ ⎟⎠ i =1 ⎝ ⎝ n

2

Turunan parsial terhadap a0 n ⎛ ⎛ ⎛ a1 ⎞ ⎞ a1 ⎞ ∂S = −2∑ ⎜⎜ yi − exp ⎜ a0 + ⎟ ⎟⎟ exp ⎜ a0 + ⎟=0 ( xi + a2 ) ⎠ ⎠ ( xi + a2 ) ⎠ ∂a0 i =1 ⎝ ⎝ ⎝

Turunan parsial terhadap a1 n ⎛⎛ ⎛ ⎛ a1 ⎞ ⎞ a1 ⎞ ⎛ 1 ⎞ ⎞ ∂S exp ⎜ a0 + = −2∑ ⎜ ⎜⎜ yi − exp ⎜ a0 + ⎟ ⎟ ⎟⎜ ⎟⎟ = 0 ⎜ ( xi + a2 ) ⎠ ⎠⎟ ( xi + a2 ) ⎠ ⎝ xi + a2 ⎠ ⎠⎟ ∂a1 i =1 ⎝ ⎝ ⎝ ⎝

Turunan parsial terhadap a2 n ⎛⎛ ⎞ ⎛ ⎛ a1 ⎞ ⎞ a1 ⎞ ∂S exp ⎜ a0 + a1 ln( xi + a2 ) ) ⎟ = 0 = −2∑ ⎜ ⎜⎜ yi − exp ⎜ a0 + ( ⎟ ⎟ ⎟ ⎜ ⎟ ( xi + a2 ) ⎠ ⎠⎟ ( xi + a2 ) ⎠ ∂a2 i =1 ⎝ ⎝ ⎝ ⎝ ⎠

Pada akhirnya diperoleh sistem persamaan tak linier yang terdiri atas 3 buah persamaan tak linier. Sistem persamaan tak linier dapat diselesaikan secara simultan menggunakan metode Newton seperti yang telah dibahas pada bab 3 persamaan tak linier.

5.5 Subrutin MATLAB: nlinfit

[beta,R] = nlinfit(x,y,modelfun,beta0)

Kasus Tabel 3.3 Tekanan uap dari Benzena (Perry) Temperatur, T o C -36.7 -19.6 -11.5 -2.6 7.6 15.4 26.1

Tekanan, P (mmHg) 1 5 10 20 40 60 100

Halaman 72 dari 101


42.2 60.6 80.1

200 400 760

Persamaan polynomial

P ( x) = a0 + a1 x + a2 x 2 + a3 x3 + ... + an x n Persamaan Clapeyron log( P ) = A +

B T

Persamaan Riedel

log( P) = A +

B + C log(T ) + DT β T

dengan harga β = 2 . a. Korelasikan data dengan berbagai orde persamaan polynomial dengan menganggap temperatur absolut (Kelvin) adalah variabel bebas dan P adalah variabel terikat. b. Korelasikan data dengan menggunakan persamaan Clapeyron c. Korelasikan data menggunakan persamaan Riedel d. Diskusikan persamaan manakah yang terbaik mewakili data-data eksperimental tersebut. Jawab: a. Polynom orde 2, 3, 4, 5 %polynom clear clc T =[-36.7,-19.6,-11.5,-2.6,7.6,15.4,26.1,42.2,60.6,80.1]+273;%K P = [1 5 10 20 40 60 100 200 400 760];%mmHg N =length(P); P2 = polyfit(T,P,2) Pmod2 = polyval(P2,T); R2 = Pmod2-P; Var2 = sum(R2.^2)/(N-2) P3 = polyfit(T,P,3) Pmod3 = polyval(P3,T); R3 = Pmod2-P; Var3 = sum(R3.^2)/(N-3) P4 = polyfit(T,P,4) Pmod4 = polyval(P4,T); R4 = Pmod2-P; Var4 = sum(R4.^2)/(N-4) P5 = polyfit(T,P,5) Pmod5 = polyval(P5,T); R5 = Pmod2-P; Var5 = sum(R5.^2)/(N-5)

Hasil di Command Window

Halaman 73 dari 101

Bab 5 Regresi Linier dan Non Linier P2 = 1.0e+003 * 0.0001

-0.0450

5.8560

Var2 = 1.0647e+003 Warning: Polynomial is badly conditioned. Remove repeated data points or try centering and scaling as described in HELP POLYFIT. > In polyfit at 81 In polinom at 9 P3 = 1.0e+004 * 0.0000

-0.0001

0.0146

-1.2519

Var3 = 1.2168e+003 Warning: Polynomial is badly conditioned. Remove repeated data points or try centering and scaling as described in HELP POLYFIT. > In polyfit at 81 In polinom at 11 P4 = 1.0e+004 * 0.0000

-0.0000

0.0001

-0.0248

1.5881

Var4 = 1.4196e+003 Warning: Polynomial is badly conditioned. Remove repeated data points or try centering and scaling as described in HELP POLYFIT. > In polyfit at 81 In polinom at 13

P5 = 1.0e+004 * -0.0000

0.0000

-0.0000

0.0002

-0.0339

2.1109

Var5 =

Halaman 74 dari 101

Bab 5 Regresi Linier dan Non Linier 1.7035e+003

b. Persamaan Clapeyron %clapeyron.m function logP = clapeyron(a,T); logP = a(1)+a(2)./T; %Persamaan Clapeyron

clear clc T =[-36.7,-19.6,-11.5,-2.6,7.6,15.4,26.1,42.2,60.6,80.1]; P = [1 5 10 20 40 60 100 200 400 760]; logP = log10(P); a0 = [0.1 0.3]; [a R]=nlinfit(T,logP,'clapeyron',a0) N = length(P); z = length(a0); Variance = sum(R.^2)/(N-z)

Hasil di Command Window a = 1.6667

1.9173

R = Columns 1 through 8 -1.6145 0.2598

-0.8699

-0.5000

0.3718

-0.3169

-0.0131

0.5889

Columns 9 through 10 0.9037

1.1902

Variance = 0.8124

c. Persamaan Riedel %riedel.m function logP = riedel(a,T) logP = a(1)+a(2)./T+a(3)*log10(T)+a(4)*T.^2;

%run_riedel.m clear clc T =[-36.7,-19.6,-11.5,-2.6,7.6,15.4,26.1,42.2,60.6,80.1]+273;%K P =[1 5 10 20 40 60 100 200 400 760];%mmHg

Halaman 75 dari 101

Bab 5 Regresi Linier dan Non Linier logP = log10(P);

a0 = [1 1 1 1];

[a R] = nlinfit(T,logP,'riedel',a0) N = length(P); z = length(a0); Variance = sum(R.^2)/(N-4)

Hasil di command window a = 1.0e+003 * 0.2162

-9.2955

-0.0756

0.0000

-0.0113

0.0080

-0.0046

0.0054

R = Columns 1 through 6 0.0107

-0.0236

0.0251

0.0081

Columns 7 through 10 -0.0065

-0.0112

Variance = 2.9689e-004

d. Pembahasan

Tugas

Nomor 1 Harga viskositas air (centipoise) telah diukur pada berbagai temperatur. Hasil dari eksperimen disajikan dibawah ini. Menggunakan regresi linier ganda (multiple regresi linier), carilah konstanta-konstanta yang sesuan dengan persamaan berikut: 1

μ

= k1 + k2T + k3T 2

T(oC)

10

20

30

40

50

60

70

μ (cp)

1.308

1.005

0.801

0.656

0.549

0.469

0.406

Nomor 2 Sebuah reaksi heterogen diketahui terjadi pada laju yang dapat digambarkan oleh model Langmuir-Hinshelwood berikut ini:

Halaman 76 dari 101


r=

k1 PA (1 + K A PA + K R PR ) 2

Dari pengukuran laju awal, k1 ditentukan sebagai 0.015 mol/s.g-cat.atm, pada 400 K. Dengan menggunakan data laju reaksi pada 400 K, perkirakan nilai dari KA dan KR. PA

1

0.9

0.8

0.7

0.6

0.5

0.4

PR

0

0.1

0.2

0.3

0.4

0.5

0.6

r

3.4x10-5 3.6x10-5 3.7x10-5 3.9x10-5 4.0x10-5 4.1x10-5 4.2x10-5

________________________________o0o______________________________

Halaman 77 dari 101

Bab 6 Integrasi

Bab 6 Integrasi (Under construction) 6.1 Metode Trapesium

f(x)

f(b)

f(a)

a

b

x

Gambar 5.1 Metode trapesium

f (a) + f (b) 2

b

∫ f ( x)dx = (b − a) a

6.2 Metode Simpson

x2

h

∫ f ( x)dx = 3 ( f ( x ) + 4 f ( x ) + f ( x ) ) 0

1

2

x0

6.3 Metode Romberg

6.4 Metode Gauss

Halaman 78 dari 101

Bab 6 Integrasi

6.5 Subrutin MATLAB: trapz

Subrutin trapz menghitung harga integral dari nilai-nilai diskrit x dan y dengan menggunakan metode Trapezoidal.

Z = trapz(x,y) Keterangan: x dan y adalah vektor

Kasus 1

Dua buah besaran yang sangat penting dalam pembelajaran proses-proses fermentasi adalan laju pembebasan CO2 dan laju pengambilan O2. Hal tersebut dihitung dari analisis eksperimental dari gas masuk dan gas keluar fermentor, dan laju alir, temperatur dan tekanan dari gas-gas ini. Rasio pembebasan CO2 terhadap pengambilan O2 menghasilkan RQ (Respiratory Quotient) yang merupakan barometer aktivitas metabolik dari mikroorganisme. Laju di atas dapat dintegrasikan untuk memperoleh jumlah keseluruhan dari CO2 yang diproduksi dan oksigen yang dikonsumsi selama fermentasi. Tabel berikut menunjukan laju respirasi dihtung dari fermentasi Penicillium chrysogenum yang menghasilkan antibiotik penicilin.

Tabel 5.1 Laju pembebasan CO2 dan laju pengambilan O2 Waktu fermentasi (jam) 140 141 142 143 144 145 146 147 148 149 150

Laju Pembebasan CO2 (g/jam) 15.72 15.53 15.19 16.56 16.21 17.39 17.36 17.42 17.60 17.75 18.95

Laju pengambilan O2 (g/jam) 15.49 16.16 15.35 15.13 14.20 14.23 14.29 12.74 14.74 13.68 14.51

Halaman 79 dari 101

Bab 6 Integrasi

Hitunglah jumlah keseluruhan CO2 yang dihasilkan dan Oksigen yang dikonsumsi selama fermentasi berlangsung.

%Respiratory_Quotient clear clc

t = [140:150];

% Waktu fermentasi

dCO2dt = [15.72,15.53,15.19,16.56,16.21,17.39,17.36,17.42,... 17.60,17.75,18.95];

% Laju pembebasan CO2

dO2dt = [15.49,16.16,15.35,15.13,14.20,14.23,14.29,12.74,... 14.74,13.68,14.51];

% Laju pengambilan O2

CO2 = trapz(t,dCO2dt) O2 = trapz(t,dO2dt)

RQ = CO2/O2

6.6 Subrutin MATLAB: quad

Q = quad(fungsi,A,B) Kasus 2 Harga kapasitas panas suatu material dapat dievaluasi dengan menggunakan persamaan sbb:

c p (T ) = 0.132 + 1.56x10−4 T + 2.64x10−7 T 2

cal/goC

Hitunglah besar entalpi material sebanyak 1 gram pada rentang temperatur -100 o

C s.d 200 oC dengan rumus sbb: T2

ΔH = m ∫ c(T )dT T1

Halaman 80 dari 101

Bab 6 Integrasi %kapasitas.m function cp = kapasitas(T)

cp=0.132+1.56e-4.*T+2.64e-7*T.^2;

%runkapasitas.m Q = quad('kapasitas',-100,200)

>> runkapasitas Q = 42.7320

6.7 Subrutin MATLAB: dblquad

Subrutin dblquad digunakan untuk menghitung integral lipat dua.

q = dblquad(fun,xmin,xmax,ymin,ymax) %integrnd.m function z = integrnd(x, y) z = y*sin(x)+x*cos(y);

%run_integrnd.m Q = dblquad(@integrnd,pi,2*pi,0,pi);

6.8 Subrutin MATLAB: triplequad

Subrutin triplequad digunakan untuk menghitung integral lipat tiga.

triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax)

%intgrnd3 function f = integrnd3(x,y,z) f = y*sin(x)+z*cos(x);

%run_integrnd3 Q = triplequad('integrnd3', 0, pi, 0, 1, -1, 1)

Halaman 81 dari 101

Bab 6 Integrasi

Tugas

Nomor 1 Lakukan komputasi yang sama seperti pada kasus 1, namun dengan massa material 2000 gram dan temperatur -200 s.d 100 oC.

Nomor 2 Profil kecepatan dari partikel pasir unggun terfluidakan dengan udara pada kecepatan 1 m/s diberikan pada tabel Tabel 1 dan Tabel 2. Hitunglah kecepatan gradien aksial ( ∂Vz / ∂z , ∂Vr / ∂z ). Plot rata-rata gradien z versus posisi radial dan bandingkanlah besar ordenya.

Posisi aksial, mm Posisi aksial, mm

25 75 125 175 225 275 325 375 425 475

4.7663 -13.09 -15.81 1.77 1.43 -5.07 13.11 11.7 8.18 3.35 -27.05

25 75 125 175 225 275 325 375 425 475

4.7663 93.33 244.73 304.34 308.81 379.66 416.08 184.46 55.74 -67.81 -136.25

14.2988 -37.66 -15.99 1.17 -0.57 -7.26 16.51 34.5 25.29 -0.39 -22.25

14.2988 74.12 217.07 260.58 281.67 328.52 366.96 157.25 -12.28 -118.77 -32.33

Posisi radial (mm) 33.3638 42.8962 -54.44 -58.21 -25.37 -22.3 5.5 1.63 2.44 0.2 -18.17 -17.3 21 20.29 71.44 73.49 37.07 30.05 -42.22 -57.42 -79.45 -110.08

52.4288 -41.35 -11.1 -1.79 -0.65 -10 15.64 64.88 2.61 -82.36 -116.62

61.9612 -23.97 -2.26 -0.26 0.35 -2.65 0.98 50.91 -17.06 -69.34 -128.25

Posisi radial (mm) 23.8313 33.3638 42.8962 69.35 43.68 18.8 177.09 103.79 16.87 201.15 118.82 22.76 209.18 133.9 53.88 279.3 165.61 53.25 314.09 203.08 44.97 111.99 63.23 1.03 -18.74 -47.26 -42.1 -108.46 -89.68 9.24 -65.5 -111.72 38.74

52.4288 -6.9 -39.74 -52.23 -51.92 -65.97 -76.93 -63.66 -9.95 61.78 115.6

61.9612 -21.56 -74.91 -82.86 -98.47 -133.92 -160.04 -71.23 125.57 175.43 84.88

23.8313 -52.41 -27.81 3.45 4.86 -18.43 19.32 58.3 31.18 -18 -49.45

Halaman 82 dari 101

71.4938 -7.21 1.63 1.09 2.21 0.29 -9.81 19.14 -15.88 -17.35 -76.49

Bab 7 Persamaan Diferensial Biasa

Bab 7 Persamaan Diferensial Biasa (PDB) Definisi PDB Persamaan diferensial biasa adalah persamaan diferensial yang terdiri atas fungsi turunan satu buah variabel bebas. Contoh: Persamaan gaya geser (shear stress) pada aliran fluida dirumuskan sbb.

dτ xz = ρg dx Perhatikan PDB hanya memiliki satu buah variabel bebas yaitu x dan satu variabel terikat yaitu τxz.

Aplikasi PDB PDB banyak ditemukan pada pemodelan-pemodelan teknik reaktor, kinetika reaksi kimia, peristiwa-peristiwa perpindahan dll.

Klasifikasi PDB Berdasarkan ordenya PDB terdiri atas tiga jenis (paling umum ditemukan dalam permasalahan teknik kimia). Orde 1

dy + y = kx dx

Orde 2

d2y dy +y = kx 2 dx dx

Orde 3

d3y d2y ⎛ dy ⎞ a + + b ⎜ ⎟ = kx 3 2 dx dx ⎝ dx ⎠

2

Berdasarkan ordenya PDB terdiri atas dua jenis. 1. Linier Persamaan umum PDB linier dirumuskan sbb:

Halaman 83 dari 101


bo ( x )

dny d n −1 y dy b x + + ... + bn −1 ( x ) + bn ( x ) y = R ( x ) ( ) 1 n n −1 dx dx dx

2. Taklinier PDB yang tidak memenuhi persamaan umum PDB linier di muka dikelompokan ke dalam PDB tak linier. Berdasarkan kondisi batasnya PDB terdiri atas dua jenis. 1. PDB bernilai awal ∂2 y = − yx ∂x 2 y (0) = 2,

∂y (0) = 1 ∂x

harga x sama

2. PDB bernilai batas

∂2 y = − yx ∂x 2 y (0) = 2, y (1) = 1 harga x berbeda

Transformasi ke Dalam Bentuk Kanonikal

Persamaan diferensial biasa linier orde 1 bernilai awal dapat diselesaikan dengan menggunakan metode matrik eksponensial dan metode eigen yang akan dibahas di depan nanti. PDB linier orde 2, 3 bernilai awal dapat pula diselesaikan dengan metode-metode tersebut, asalkan PDB tersebut ditransformsikan terlebih dahulu ke dalam PDB orde 1. Berikut ini penjelasan teknik transformasi dari PDB berorde tinggi menjadi PDB berorde 1. Contoh 1: d4z d 3z d 2z dz + 5 − 2 − 6 + 3z = 0 4 3 2 dt dt dt dt

Halaman 84 dari 101


Transformasi PDB orde 4 linier tersebut akan menghasilkan 4 buah PDB linier orde 1. Misalkan: z = y1 dz dy1 = = y2 dt dt d 2 z dy2 = = y3 dt 2 dt d 3 z dy3 = = y4 dt 3 dt d 4 z dy4 = dt 4 dt

Maka PDB orde 4 dapat dituliskan sbb:

dy1 dt dy2 dt dy3 dt dy4 dt

= y2 = y3

Penulisan dalam bentuk matrik sbb:

= y4

⎡ dydt1 ⎤ ⎡ 0 ⎢ dy2 ⎥ ⎢ ⎢ dt ⎥ = ⎢ 0 ⎢ dy3 ⎥ ⎢ 0 ⎢ dydt ⎥ ⎢ ⎢⎣ dt4 ⎥⎦ ⎣ −3

1 0 0 6

= −3 y1 + 6 y2 + 2 y3 − 5 y4

Contoh 2: d4z d 3z d 2z dz + 5 − 2 − 6 + 3z = e−t 4 3 2 dt dt dt dt Transformasi PDB orde 4 linier tersebut akan menghasilkan 5 buah PDB linier orde 1. Misalkan: z = y1 dz dy1 = = y2 dt dt d 2 z dy2 = = y3 dt 2 dt d 3 z dy3 = = y4 dt 3 dt d 4 z dy4 = dt 4 dt −t y5 = e dy5 = −e − t = − y5 dt

Maka PDB orde 4 dapat dituliskan sbb:

dy1 dt dy2 dt dy3 dt dy4 dt dy5 dt

= y2 = y3 = y4 = −3 y1 + 6 y2 + 2 y3 − 5 y4 + y5 = − y5

Penulisan dalam bentuk matrik sbb:

Halaman 85 dari 101

0 0 ⎤ ⎡ y1 ⎤ 1 0 ⎥⎥ ⎢⎢ y2 ⎥⎥ 0 1 ⎥ ⎢ y3 ⎥ ⎥⎢ ⎥ 2 −5⎦ ⎢⎣ y4 ⎥⎦


⎡ dydt1 ⎤ ⎡ 0 ⎢ dy2 ⎥ ⎢ ⎢ dt ⎥ ⎢ 0 ⎢ dy3 ⎥ = ⎢ 0 ⎢ dt ⎥ ⎢ ⎢ dydt4 ⎥ ⎢ −3 ⎢ dy ⎥ ⎢ ⎢⎣ dt5 ⎥⎦ ⎣ 0

1 0 0 6 0

0 0 0 ⎤ ⎡ y1 ⎤ 1 0 0 ⎥⎥ ⎢⎢ y2 ⎥⎥ 0 1 0 ⎥ ⎢ y3 ⎥ ⎥⎢ ⎥ 2 −5 1 ⎥ ⎢ y4 ⎥ 0 0 −1⎥⎦ ⎢⎣ y5 ⎥⎦

Contoh 3: 3

2 d 3z ⎛ dz ⎞ 2 d z z + − ⎜ ⎟ − 2z = 0 3 2 dx dx ⎝ dx ⎠

Transformasi PDB orde 3 taklinier. Misalkan: z = y1 dz dy1 = = y2 dx dx d 2 z dy2 = = y3 dx 2 dx d 3 z dy3 = dx 3 dx

Maka PDB orde 3 taklinier dituliskan sbb:

dy1 = y2 dx dy2 = y3 dx dy3 = 2 y1 − y12 y3 + y23 dx

PDB taklinier tidak dapat dituliskan dalam bentuk matrik.

Halaman 86 dari 101


Contoh 4:

d 3 z 3 d 2 z 2 dz +t −t + 5z = 0 dt 3 dt 2 dt Transformasi PDB orde 3 taklinier. Misalkan: z = y1 dz dy1 = = y2 dt dt d 2 z dy2 = = y3 dt 2 dt d 3 z dy3 = dt 3 dt y4 = t

Maka PDB orde 3 taklinier dituliskan sbb:

dy4 =1 dt

dy1 dt dy2 dt dy3 dt dy4 dt

= y2 = y3 = −5 y1 + y42 y2 − y43 y3 =1

PDB taklinier tidak dapat dituliskan dalam bentuk matrik.

Nilai dan Vektor Eigen Aw[ k ] = λk w[ k ] Aw[ k ] − λk w[ k ] = 0

( A − λk I ) w[ k ] = 0 ( A − λk I ) = 0

atau w[ k ] ≠ 0

det ( A − λk I ) = 0

Vektor eigen tidak bernilai nol

Keterangan: A adalah sebuah matrik kubus

λk adalah nilai eigen w[ k ] adalah vektor eigen

Halaman 87 dari 101


Berikut ini akan dipaparkan cara menghitung nilai dan vektor eigen secara analitik.

Kasus 6

Tentukanlah vektor dan nilai eigen dari matrik A berikut ini dengan menggunakan cara analitik.

⎡1 0 3 ⎤ A = ⎢⎢0 2 1 ⎥⎥ ⎣⎢ 3 1 −1⎦⎥ Jawaban:

( A − λk I ) w[ k ] = 0 det ( A − λk I ) = 0 (1 − λ )

0

3

0

(2 − λ )

1

3

−1

(−1 − λ )

(1 − λ )

0

3

(1 − λ )

0

(2 − λ )

1

0

3

−1

(−1 − λ )

3

=0

0 (2 − λ ) = 0 −1

(1 − λ )(2 − λ )(−1 − λ ) − (3)(2 − λ )(3) − (1 − λ ) = 0 (1 − λ 2 )(λ − 2) − 9(2 − λ ) − (1 − λ ) = 0

λ − 2 − λ 3 + 2λ 2 − 18 + 9λ − 1 + λ = 0 −λ 3 + 2λ 2 + 11λ − 21 = 0 Dengan menggunakan subrutin roots MATLAB diperoleh harga akar-akar polinom pangkat 3 (nilai eigen) tersebut, yaitu:

λ1 = 3.4211 λ2 = −3.2880 λ3 = 1.8669 Kembali ke persamaan awal.

( A − λk I ) w[ k ] = 0 0 3 ⎤ ⎡ w1 ⎤ ⎡ (1 − λ ) ⎢ 0 (2 − λ ) 1 ⎥⎥ ⎢⎢ w2 ⎥⎥ = 0 ⎢ ⎢⎣ 3 1 (−1 − λ ) ⎥⎦ ⎢⎣ w3 ⎥⎦

Halaman 88 dari 101


Karena vektor eigen (w) tidak bernilai nol, maka kita misalkan harga w3 sebagai basis bernilai 1.

(1 − λ ) w1 + 3w3 = 0 (2 − λ ) w2 + w3 = 0 Misalkan w3 = 1, maka system persamaan linier menjadi −3 (1 − λ ) −1 w2 = (2 − λ ) w3 = 1 w1 =

(1 − λ ) w1 = −3 (2 − λ ) w2 = −1 Masukan harga nilai eigen Untuk:

λ1 = 3.4211

λ2 = −3.2880

⎡1.2391 ⎤ w = ⎢⎢0.7037 ⎥⎥ ⎢⎣ 1 ⎥⎦ [1]

[2]

w

⎡ −0.6996 ⎤ = ⎢⎢ −0.1891⎥⎥ ⎢⎣ 1 ⎥⎦

λ3 = 1.8669 [3]

w

⎡ 3.4607 ⎤ = ⎢⎢ −7.5131⎥⎥ ⎢⎣ 1 ⎥⎦

Normalisasi vektor-vektor eigen tersebut dengan menggunakan norma ke-2. w[1]

2

=

(1.2391

2

+ 0.7037 2 + 12 ) = 1.741

⎡ 1.2391 ⎤ 1.741 ⎥ ⎡0.7117 ⎤ ⎢ ⎥ = ⎢ 0.4042 ⎥ w[1] = ⎢ 0.7037 ⎥ 1.741⎥ ⎢ ⎢ ⎢ ⎥⎦ 0.5744 ⎢ 1 ⎥ ⎣ ⎢⎣ ⎥ 1.741 ⎦

w[2]

w[2]

2

=

( 0.6996

2

+ 0.18912 + 12 ) = 1.235

⎡ −0.6996 ⎤ 1.235⎥ ⎡ −0.5665⎤ ⎢ ⎥ = ⎢ −0.1531⎥ = ⎢ −0.1891 ⎥ 1.235 ⎥ ⎢ ⎢ ⎢ ⎥ ⎢⎣ 0.8097 ⎥⎦ 1 ⎢⎣ 1.235 ⎥⎦

Halaman 89 dari 101


w[3]

2

=

( 3.4607

2

+ 7.51312 + 12 ) = 8.332

⎡ 3.4607 ⎤ 8.332 ⎥ ⎡ 0.4153 ⎤ ⎢ ⎥ = ⎢ −0.9017 ⎥ w[3] = ⎢ −7.5131 ⎥ 8.332 ⎥ ⎢ ⎢ ⎢ ⎥⎦ 0.1200 ⎢ ⎥ ⎣ 1 8.332 ⎦⎥ ⎣⎢ Jadi nilai dan vektor eigen matrik A adalah ⎡ 3.4211 ⎤ λ = ⎢⎢ −3.2880 ⎥⎥ ⎢⎣ 1.8669 ⎥⎦

⎡0.7117 −0.5665 0.4153 ⎤ w = ⎢⎢ 0.4042 −0.1531 −0.9017 ⎥⎥ ⎢⎣ 0.5744 0.8097 0.1200 ⎥⎦

Catatan: Perkalian konstanta dengan vektor eigen tidak akan mengubah esensi dari vektor eigen tersebut. Untuk persoalan ini harga vektor eigen yang diperoleh menggunakan MATLAB (sekejap lagi akan dibahas) adalah hasil perkalian antara -1 dengan vektor eigen yang telah diperoleh pada perhitungan secara analitik.

MATLAB telah menyediakan rutin untuk menghitung nilai dan vektor eigen matriks A yaitu eig. Penulisan perintahnya pada MATLAB command window sbb:

[V , D ] = eig ( A) Vektor eigen

Nilai eigen

Sebagai contoh berikut ini akan ditampilkan perintah pada command window untuk menghitung nilai dan vektor eigen dari matrik A yang telah diselesaikan secara analitik sebelumnya. >> [V,D]=eig(A) V= -0.5665 -0.4153 -0.7118 -0.1531 0.9018 -0.4042 0.8097 -0.1200 -0.5744

Halaman 90 dari 101


D= -3.2880

0

0

0 1.8669

0

0

0 3.4211

Tugas 6 Transformasi kanonikal PDB dan analisis eigen

Nomor 1 Hitunglah nilai dan vektor eigen dari matrik A berikut ini. Bandingkan hasilnya dengan menggunakan subrutin eig di MATLAB. ⎡1 2 3⎤ A = ⎢⎢ 2 5 1 ⎥⎥ ⎢⎣ 3 1 4 ⎥⎦

Nomor 2 Ubahlah persamaan differensial berikut ke dalam bentuk kanonikal. a.

d 2x dx − 3 − 10 x = 0 2 dt dt

b.

d 3T 3 d 2T 2 dT +t −t − 10T = 0 dt 3 dt 2 dt

c.

2 d3y ⎛ dy ⎞ 3 d y y − −⎜ ⎟ + 9y = 0 3 2 dx dx ⎝ dx ⎠

2

Halaman 91 dari 101


Solusi Persamaan Differensial Biasa Linier bernilai awal

1. Metode matriks eksponensial

dy = Ay dt

y ( 0 ) = y0

A adalah matriks persegi (m x m) dan y adalah vektor kolom (m x 1) Integrasikan persamaan diferensial linier tersebut. y

t

dy ∫y y = A∫0 dt 0

ln

y = At y0

y = e At y0 Fungsi matriks eksponensial dapat dituliskan sebagai berikut: exp(At ) = I + At +

A 2 t 2 A 3t 3 + + ... 2! 3!

Contoh soal: Kasus 7

Berikut ini adalah PDB linier orde 2. d 2x dx − 3 − 10 x = 0 2 dt dt Dengan nilai awal pada t = 0, sbb: x(0) = 3 dx 0 = 15 dt Selesaikan PDB tercetak menggunakan metode matrik eksponensial dalam interval 0 ≤ t ≤ 1.0 (Langkah integrasi 0.1).

Halaman 92 dari 101


Jawaban : d 2x dx − 3 − 10 x = 0 2 dt dt x = y1

⎡ dy1 ⎤ ⎢ dt ⎥ ⎡ 0 1⎤ ⎡ y1 ⎤ ⎢ ⎥=⎢ ⎥⎢ ⎥ ⎢ dy2 ⎥ ⎣1 3⎦ ⎣ y2 ⎦ ⎣⎢ dt ⎦⎥

dx dy1 = = y2 dt dt d 2 x dy2 = = 3 y2 + y1 dt 2 dt

A

dy dt

y

dy = Ay dt Integrasikan.

⎡3⎤ y0 = ⎢ ⎥ ⎣15⎦

y = e At y0

Rentang integrasi 0 ≤ t ≤ 1.0. Langkah integrasi 0.1 t

0 0.1

⎡0 ⎡ y1 ⎤ ⎛ ⎢⎣ t ⎢ y ⎥ = ⎜⎜ e ⎣ 2⎦ ⎝

t⎤ ⎥ 3t ⎦

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9 1.0

⎞⎡ 3 ⎤ ⎟ ⎟ ⎢⎣15⎥⎦ ⎠

Dengan mensubstitusikan t = 0 s.d 1 (langkah integrasi 0.1) selesailah persoalan ini. Berikut ini pemrograman MATLAB-nya. kasus7.m

Halaman 93 dari 101


clear clc A = [0 1 1 3]; % Nilai awal yo = [3;15]; t = [0:0.1:1]'; for i = 1:length(t) y(i,:) = expm(A*t(i))*yo; end %kurva t-x x = y(:,1) plot(t,x) xlabel('t') ylabel('x') grid on

>>kasus7a t= 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000

x= 3.0000 4.7688 7.2122 10.5945

Halaman 94 dari 101

Bab 7 Persamaan Diferensial Biasa 15.2839 21.7922 30.8319 43.3941 60.8578 85.1416 118.9150

120

100

x

80

60

40

20

0

0

0.1

0.2

0.3

0.4

0.5 t

0.6

0.7

0.8

0.9

1

2. Metode nilai-vektor eigen Harga e At dapat dihitung dengan menggunakan bantuan nilai dan vektor eigen.

e At = Ve Dt V −1 Sehingga solusi PDB linier menjadi.

y = ⎡⎣ Ve Dt V −1 ⎤⎦ y0 Untuk lebih memahami metode nilai-vektor eigen berikut ini disajikan penyelesaian kasus 7 dengan menggunakan metode nilai-vektor eigen. Langkah awal sama dengan metode matriks eksponensial. ⎡ 0 1⎤ A=⎢ ⎥ ⎣ 1 3⎦

Halaman 95 dari 101


Dengan menggunakan rutin eig MATLAB diperoleh harga nilai (D) dan vektor eigen (V). >> A=[0 1 1 3] A= 0 1 1 3 >> [V,D]=eig(A) V= -0.9571 0.2898 0.2898 0.9571 D= -0.3028 0 0 3.3028

⎡ −0.9571 0.2898⎤ V =⎢ ⎥ ⎣ 0.2898 0.9571⎦

0 ⎤ ⎡ −0.3028 dan D = ⎢ 3.3028⎥⎦ ⎣ 0

Substitusikan matriks V dan D ke dalam persamaan

y = ⎡⎣ Ve Dt V −1 ⎤⎦ y0 ⎡ ⎡ −0.9571 0.2898⎤ ⎡e−0.3028t y = ⎢⎢ ⎥⎢ ⎢⎣ ⎣ 0.2898 0.9571⎦ ⎣ 0

−1 0 ⎤ ⎡ −0.9571 0.2898⎤ ⎤ ⎡ 3 ⎤ ⎥ ⎢ ⎥ ⎥⎢ ⎥ e0.3028t ⎦ ⎣ 0.2898 0.9571⎦ ⎥⎦ ⎣15⎦

Dengan mensubstitusikan t = 0 s.d 1 (langkah integrasi 0.1) selesailah persoalan ini. Berikut ini pemrograman MATLAB-nya.

Halaman 96 dari 101


kasus7.m clear clc A = [0 1 1 3]; % Nilai awal yo = [3;15]; a=length(yo); % Vektor dan Nilai eigen [V,D]=eig(A); % Rentang integrasi t=[0:0.1:1]' x =zeros(length(t),a); for i = 1 : length(t) y = (V*diag(exp(diag(D)*t(i)))*inv(V))*yo; x(i,:) = y; end x % kurva t-x plot(t,x(:,1)) xlabel('t') ylabel('x') grid on

Eksekusi program kasus7.m (lanjutan) Hasil di Command Window : >>kasus7 t= 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000

x= 3.0000 15.0000 4.7688 20.6902 7.2122 28.6127 10.5945 39.6409 15.2839 54.9901 21.7922 76.3512 30.8319 106.0768 43.3941 147.4403 60.8578 204.9959 85.1416 285.0806 118.9150 396.5110

Halaman 97 dari 101


120

100

x

80

60

40

20

0

0

0.1

0.2

0.3

0.4

0.5 t

0.6

0.7

0.8

0.9

1

Tugas 7 Metode eigen untuk menyelesaikan sistem persamaan diferensial biasa linier

Suatu bahan radioaktif meluruh berdasarkan mekanisme reaksi berantai sbb: k1 k2 A ⎯⎯ → B ⎯⎯ →C

k1 dan k2 adalah konstanta laju reaksi. B adalah produk intermediate dan C adalah produk akhir. Persamaan laju reaksinya sbb:

dC A = − k1C A CA, CB, dan CC adalah konsentrasi dt bahan A, B, dan C. dCB = k1C A − k2CB -1 dt k1= 3 s , k2= 1 s-1. dCC = k2CB Konsentrasi mula-mula bahan sbb: dt CA(0)=1 mol/m3 CB(0)=0 CC(0)=0

Halaman 98 dari 101


a) Menggunakan metode matriks eksponensial dan metode eigen tentukan konsentrasi CA, CB, dan CC sebagai fungsi waktu. b) Hitunglah konsentrasi CA, CB, dan CC saat t = 1 s dan t = 10 s. c) Buatlah profil konsentrasi A, B, dan C.

Sistem PDB tak linier bernilai awal

Sistem PDB tak linier bernilai awal banyak ditemukan pada kajian-kajian teknik kimia berikut ini: 1. Neraca massa dan energi yang melibatkan reaksi kimia. 2. Sistem dinamik nyata (tak ideal). 3. Hampir seluruh peristiwa perpindahan dalam teknologi proses.

Solusi PDB tak linier bernilai awal

1. Metode Euler

y −y y −y dy = f ( x, y ) = i +1 i = i +1 i dx xi +1 − xi h yi +1 = yi + hf ( xi , yi ) yi +1 = yi + hf ( xi , yi ) + O ( h 2 ) Kesalahan pemotongan lokal

Harga baru = harga lama + ukuran langkah x slope 2. Metode Runge-Kutta orde ke-2 (Crank-Nicholson)

yi +1 = yi + 12 ( k1 + k2 ) + O ( h3 ) k1 = hf ( xi , yi ) k2 = hf ( xi + h, yi + k1 )

3. Metode Runge-Kutta orde ke-3

Halaman 99 dari 101


yi +1 = yi + 16 ( k1 + 4k2 + k3 ) + O ( h 4 ) k1 = hf ( xi , yi ) h k ⎞ ⎛ k2 = hf ⎜ xi + , yi + 1 ⎟ 2 2⎠ ⎝ k3 = hf ( xi + h, yi + 2k2 − k1 ) 3. Metode Runge-Kutta orde ke-4 yi +1 = yi + 16 ( k1 + 2k2 + 2k3 + k4 ) + O ( h5 ) k1 = hf ( xi , yi ) k ⎞ h ⎛ k2 = hf ⎜ xi + , yi + 1 ⎟ 2 2⎠ ⎝ h k ⎞ ⎛ k3 = hf ⎜ xi + , yi + 2 ⎟ 2 2⎠ ⎝ k4 = hf ( xi + h, yi + k3 )

4. Metode Runge-Kutta orde ke-5 yi +1 = yi + 901 ( 7 k1 + 32k3 + 12k4 + 32k5 + 7 k6 ) + O ( h6 ) k1 = hf ( xi , yi ) h k ⎞ ⎛ k2 = hf ⎜ xi + , yi + 1 ⎟ 2 2⎠ ⎝ h 3k k ⎞ ⎛ k3 = hf ⎜ xi + , yi + 1 + 2 ⎟ 4 16 16 ⎠ ⎝ k ⎞ h ⎛ k4 = hf ⎜ xi + , yi + 3 ⎟ 2 2⎠ ⎝ 3h 3k 6k 9k ⎞ ⎛ k5 = hf ⎜ xi + , yi − 2 + 3 + 4 ⎟ 4 16 16 16 ⎠ ⎝ k1 4k 6k 12k4 8k5 ⎞ ⎛ k6 = hf ⎜ xi + h, yi + + 2 + 3 − + ⎟ 7 7 7 7 7 ⎠ ⎝

Halaman 100 dari 101


5. Metode Runge-Kutta-Fehlberg 16 6656 28561 yi +1 = yi + h ( 135 k1 + 12825 k3 + 56430 k4 − 509 k5 + 552 k6 )

k1 = f ( xi , yi ) k2 = f ( xi + h4 , yi + h4 k1 ) k3 = f ( xi + 38h , yi + 38h k1 + 932h k2 ) 7200 h 7296 h h k4 = f ( xi + 1213h , yi + 1932 2197 k1 − 2197 k 2 + 2197 k3 ) 3680 h 845 h h k5 = f ( xi + h, yi + 439 216 k1 − 8hk 2 + 513 k3 − 4104 k 4 ) 1859 h 11h h k6 = f ( xi + h2 , yi − 827h k1 + 2hk2 − 3544 2565 k3 + 4104 k 4 − 40 k5 ) 128 2197 1 d i +1 = h ( 360 k1 − 4275 k3 − 75240 k4 + 501 k5 + 552 k6 )

Algoritma dan pemrograman metode euler

Algoritma 1. Mulai 2. Definisikan fungsi f(x,y) [turunan ke-1] 3. Tetapkan nilai awal, x awal, langkah integrasi (h), dan x akhir. 4. Hitung panjang vektor x 5. Hitung yi +1 = yi − hf ( xi , yi ) untuk i = 1 s.d (panjang vektor-1) 6. selesai



Algorima dalam diagram alir

Mulai

Definisikan fungsi f(x,y)

Masukan harga yo, xo, h, xn

Hitung panjang vector x

yi +1 = yi − hf ( xi , yi ) untuk i = 1 s.d (panjang vektor x -1)

Selesai

eulerku.m function [x y] = eulerku(dfungsi,yo,xo,dx,xn) x = xo:dx:xn a = length(x) y(1) = yo for i = 1 : (a-1) y(i+1)= y(i) + dx*feval(dfungsi,x(i),y(i)) end



dcoba.m function f = dcoba(x,y) f = 2*x

>> [x y] = eulerku('dcoba',0,0,1,10) x= Columns 1 through 9 0

1

2

3

4

5

6

7

8

Columns 10 through 11 9 10

y= Columns 1 through 9 0

0

2

6 12 20 30 42 56

Columns 10 through 11 72 90

Subrutin dalam MATLAB untuk solusi PDB bernilai awal

Pada bagian ini akan dijelaskan subrutin ode23 dalam MATLAB untuk menyelesaikan PDB bernilai awal dengan karakter linier, taklinier, tunggal maupun jamak (sistem).



Cara penulisan sintaks ode23 [t,y] = ode23(‘fungsiPDB’,rentang_t,y0)

Fungsi PDB yang akan dievaluasi

Rentang integrasi

Harga awal

Misal: dy = −2 x 3 + 12 x 2 − 20 x + 8.5 dengan kondisi awal y = 1 pada x = 0 dan rentang dx integrasi dari x = 0 s.d x = 4. Berikut ini pemrograman MATLAB-nya. %pdb.m function dydx = pdb(x,y) dydx = -2*x^3+12*x^2-20*x+8.5; %runpdb.m clear clc rentang_x = [0 4]; y0 = 1; [x,y] = ode23('pdb',rentang_x,y0) plot(x,y) xlabel('x') ylabel('y')



Eksekusi di command window x=

y= 1.0000 1.0791 1.4488 2.1817 2.7701 3.1483 3.3031 3.2609 3.0709 2.7893 2.4787 2.2006 2.0131 1.9808 2.2494 2.6978 3.2901 3.8780 4.3716 4.6749 4.6862 4.3176 3.4933 3.0015

0 0.0094 0.0565 0.1712 0.3046 0.4524 0.6111 0.7777 0.9511 1.1319 1.3196 1.5150 1.7217 1.9512 2.2503 2.4902 2.7301 2.9516 3.1622 3.3655 3.5614 3.7483 3.9277 4.0000

5 4.5 4

y

3.5 3 2.5 2 1.5 1

0

0.5

1

1.5

2 x

2.5

3

3.5

4



Kasus 8 Studi terhadap kinetika proses fermentasi berhasil dimodelkan secara matematis sbb: ⎛ dy1 y ⎞ = k1 y1 ⎜1 − 1 ⎟ dt ⎝ k2 ⎠ dy2 = k3 y1 − k4 y2 dt

Dengan k1 = 0.03120; k2 = 47.70; k3 = 3.374 ;k4 = 0.01268 serta nilai pada t = 0, y1=5, y2=0. Evaluasi harga y1 dan y2 dalam interval waktu 0 s.d 10 jam setiap jamnya!. Berikut ini pemrograman MATLAB-nya. %fermen.m function dydt = fermen(t,y) k1=0.03120; k2=47.70; k3=3.374; k4=0.01268; dydt=[k1*y(1)*(1-y(1)/k2) k3*y(1)-k4*y(2)];

%kasus8 clear clc tspan = [0:1:10]; yo = [5 0]; [t,y]=ode23('fermen',tspan,yo)



Eksekusi di command window t= 0 1 2 3 4 5 6 7 8 9 10

y= 5.0000 5.1414 5.2863 5.4347 5.5868 5.7425 5.9020 6.0652 6.2323 6.4033 6.5783

0 17.0000 34.2657 51.8056 69.6282 87.7422 106.1564 124.8796 143.9206 163.2886 182.9924

Tugas 8 Solusi PDB tak linier menggunakan subrutin MATLAB ode23

Nomor 1 Konversi glukosa menjadi asam glukonik merupakan reaksi oksidasi sederhana dari gugus aldehid gula menjadi gugus karboksil. Enzim glukosa oksidase, terbentuk

dalam

mikroorganisme

untuk

mengubah

glukosa

menjadi

glukonolaktona. Kemudian glukonolaktona bereaksi dengan air membentuk asam glukonik. Mekanisme reaksi secara keseluruhan proses fermentasi dapat dituliskan sbb: Pertumbuhan sel: Glukosa + sel → Oksidasi glukosa:

sel

Glukosa oksidase Glukosa + O 2 ⎯⎯⎯⎯⎯→ Glukonolactona + H 2 O 2

Hidrolisis glukonolactona: Glukonolactona + H 2 O → Asam glukonik Dekomposisi peroksida: 1 Katalis H 2 O 2 ⎯⎯⎯ → H 2O + O2 2

Model matematika untuk fermentasi bakteri Pseudomonas ovalis yang memproduksi asam glukonik telah dirumuskan oleh Rai dan Constantinides. Model yang menggambarkan dinamika pertumbuhan fasa logaritmik ini dapat dituliskan sbb:



Laju pertumbuhan sel: ⎛ y ⎞ dy1 = b1 y1 ⎜1 − 1 ⎟ dt ⎝ b2 ⎠ Laju pembentukan glukonolaktona: dy2 b3 y1 y4 = − 0.9082b5 y2 dt b4 + y4 Laju pembentukan asam glukonik: dy3 = b5 y2 dt Laju konsumsi glukosa: ⎛byy ⎞ dy4 = −1.011⎜ 3 1 4 ⎟ dt ⎝ b4 + y4 ⎠

Keterangan: y1 = konsentrasi sel y 2 = konsentrasi glukonolaktona y3 = konsentrasi asam glukonik y 4 = konsentrasi glukosa b1 s.d b5 = parameter, f(T,pH) Pada kondisi operasi T=30 oC dan pH 6.6 harga dari lima parameter yang diperoleh secara eksperimental sbb: b1 = 0.949 b2 = 3.439 b3 = 18.72 b4 = 37.51 b5 = 1.169 Pada kondisi tersebut buatlah profil y1,y2,y3, dan y4 terhadap waktu selama 0 ≤ t ≤ 9 jam!. Nilai-nilai awal (pada saat t=0) adalah sbb: y1(0) = 0.5 U.O.D./mL y2(0) = 0.0 mg/mL y3(0) = 0.0 mg/mL y4(0) = 50.0 mg/mL Petunjuk: soal ini mudah….!!!

_______________________________o0o________________________________


Bab 8 Persamaan Diferensial Parsial

Bab 8 Persamaan Diferensial Parsial (PDP) Definisi PDP Persamaan diferensial parsial adalah persamaan diferensial yang terdiri atas fungsi turunan lebih dari satu variabel bebas. Contoh: Persamaan konduksi tak tunak satu dimensi pada lembaran suatu material dirumuskan dalam PDP sbb.

∂T k ∂ 2T = ∂t ρc p ∂x 2 PDP tersebut terdiri dari dua buah variabel bebas, yaitu x (tebal lembaran) dan t (waktu). Sedangkan variabel terikatnya adalah T (Temperatur). Jika digambarkan pada koordinat kartesius akan diperoleh gambar 3 dimensi.

Aplikasi PDP Pemodelan kimia dan fisika atas suatu fenomena dalam bidang teknik kimia seringkali menghasilkan rumusan matematis dalam bentuk PDP khususnya pada berbagai fenomena perpindahan (momentum, massa, dan panas). Klasifikasi PDP Berdasarkan ordenya PDP terdiri atas tiga jenis. Orde 1 ∂u ∂u −α =0 ∂x ∂y Orde 2

∂ 2u ∂u +u =0 2 ∂x ∂y Orde 3


Bab 8 Persamaan Diferensial Parsial 2

⎛ ∂ 3u ⎞ ∂ 2u ∂u + + =0 ⎜ 3⎟ ⎝ ∂x ⎠ ∂x∂y ∂y

Umumnya PDP dalam teknik kimia berorde 2 dengan 2 sampai 4 variabel bebas. Berdasarkan kelinierannya PDP terdiri atas tiga jenis. 1. Linier 2. Quasilinier 3. Taklinier Kajian PDP dalam diktat ini terbatas hanya untuk PDP orde 2 linier saja. PDP linier orde 2 memiliki persamaan umum sbb: a

∂ 2u ∂ 2u ∂ 2u ∂u ∂u 2 b c + + +d + e + fu + g = 0 2 2 ∂x ∂x∂y ∂y ∂x ∂y

Berdasarkan harga b 2 − ac PDP orde 2 linier terbagi atas tiga bagian, sbb: 1. Eliptik

b 2 − ac < 0

2. Parabolik

b 2 − ac = 0

3. hiperbolik

b 2 − ac > 0

Kondisi Batas

Untuk lebih memahami kondisi batas pada PDP perhatikan contoh persamaan konduksi tak tunak satu dimensi sbb:

∂T k ∂ 2T = ∂t ρc p ∂x 2 1. Kondisi Dirichlet Nilai variabel terikat (T) diketahui pada nilai variabel bebas (x,t) tertentu



kondisi awal T = f ( x)

pada t = 0 dan 0 ≤ x ≤ 1, atau

T = T0

pada t = 0 dan 0 ≤ x ≤ 1

kondisi batas T = f (t )

pada x = 0 dan t > 0, dan

T = T1

pada x = 0 dan t > 0

Contoh :

∂ 2T ∂T = ∂x 2 ∂t T ( 0, t ) = 1200 C

α

T (1, t ) = 1200 C T ( x = 0 :1, 0 ) = 250 C

2. Kondisi Neuman Turunan variabel terikat (T) diketahui pada nilai tertentu atau sebagai fungsi dari

variabel bebas (x,t). ∂T =0 ∂x

pada x = 1 dan t ≥ 0

3. Kondisi Robbin k

∂T = h(T − T f ) ∂x

pada x = 0 dan t ≥ 0

Turunan variabel terikat (T) diketahui sebagai fungsi dari variabel terikat itu sendiri. Solusi PDP

Solusi yang paling sederhana dan gampang untuk diterapkan yaitu dengang metode penghampiran selisih terhingga (finite difference approximationi).



Penghampiran Selisih Terhingga

1. Penghampiran selisih terpusat ∂ui 1 = ( ui +1 − ui −1 ) ∂x 2Δx ∂ 2ui 1 = 2 ( ui +1 − 2ui + ui −1 ) 2 ∂x Δx 2 ∂ ui , j 1 = ( ui+1, j +1 − ui −1, j +1 − ui+1, j −1 + ui−1, j −1 ) ∂y∂x 4ΔxΔy

2. Penghampiran selisih maju ∂ui 1 = ( ui +1 − ui ) ∂x Δx ∂ 2ui 1 = 2 ( ui + 2 − 2ui +1 + ui ) 2 ∂x Δx 2 ∂ ui , j 1 = ( ui +1, j +1 − ui, j +1 − ui+1, j + ui, j ) ∂y∂x 4ΔxΔy

3. Penghampiran selisih mundur ∂ui 1 = ( ui − ui −1 ) ∂x Δx ∂ 2ui 1 = 2 ( ui − 2ui −1 + ui − 2 ) 2 ∂x Δx 2 ∂ ui , j 1 = ( ui, j − ui, j −1 − ui−1, j + ui −1, j −1 ) ∂y∂x ΔxΔy

Diskretisasi Persamaan Diferensial

Dalam menyelesaikan persamaan diferensial menggunakan penghampiran selisih terhingga dikenal teknik diskretisasi. Penjelasannya diberikan berdasarkan contoh soal sebagai berikut : Kasus 9: Selesaikan persamaan differensial di bawah ini, kemudian petakan harga x dan y pada koordinat kartesius.



d2y = 6x + 4 dx 2 x = 0 → y =1 x =1→ y =1 Rentang Integrasi x = 0 s/d 1

Jawaban: Rentang integrasi x = 0 s.d 1 didiskretisasikan menjadi 10 bagian (N = 10). N = 10 1 Δx = = 0.1 N

∆x 1

2

3

4

5

6

7

8

9

10

x0

x1

x2

x3

x4

x5

x6

x7

x8

x9

x10

0

1

2

3

4

5

6

7

8

9

10

Menggunakan penghampiran selisih terpusat. ∂ 2 yi 1 = 2 ( yi +1 − 2 yi + yi −1 ) 2 ∂x Δx Substitusikan penghampiran selisih terpusat itu ke persamaan diferensial. ∂2 y − 6x + 4 = 0 ∂x 2 menghasilkan 1 ( yi +1 − 2 yi + yi −1 ) − 6 xi − 4 = 0 Δx 2



Untuk i = 1 1 ( y2 − 2 y1 + 1) − 6(0.1) − 4 = 0 Δx 2

Untuk i = 2 s.d. 8 1 ( yi +1 − 2 yi + yi −1 ) − 6 xi − 4 = 0 Δx 2

SISTEM PERSAMAAN LINIER

Untuk i = 9 1 (1 − 2 y9 + y8 ) − 6(0.9) − 4 = 0 Δx 2

Transformasi sistem persamaan linier di atas menjadi bentuk matrik sbb: 2 ⎡ −2 1 0 0 0 0 0 0 0 ⎤ ⎡ y1 ⎤ ⎡{6(0.1) + 4} Δx − 1⎤ ⎢ 1 −2 1 0 0 0 0 0 0 ⎥ ⎢ y ⎥ ⎢ 6(0.2) + 4 Δx 2 ⎥ } ⎥ ⎢ ⎥⎢ 2⎥ ⎢ { ⎢ 0 1 −2 1 0 0 0 0 0 ⎥ ⎢ y3 ⎥ ⎢ {6(0.3) + 4} Δx 2 ⎥ ⎢ ⎥⎢ ⎥ ⎢ 2 ⎥ ⎢ 0 0 1 −2 1 0 0 0 0 ⎥ ⎢ y4 ⎥ ⎢ {6(0.4) + 4} Δx ⎥ ⎢ 0 0 0 1 −2 1 0 0 0 ⎥ ⎢ y5 ⎥ = ⎢ {6(0.5) + 4} Δx 2 ⎥ ⎥ ⎢ ⎥⎢ ⎥ ⎢ 2 ⎢ 0 0 0 0 1 −2 1 0 0 ⎥ ⎢ y6 ⎥ ⎢ {6(0.6) + 4} Δx ⎥ ⎢ 0 0 0 0 0 1 −2 1 0 ⎥ ⎢ y ⎥ ⎢ {6(0.7) + 4} Δx 2 ⎥ ⎥ ⎢ ⎥⎢ 7⎥ ⎢ 2 ⎢ 0 0 0 0 0 0 1 −2 1 ⎥ ⎢ y8 ⎥ ⎢ {6(0.8) + 4} Δx ⎥ ⎢⎢ 0 0 0 0 0 0 0 1 −2 ⎥⎥ ⎢ y ⎥ ⎢ 6(0.9) + 4 Δx 2 − 1⎥ } ⎥⎦ ⎣ ⎦ ⎣⎢ 9 ⎦⎥ ⎢⎣{

A

y

C

Harga y dapat dihitung dengan metode yang telah dipelajari pada bagian sistem persamaan linier.

Ay = C y = A −1C Berikut ini kita hitung harga vektor y dalam m-file MATLAB.



kasus9.m clear clc %Diskretisasi terhadap x N=10; dx=1/N; x =[0:dx:1]' %Membuat matrik A koefisien y A = diag(-2*ones(N-1,1))+diag(ones(N-2,1),1) + diag(ones(N-2,1),-1) %Vektor konstanta C C = (6*[dx:dx:x(end)-dx]+4)*dx^2 C(1)=C(1)-1 C(end)=C(end)-1 %Menghitung harga y y=inv(A)*C' y=[1;y;1] %Membuat kurva x-y plot(x,y) xlabel('x') ylabel('y') grid on

Eksekusi program kasus9.m Masukan dan hasil di Command Window : >> kasus9 y= 1.0000 0.7210 0.4880 0.3070 0.1840 0.1250 0.1360 0.2230 0.3920 0.6490 1.0000



Dari hasil perhitungan MATLAB di atas dapat dibuat grafik dan hasil lengkap untuk persoalan ini sbb: x

y

0

1

0.1

0.721

0.2

0.488

0.3

0.307

0.4

0.184

0.5

0.125

0.6

0.136

0.7

0.223

0.8

0.392

0.9

0.649

1

1

1 0.9 0.8 0.7

y

0.6 0.5 0.4 0.3 0.2 0.1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1



Tugas 9 Menyelesaikan persamaan differensial dengan penghampiran selisih terhingga (diskretisasi)

Nomor 1 Dengan menggunakan penghampiran selisih terhingga terpusat selesaikan persamaan diferensial sbb:

d2y = y+2 dx 2 y (0) = y (1) = 1 Rentang Integrasi = 0 s/d 1 Sertakan pula kurva x,y diagram kartesiusnya.

Semidiskretisasi Persamaan Diferensial Parsial

Di muka telah dipaparkan penjelasan mengenai diskretisasi untuk menyelesaikan persamaan diferensial. PDP dapat diselesaikan juga dengan menggunakan teknik diskretisasi, sayangnya penerapan pada PDP lebih rumit dan melibatkan banyak angka. Sebagai penyederhanaan dapat digunakan teknik semidiskretisasi. Penjelasannya akan diberikan berdasarkan contoh soal sbb: kasus10 Sebuah lembaran plastik dengan tebal 1 cm mula-mula bertemperatur 25 oC diletakan diantara pelat yang bertemperatur 120 oC. Diketahui massa jenis plastik 900 kg/m3, konduktivitas termal 0.13 W/moC, dan kalor jenis 1670 J/kgoC . Pemodelan matematis untuk kasus konduksi tak tunak adalah sbb: ∂ 2T ∂T α 2 = ∂x ∂t

α=

k ρcp

Dengan metode numeris evaluasi temperatur rata-rata plastik 10 menit kemudian?



Jawaban: ∂ 2T ∂T = ∂x 2 ∂t T ( 0, t ) = 1200 C

α

Kondisi Dirichlet

T (1, t ) = 1200 C T ( x = 0 :1, 0 ) = 250 C

α=

k (0.13) = 8.6494x10−8 m 2 / s = ρ c p (900)(1670)

= 5.1896x10−2 cm 2 / menit

Diskretisasi dilakukan terhadap x.

N = 20 1− 0 Δx = = 0.05 20

∆x 1 0

0.05

3

2 0.1

0.15

4

5 0.2

6 0.25

8

7 0.3

0.35

9 0.4

10 0.45

11 0.5

0.55

12

13 0.6

x 15

14 0.65

0.7

0.75

17

16 0.8

18

0.85

Menggunakan penghampiran selisih terpusat turunan kedua T terhadap x, ∂ 2Ti 1 = 2 (Ti +1 − 2Ti + Ti −1 ) 2 ∂x Δx PDP tersebut akan menjadi. ∂T α = 2 (Ti +1 − 2Ti + Ti −1 ) ∂t Δx


19 0.9

0.95

20 1


Untuk i = 1

∂T α = 2 (T2 − 2T1 + 120 ) ∂t Δx

Untuk i = 2 : 18 ∂T α = 2 (Ti +1 − 2Ti + Ti −1 ) ∂t Δx

Untuk i = 19

SISTEM PERSAMAAN DIFERENSIAL BIASA

∂T α = 2 (120 − 2T19 + T18 ) ∂t Δx

Solusi sistem persamaan diferensial biasa telah dibahas pada bagian sebelumnya. Rutin yang akan digunakan untuk sistem PDB di termaksud adalah ode23 MATLAB. Berikut ini pemrogramannya.

taktunak.m

Pemrograman MATLAB

function dTdt=taktunak(t,T) N=20; dx=1/N; a=5.1896e-2; %diffusivitas termal,cm2/menit % untuk i = 1 dTdt(1,:) = a/(dx^2)*(T(2)-2*T(1)+120); % untuk i = 2 s.d 18 for i = 2:18 dTdt(i,:) = a/(dx^2)*(T(i+1)-2*T(i)+T(i-1)); end % untuk i = 19 dTdt(19,:) = a/(dx^2)*(120-2*T(19)+T(18));



kasus10.m

Pemrograman MATLAB

clear clc % Nilai awal dan rentang integrasi. tspan=[0:1:10]; To = 25*ones(1,19); % Fungsi pengintegrasi. [t,T] = ode23(‘taktunak',tspan,To); % Menampilkan hasil perhitungan. x=[0:1/20:1]' t T0=120*ones(length(tspan),1); T0(1)=(25+120)/2 % Harga T pada t=0,x=0, diambil rata-rata. T20=120*ones(length(tspan),1); T20(1)=(25+120)/2 % Harga T pada t=0,x=1, diambil rata-rata. T=[T0 T T20] % Memetakan hasil pd diagram kartesius 3-D. figure(1) surf(x,t,T) xlabel('x') ylabel('t') zlabel('T') colorbar shading interp figure(2) pcolor(x,t,T) xlabel('x') ylabel('t') zlabel('T') colorbar shading interp % Menghitung rata-rata suhu plastik pd menit ke 10 [b k]=size(T); z=T(b,:); Tslab=mean(z)

Eksekusi program kasus10.m &runkasus10.m. Masukan dan hasil di Command Window : >>runkasus10 Tslab = 119.5598



Gambar 3-D. Profil temperatur plastik sepanjang x pada interval waktu t



Gambar 3-D semu. Profil temperatur plastik sepanjang x pada interval waktu t

Tugas 10: Menyelesaikan PDP dengan penghampiran selisih terhingga terpusat (semi diskretisasi)

Nomor 1 Suatu fenomena difusi-konveksi dapat dideskripsikan dengan PDP berikut ini: ∂θ ∂ 2θ ∂θ 0 < x < 1, = 2 −λ ∂t ∂x ∂x t >0 θ (0, t ) = 1, ∂θ (1, t ) = 0, t > 0 ∂x θ ( x, 0) = 0, 0 < x <1

t>0

Jika harga λ = 25, selesaikan PDP diatas untuk rentang t = 0 s.d 5. Buat pula gambar 3-D θ,t,x pada koordinat kubus (kartesius).


Daftar Pustaka

Daftar Pustaka 1. Chapra, Steven C. & Canale, Raymond P., “Numerical Methods for Engineers”, 1985, diterjemahkan ke dalam bahasa Indonesia dengan judul “Metode Numerik untuk Teknik” oleh UI-Press, Jakarta, 1991. 2. Constantinides A. & Mostoufi N., “Numerical Methods for Chemical Engineers with MATLAB Applications”, Prentice Hall, New Jersey, 1983. 3. Cutlip, Michael B. & Shacham Mordechai, “Problem Solving in Chemical Engineering with Numerical Methods”, Prentice Hall, New Jersey, 1999. 4. Hanna, Owen T. & Sandal, Orville C., “Computaional Methods In Chemical Engineering”, Prentice Hall, New Jersey, 1999. 5. Lindfield G. & Penny J., “Numerical Method Using MATLAB”, Elis Horwood, London, 1995. 6. Riggs, James B., “An Introduction to Numerical Methods for Chemical Engineers”, Texas Tech University Press, Texas, 1988. 7. Sediawan, W.B. & Prasetya A., “Pemodelan Matematis dan Penyelesaian Numeris dalam Teknik Kimia”, Andi, Yogyakarta, 1997.

KOMPUTASI PROSES TEKNIK KIMIA MENGGUNAKAN MATLAB

Recommend Documents