Metode Numerik

Metode Numerik Imam Fachruddin (Departemen Fisika, Universitas Indonesia)

Daftar Pustaka: •

P. L. DeVries, A First Course in Computational Physics (John Wiley & Sons, Inc., New York, 1994)

•

W. H. Press, et. al., Numerical Recipes in Fortran 77, 2nd Ed. (Cambridge University Press, New York, 1992) (online / free download: http://www.nrbook.com/a/bookfpdf.php)

•

R. H. Landau & M. J. Páez, Computational Physics: Problem Solving with Computers (John Wiley & Sons, Inc., New York, 1997)

•

S. E. Koonin, Computational Physics (Addison-Wesley Publishing Co., Inc., Redwood City, 1986)

Isi •

akar fungsi

•

solusi sistem persamaan linear

•

fitting dengan least square

•

interpolasi

•

integrasi

•

persamaan differensial

Akar Fungsi f(x) = 0

x=? akar fungsi f(x)

Contoh:

x-

1 =0 x

x2 = 1

3x2 = 6 - 7x

x = 1 dan -1

3x2 + 7x - 6 = (3x - 2)(x + 3) = 0

Pada dua contoh di atas akar fungsi dapat dicari secara analitik. Secara umum, tidak selalu begitu keadaannya.

x = 2/3 dan -3

Problem: Sebuah lampu dipasang di pinggir sebuah piringan berjari-jari 10 cm. Sebuah plat bercelah sempit diletakkan di dekat piringan itu. Tepat di belakang celah itu dipasang sebuah sensor cahaya yang menghadap tegak lurus ke celah. Piringan diputar konstan 1 rad/s dan plat beserta sensor digeser lurus konstan 10 cm/s. Saat ini posisi celah dan lampu seperti pada gambar 1. Kapan sensor cahaya menerima cahaya terbanyak? Sensor menerima cahaya terbanyak pada saat posisi lampu dan celah membentuk garis tegak lurus terhadap plat, seperti pada gambar 2. x = r cos (ωt) = vt r

lampu

gambar 1

gambar 2 celah

sensor

ω plat v

cos (t) = t

?

Plot cos(x) dan x: Grafik ini menunjukkan bahwa cos(x) = x pada x sedikit kurang dari 0.75.

Bisakah lebih akurat lagi?

Cari secara numerik akar fungsi dari f(x) = cos(x) - x

Bisection Prinsip: Kurung akar fungsi di antara dua batas, lalu paruh batas itu terus menerus sampai batas itu sedemikian sempit dan dengan demikian lokasi akar fungsi diketahui dengan keakuratan tertentu. Langkah: 1. Perkirakan akar fungsi (bisa dengan cara memplot fungsi). 2. Tentukan batas awal yang mengurung akar fungsi.

akar fungsi

a

df e c

b

3. Belah dua daerah berisi akar fungsi itu. 4. Tentukan daerah yang berisi akar fungsi. 5. Ulangi langkah 3 dan 4 sampai dianggap cukup. 6. Tentukan akar fungsi.

Batas e, f atau nilai di tengahnya bisa dipilih sebagai akar fungsi.

•

Menentukan daerah yang berisi akar fungsi: Jika z merupakan akar fungsi, maka f(x < z) dan f(x > z) saling berbeda tanda.

a

c

b

f(x)

f(a)*f(c) negatif, berarti di antara a & c ada akar fungsi.

z

f(b)*f(c) positif, berarti di antara b & c tidak ada akar fungsi •

Menentukan kapan proses pencarian akar fungsi berhenti: Proses pencarian akar fungsi dihentikan setelah keakuratan yang diinginkan dicapai, yang dapat diketahui dari kesalahan relatif semu.

kesalahan relatif semu =

perkiraan sebelum - perkiraan berikut perkiraan berikut

x

Kesalahan kesalahan mutlak = | perkiraan – nilai sebenarnya |

kesalahan relatif =

perkiraan – nilai sebenarnya nilai sebenarnya

Dalam perhitungan numerik, nilai sebenarnya justru sering tidak diketahui, yang didapat hanya perkiraan terbaik. Karena perkiraan langkah berikut dianggap lebih akurat, yaitu lebih mendekati nilai sebenarnya, maka kesalahan yang dihitung yaitu: kesalahan mutlak semu = | perkiraan sebelum – perkiraan berikut |

kesalahan relatif semu =

perkiraan sebelum - perkiraan berikut perkiraan berikut

False Position Prinsip: Di sekitar akar fungsi yang diperkirakan, anggap fungsi merupakan garis lurus. Titik tempat garis lurus itu memotong garis nol ditentukan sebagai akar fungsi.

f(x)

akar fungsi sebenarnya

x

akar fungsi yang diperoleh

garis lurus sebagai pengganti f(x)

f(x) b x

c a

p(x)

Diperoleh:

 x −b  x−a p(x) =  f(a) +  f(b)  a −b  b−a 

p(c) = 0

c=

af(b) − bf(a) f(b) − f(a)

Langkah: 1. Perkirakan akar fungsi (bisa dengan cara memplot fungsi).

akar fungsi

f(x)

b

2. Tentukan batas awal yang mengurung akar fungsi. 3. Tarik garis lurus penghubung nilai fungsi pada kedua batas, lalu cari titik potongnya dengan garis nol. 4. Geser salah satu batas ke titik potong itu, sementara batas lain tidak berubah. Ulangi langkah 3.

a

f(x) b

5. Ulangi langkah 4 sampai dianggap cukup. 6. Titik potong garis nol dan garis lurus yang terakhir dinyatakan sebagai akar fungsi.

x

c

c a

c=

af(b) − bf(a) f(b) − f(a)

x

f(x)

c=

af(b) − bf(a) f(b) − f(a)

b c Metode false position juga menggunakan dua batas seperti metode bisection. Namun, berbeda dari metode bisection, pada metoda false position hanya satu batas yang berubah.

x

a

f(x)

Pada contoh sebelum ini, batas a berubah sementara batas b tetap. Pada contoh berikut terjadi sebaliknya.

b c a

x

Menghitung akar fungsi dengan metode false position, menggunakan a dan b sebagai batas awal: •

jika batas a tetap, batas b berubah:

xi+1 = •

(i = 0, 1, 2, ...;

x0 = b)

jika batas b tetap, batas a berubah:

xi+1 = •

af(xi ) − xif(a) f(xi ) − f(a)

bf(xi ) − xif(b) f(xi ) − f(b)

(i = 0, 1, 2, ...;

x0 = a )

kesalahan relatif semu:

∆ rel =

xi − xi+1 xi+1

Penghitungan dihentikan jika kesalahan relatif semu sudah mencapai / melampaui batas yang diinginkan.

Newton-Raphson Prinsip: Buat garis singgung kurva f(x) di titik di sekitar akar fungsi. Titik tempat garis singgung itu memotong garis nol ditentukan sebagai akar fungsi.

akar fungsi yang diperoleh

f(x)

x a akar fungsi sebenarnya

p(x) = garis singgung kurva f(x) di titik f(a)

f(x) c

x

a

Diperoleh:

p(x)

p(x) = f(a) + (x − a)f'(a) (f’(a) turunan pertama f(x) pada x = a)

p(c) = 0

c=a−

f(a) f'(a)

Langkah:

f(x)

1. Perkirakan akar fungsi.

c

2. Buat garis singgung pada titik sesuai akar fungsi yang diperkirakan itu, lalu cari titik potongnya dengan garis nol.

x

a akar fungsi sebenarnya

3. Titik potong itu merupakan perkiraan akar fungsi baru. 4. Ulangi langkah 2 dan 3 sampai dianggap cukup. 5. Titik potong garis nol dan garis singgung kurva yang terakhir dinyatakan sebagai akar fungsi.

f(x) a c

c=a−

f(a) f'(a)

x

1

f(x) Contoh perkiraan akar fungsi awal yang baik
perkiraan akar fungsi makin mendekati akar fungsi sebenarnya.

x

2 1 f(x) Contoh perkiraan akar fungsi awal yang buruk
perkiraan akar fungsi makin menjauhi akar fungsi sebenarnya.

x

2

Menghitung akar fungsi dengan metode Newton-Raphson:

xi+1 = xi −

f(xi ) f'(xi )

(i = 0, 1, 2, ...;

x0 = a )

kesalahan relatif semu:

∆ rel =

xi − xi+1 xi+1

Penghitungan dihentikan jika kesalahan relatif semu sudah mencapai / melampaui batas yang diinginkan.

Secant Kembali ke metode False Position, untuk contoh batas b tetap, akar fungsi dicari sebagai berikut:

x1 =

bf(x0 ) − x0 f(b) f(x0 ) − f(b)

x2 =

bf(x1 ) − x1f(b) f(x1 ) − f(b)

x3 =

bf(x2 ) − x2f(b) f(x2 ) − f(b)

Pada metode Secant, batas tidak dijaga tetap, melainkan berubah. Akar fungsi dicari sebagai berikut:

x1 =

bf(x0 ) − x0 f(b) f(x0 ) − f(b)

x2 =

bf(x1 ) − x1f(b) f(x1 ) − f(b)

x3 =

Jadi, mulai dari i = 3, akar fungsi dihitung dengan: xi =

x1f(x2 ) − x2f(x1 ) f(x2 ) − f(x1 ) xi-2f(xi-1 ) − xi-1f(xi-2 ) f(xi-1 ) − f(xi-2 )

f(x)

f(x) b

b x

x

x1

a

x1

x2 II

I

secant

false position f(x)

x3

f(x) b

x2

x

x

x1

x3 III

III

x2

Akar fungsi pada metode Secant untuk i = 1, 2 bisa dihitung dengan metode yang lain atau ditebak. Mulai i = 3, akar fungsi dihitung dengan rumus:

x f(x ) − xi-1f(xi-2 ) xi = i-2 i-1 f(xi-1 ) − f(xi-2 )

-1

 f(xi-1 ) − f(xi-2 )   f(xi-1 ) xi = xi-1 −   xi-1 − xi-2 

Yang menarik, jika i makin besar, maka beda antar dua akar fungsi yang berturutan semakin kecil, sehingga

f(xi-1 ) − f(xi-2 ) df(xi-1 ) ≅ = f'(xi-1 ) xi-1 − xi-2 dxi-1

xi ≅ xi-1 −

f(xi-1 ) f'(xi-1 )

Dengan begitu, metode Secant menyerupai metode Newton-Raphson. Jika turunan fungsi f(x) sulit diperoleh / dihitung, maka metode Secant menjadi alternatif yang baik bagi metode Newton-Raphson. Kesalahan relatif semu dihitung sama seperti pada metode False Position atau Newton-Raphson.

Kecepatan Konvergensi Pencarian akar fungsi dimulai dengan perkiraan akar fungsi yang pertama, lalu diikuti oleh perkiraan berikutnya dan seterusnya sampai perkiraan yang terakhir, yang kemudian dinyatakan sebagai akar fungsi hasil perhitungan tersebut. Proses itu harus bersifat konvergen yaitu, selisih perkiraan sebelum dari yang setelahnya makin lama makin kecil. Setelah dianggap cukup, proses pencarian akar fungsi berhenti.

x2 − x1 > x3 − x2 > x4 − x3 ... xn − xn −1 ≤ ε (ε = bilangan kecil) Kecepatan konvergensi sebuah proses yaitu, kecepatan proses itu untuk sampai pada hasil akhir.

Contoh pencarian akar fungsi dengan metode Bisection: akar fungsi

x a

b

x1 x3

x2

x4

Jika εi ≡ xi+1 − xi , maka dari gambar diperoleh:

ε1 = x2 − x1 , ε2 = x3 − x2 , ε3 = x4 − x3 ε2 = 21 ε1 , ε3 = 21 ε2 Kecepatan konvergensi bersifat linear:

εi+1 = 21 εi

Pada metode False Position, Newton-Raphson dan Secant akar fungsi dicari dengan rumus yang bentuknya serupa: -1

 f(xi ) − f(a)   f(xi ) = xi −   xi − a 

False Position:

xi+1

Newton-Raphson:

xi+1 = xi −

(atau a diganti b)

f(xi ) f'(xi ) -1

Secant:

 f(xi ) − f(xi-1 )   f(xi ) xi+1 = xi −   xi − xi-1 

Mengingat dengan berjalannya proses pencarian akar fungsi rumus pada metode False Position dan terlebih lagi Secant semakin mendekati rumus pada metode Newton-Raphson, maka akan dibahas kecepatan konvergen pada metode Newton-Raphson.

xi+1 = xi −

f(xi ) f'(xi )

εi ≡ xi − xi+1 =

f(xi ) f'(xi )

εi+1 =

f(xi+1 ) =? f'(xi+1 )

ekspansi deret Taylor:

f(xi+1 ) = f(xi − εi ) = f(xi ) − εif'(xi ) + 21 εi2f''(xi ) − ... f'(xi+1 ) = f'(xi − εi ) = f'(xi ) − εif''(xi ) + ...

f(xi ) − εif'(xi ) + 21 εi2f''(xi ) − ... εi+1 = f'(xi ) − εif''(xi ) + ... f(xi ) − εif'(xi ) + 21 εi2f''(xi ) ≅ f'(xi ) ≅

f''(xi ) 2 εi 2f'(xi )

Kecepatan konvergensi pada metode NewtonRaphson (kira-kira demikian juga False Position dan Secant) bersifat kurang lebih kuadratik:

εi+1 ≅

f''(xi ) 2 εi 2f'(xi )

Dengan begitu, metode metode Newton-Raphson, False Position dan Secant lebih cepat dari metode Bisection.

Contoh hasil pencarian akar fungsi untuk soal cos(x) = x:

metode

akar

f(akar)

jumlah langkah

Bisection

0.7390795

9.3692161E-06

12

False Position

0.7390851

-7.7470244E-09

3

Newton-Raphson

0.7390851

-7.7470244E-09

4

Secant

0.7390851

-7.7470244E-09

3

Keterangan:

•

Pencarian akar berhenti jika kesalahan relatif semu sama atau kurang dari 1.0E-05.

•

Batas awal kiri dan kanan untuk metode Bisection, False Position dan Secant 0.72 dan 0.75.

•

Perkiraan akar fungsi pertama untuk metode Newton-Raphson 0.72.

Solusi Sistem Persamaan Linear Sistem persamaan linear:

a11x1

+ a12x2

+ a13x3 ⋯ + a1n xn

= b1

a21x1 a31x1

+ a22x2 + a32x2

+ a23x3 ⋯ + a2n xn + a33x3 ⋯ + a3n xn

= b2 = b3

⋮

⋮

an1x1

+ an2x2

⋮

⋱

⋮

+ an3x3 ⋯ + ann xn

⋮ = bn

n buah persamaan dengan n buah unknown

xj

aij dan bi

diketahui i, j = 1, 2, …, n

xj = ?

Soal:

Jawab:

2x − 3y + 2z = −6

(1)

− x + 2y − 3z = 2

(2)

x+ y− z =0

(3)

2x −

2x −

3y + 2z = −6

(1)

0.5y − 2z = −1

(2)

2.5y − 2z = 3

(3)

3y + 2z = −6

(1)

0.5y − 2z = −1

(2)

8z = 8

eliminasi x: pers. (2) + 0.5 pers. (1) pers. (3) – 0.5 pers. (1)

eliminasi y: pers. (3) – 5 pers. (2)

(3)

z =1 − 1 + 2z =2 0.5 − 6 + 3y − 2z = −1 x= 2 y=

3 persamaan dan 3 unknown

substitusi mundur: pers. (3)
mencari z pers. (2)
mencari y pers. (1)
mencari x

Dalam bentuk matriks: Soal:

2  x   − 6   2 −3      2 − 3  y  =  2   −1  1 1 − 1  z   0  

Jawab:

2  x   − 6  2 −3       0 0.5 − 2  y  =  − 1   0 2.5 − 2  z   3       2  x   − 6  2 −3       0 0.5 − 2  y  =  − 1  0 0 8  z   8  

z =1 − 1 + 2z =2 y= 0.5 − 6 + 3y − 2z = −1 x= 2

Eliminasi Gauss Metode Eliminasi Gauss mencari solusi sebuah sistem persamaan linear dengan cara seperti ditunjukkan pada contoh sebelum ini:

 a11   a21 a  31  ⋮ a  n1

a12 a22 a32 ⋮ an2

a13 ⋯ a1n  x1   b1      a23 ⋯ a2n  x2   b2  a33 ⋯ a3n  x3  =  b3      ⋮ ⋱ ⋮  ⋮   ⋮  an3 ⋯ ann  xn   bn 

aij(0) = aij , bi(0) = bi (k) ij

a

(k) i

b

(k -1) ij

aik(k-1) (k-1) − (k-1) akj akk

(k −1) i

aik(k-1) (k −1) − (k-1) bk akk

=a

=b

(k = 1, ..., n − 1;

a   0   0  ⋮   0

(0) 11

(0) 12 (1) 22

a a 0 ⋮ 0

(0) 13 (1) 23 (2) 33

a a a ⋮ 0

⋯ a  x1   b      ⋯ a  x2   b     ⋯ a  x3  =  b    ⋱ ⋮  ⋮   ⋮    (n-1)  (n-1)  x ⋯ ann  n   bn  (0) 1n (1) 2n (2) 3n

(0) 1 (1) 2 (2) 3

i = k + 1, ..., n; j = k, ..., n)

aij(m) , bi(m) ≡ aij , bi pada langkah ke m halaman berikut

Substitusi mundur:

 a11(0)   0   0  ⋮   0

(0) a12 (1) a22 0

⋮ 0

(0) a13 ⋯ a1n(0)  x1   b1(0)     (1)  (1) (1) a23 ⋯ a2n  x2   b2  (2) (2)  x3  =  b3(2)  a33 ⋯ a3n     ⋮ ⋱ ⋮  ⋮   ⋮    (n-1)  (n-1)  x 0 ⋯ ann  n   bn 

bn(n-1) xn = (n-1) ann (n- j-1) n-j

b xn − j =

n

∑a

−

k =n- j+1 (n- j-1) n - j,n - j

a

(n- j-1) n - j,k k

x

(j = 1, ..., n − 1)

A

X = B

atau

AX = B

Jadi, metode Eliminasi Gauss terdiri dari dua tahap: 1. triangulasi: mengubah matriks A menjadi matriks segitiga (matriks B dengan begitu juga berubah)

=

=

2. substitusi mundur: menghitung x mengikuti urutan terbalik, dari yang terakhir ( xn ) sampai yang pertama ( x1 )

LU Decomposition

A

X = B

atau

AX = B

Pada metode LU Decomposition, matriks A ditulis ulang sebagai perkalian matriks L dan U (matriks A diurai menjadi matriks L dan U). Matriks L dan U merupakan matriks segitiga. Matriks B tidak berubah, karena matriks A tidak berubah, melainkan hanya ditulis ulang.

A

X = B

U L

X = B

Langkah: 1. Cari matriks L dan U sehingga A = LU. Matriks B tetap.

A

=

U

U

L

L

X = B

2. Definisikan sebuah matriks kolom baru, misalnya Y, yaitu Y = UX, sehingga LY = B. Lalu hitung y dengan substitusi maju (mulai dari y1 sampai yn ).

U

Y =

X

L

Y = B

3. Hitung x dengan substitusi mundur (mulai dari xn sampai x1 ). U

X = Y

Jelas bahwa metode LU Decomposition pada prinsipnya sama dengan metode Eliminasi Gauss: matriks U merupakan hasil triangulasi matriks A, yang juga mengakibatkan B berubah menjadi Y.

Mencari matriks L dan U:

 l11 0 0  l21 l22 0 l l l  31 32 33 ⋮ ⋮ ⋮   ln1 ln2 ln3

⋯ 0  ⋯ 0 ⋯ 0  ⋱ ⋮  ⋯ lnn 

 1 u12 u13   0 1 u23 0 0 1  ⋮ ⋮ ⋮ 0 0 0 

⋯ u1n   a11   ⋯ u2n   a21 ⋯ u3n  =  a31   ⋱ ⋮   ⋮  ⋯ 1   an1

a12 a22 a32 ⋮ an2

a13 a23 a33 ⋮ an3

⋯ ⋯ ⋯ ⋱ ⋯

Diperoleh:

ai1 = li1 a1j = l11u1j ai2 = li1u12 + li2 a2j = l21u1j + l22u2j ai3 = li1u13 + li2u23 + li3 a3j = l31u1j + l32u2j + l33u3j …

→ li1 = ai1 a1j → u1j = l11 → li2 = ai2 − li1u12 a2j − l21u1j → u2j = l22 → li3 = ai3 − li1u13 − li2u23 a3j − l31u1j − l32u2j → u3j = l33

(i = 1, ..., n) (j = 2, ..., n) (i = 2, ..., n) (j = 3, ..., n) (i = 3, ..., n) (j = 4, ..., n)

a1n   a2n  a3n   ⋮  ann 

Jadi, elemen matriks L dan U dicari menurut:

li1 = ai1 u1j =

(i = 1, ..., n)

a1j

(j = 2, ..., n)

l11 j−1

lij = aij − ∑ likukj

(i = j, ..., n; j = 2, ..., n)

k =1 i−1

aij − ∑likukj uij =

k =1

lii

(j = i + 1, ..., n; i = 2, ..., n - 1)

secara bergantian: 1. matriks L kolom 1, matriks U baris 1 2. matriks L kolom 2, matriks U baris 2 3. … 4. matriks L kolom (n-1), matriks U baris (n-1) 5. matriks L kolom n

Substitusi maju untuk menghitung y:

 l11 0 0  l21 l22 0 l l l  31 32 33 ⋮ ⋮ ⋮   ln1 ln2 ln3

⋯ 0  y1   b1      ⋯ 0  y2   b2  ⋯ 0  y3  =  b3      ⋱ ⋮  ⋮   ⋮      ⋯ lnn  yn   bn 

y1 =

b1 l11 i−1

bi − ∑ lij yj yi =

j=1

(i = 2, ..., n)

lii

Substitusi mundur untuk menghitung x:

 1 u12 u13   0 1 u23 0 0 1  ⋮ ⋮ ⋮ 0 0 0 

⋯ u1n  x1   y1      ⋯ u2n  x2   y2  ⋯ u3n  x3  =  y3      ⋱ ⋮  ⋮   ⋮      ⋯ 1  xn   yn 

xn = yn n

xn-i = yn −i −

∑u

n −i, j

j=n −i+1

xj (i = 1, ..., n − 1)

2  2 −3  − 6     AX = B Kembali ke soal , dengan A =  − 1 2 − 3 , B =  2  .  1  0 1 − 1    

Jawab:

A = LU

0 1  2 0  1 − 1.5     L =  − 1 0.5 0 , U =  0 1 − 4  1 2.5 8  0 0 1    

Y = UX, LY = B

UX = Y

 −1   X =  2  1  

 − 3   Y =  − 2  1  

Kasus Beberapa Sistem Persamaan Linear Pada kasus yang lebih umum bisa saja terdapat beberapa sistem persamaan linear dengan nilai B yang berlainan, namun memiliki nilai A yang sama. Dalam bentuk matriks sistem seperti ini dituliskan sebagai:

A

Keterangan: •

•

X

=

B

atau AX = B

A matriks n x n, X dan B matriks n x m, dengan m = jumlah sistem persamaan linear, n = jumlah persamaan / unknown dalam tiap sistem persamaan tersebut Tiap kolom matriks X merupakan solusi untuk kolom yang sama pada matriks B.

Langkah dan rumus pada metode Eliminasi Gauss dan LU Decomposition berlaku sama untuk kasus ini. Hanya saja, di sini matriks X dan B terdiri dari beberapa kolom, bukan hanya satu.

Contoh dua sistem persamaan linear yang memiliki nilai A sama tapi B berbeda.

 a11   a21 a  31  ⋮ a  n1

a12 a22 a32 ⋮ an2

a13 a23 a33 ⋮ an3

⋯ a1n  x11   b11      ⋯ a2n  x21   b21  ⋯ a3n  x31  =  b31      ⋱ ⋮  ⋮   ⋮  ⋯ ann  xn1   bn1 

 a11   a21 a  31  ⋮ a  n1

a12 a22 a32 ⋮ an2

 a11   a21 a  31  ⋮ a  n1

a13 ⋯ a1n  x11  a23 ⋯ a2n  x21 a33 ⋯ a3n  x31  ⋮ ⋱ ⋮  ⋮ an3 ⋯ ann  xn1

a12 a22 a32 ⋮ an2

a13 a23 a33 ⋮ an3

⋯ a1n  x12   b12      ⋯ a2n  x22   b22  ⋯ a3n  x32  =  b32      ⋱ ⋮  ⋮   ⋮  ⋯ ann  xn2   bn2 

x12   b11 b12     x22   b21 b22  x32  =  b31 b32     ⋮   ⋮ ⋮  xn2   bn1 bn2 

Metode Eliminasi Gauss: •

rumus triangulasi:

aij(0) = aij , bir(0) = bir

aij(k) = aij(k-1) − bir(k) •

(k -1) ik (k -1) kk

a a

(i, j = 1, ..., n; r = 1, ..., m)

akj(k-1) (k = 1, ..., n − 1;i = k + 1, ..., n;

(k -1) a = bir(k −1) − ik(k-1) bkr(k −1) akk

aij(m) , bir(m) ≡ aij , bir pada langkah ke m

j = k, ..., n;r = 1, ..., m)

rumus substitusi mundur:

bnr(n-1) xnr = (n-1) ann (n- j-1) n- j,r

b xn − j,r =

(r = 1, ..., m) n

∑a

−

k =n - j+1 (n- j-1) n - j,n - j

a

(n- j-1) n- j,k kr

x

(j = 1, ..., n − 1; r = 1, ..., m)

Metode LU Decomposition: •

rumus substitusi maju untuk menghitung y (kini Y matriks n x m):

y1r =

b1r l11

(r = 1, ..., m) i−1

bir − ∑ lij yjr yir =

•

j=1

(i = 2, ..., n; r = 1, ..., m)

lii

rumus substitusi mundur untuk menghitung x:

xnr = ynr

(r = 1, ..., m) n

xn-i,r = yn −i,r −

∑u

n −i, j

j=n −i+1

xjr (i = 1, ..., n − 1;r = 1, ..., m)

Perhatikan metode LU Decomposition, anggap matriks L dan U telah diperoleh. Jika kemudian terdapat lagi sistem persamaan linear dengan A sama dan B berbeda, maka matriks L dan U yang telah diperoleh itu bisa langsung dipakai untuk sistem persamaan yang baru tersebut. Kini perhatikan metode Eliminasi Gauss, anggap triangulasi matriks A sudah dikerjakan. Jika kemudian terdapat lagi sistem persamaan linear dengan A sama dan B berbeda, maka hasil triangulasi matriks A yang sudah diperoleh tidak dapat dipakai untuk sistem persamaan yang baru. Untuk sistem yang baru tersebut proses triangulasi matriks A harus dilakukan lagi dari awal. Hal ini disebabkan, matriks B harus berubah mengikuti proses triangulasi matriks A, sementara proses penguraian matriks A menjadi matriks L dan U tidak melibatkan matriks B.

Catatan: Dalam rumus-rumus baik pada metode Eliminasi Gauss maupun LU Decomposition terdapat pembagian oleh elemen diagonal matriks yaitu, oleh elemen diagonal matriks A pada metode Eliminasi Gauss, dan elemen diagonal matriks L pada metode LU Decomposition. Jika secara kebetulan elemen diagonal itu nol, maka akan timbul error. Karena itu, pada setiap langkah dalam proses triangulasi matriks A (metode Eliminasi Gauss) atau pencarian matriks L dan U (metode LU Decomposition) perlu dilakukan pemeriksaan, apakah elemen matriks A atau L yang bersangkutan sama dengan nol. Jika bernilai nol, maka baris berisi elemen diagonal nol itu harus ditukar dengan salah satu baris setelahnya, sehingga elemen diagonal menjadi bukan nol. Perubahan baris pada matriks A (metode Eliminasi Gauss) harus disertai perubahan baris yang sama pada matriks B. Perubahan baris pada matriks L (metode LU Decomposition) harus disertai perubahan baris yang sama pada matriks A dan B.

Soal:

Jawab:

1 3  x1   2   2 −4      2 3 − 2  x2   2   −1 =     3 −4  1 2 x3 2       1 − 3 −1 5  x4   2  

baris 2 ditukar dengan baris 3

1 3  x1   2  2 − 4      0 3.5 − 0.5  x2   3  0 =    0  2 − 0.5 − 2.5 x3 −1       0 − 1 − 1 .5  x   1 3 . 5   4   

1 3  x1   2  2 − 4      2 − 0.5 − 2.5  x2   − 1 0 = 0 0 3.5 − 0.5  x3   3       0 0 0 2  x4   2  

x4 = 1 x3 =

3 + 0.5x4 =1 3.5

1 3  x1   2  2 − 4      2 − 0.5 − 2.5  x2   − 1 0 =    0  0 3.5 − 0.5 x3 3       0 − 1 − 1 .5  x   1 3 . 5   4   

1 3  x1   2  2 − 4      2 − 0.5 − 2.5  x2   − 1  0 = 0 0 3.5 − 0.5  x3   3       0     0 − 1.75 2.25  x4   0.5  

− 1 + 0.5x3 + 2.5x4 =1 2 2 + 4x2 − x3 − 3x4 =1 x1 = 2 x2 =

1 3  2  2 −4     2 3 − 2  2  −1 B=  A=  2 3 −4 1 2      2  1 − 3 −1 5     1 3  2 −4   1 2  3 −4 A= −1 2 3 − 2    1 − 3 −1  5    2    2 B=  2    2  

 2 0 0 0    − 1 l22 0 0  L= 3 l32 l33 0     1 l l43 l44  42 

 2 0 0 0    3 2 0 0 L= − 1 0 l33 0     1 −1 l l44  43 

0.5 1.5   1 −2   1 − 0.25 − 1.25  0 U= 0 0 1 u34    0  0 0 1  

0  2 0  ¨ 2 0  3 2 s i r 3 L= ba ris − 1 0 3.5  ba  1 − 1 − 1.75 

0  0 0  2 

 1 − 2 0.5 1.5    1 u23 u24  0 U= 0 0 1 u34    0  0 0 1  

 2 0 0 0   −1 0 0 0 L= 3 2 l33 0     1 −1 l l44  43  0 0  2 0   0 0  3 2 L= − 1 0 3.5 0    1 − 1 − 1.75 l  44  

0.5 1.5   1 −2   1 − 0.25 − 1.25  0 U= 0 0 1 − 1/7    0  0 0 1  

Data Fitting dengan Metode Least Square Keterangan:

f(x) p(x)

x

•

f(xi ) mewakili data; i = 1, …, N; N = jumlah data

•

p(x) merupakan fungsi yang dicocokkan (fitted) terhadap data f( xi )

Sifat fitting: tidak selalu p(xi ) = f( xi ) untuk semua xi .

Prinsip penentuan fungsi p(x): •

p(x) merupakan polinomial orde m: m

p(x) = a0 + a1x + a2x + a3x + ... + amx = ∑ ajx j 2

3

m

j= 0

(Secara umum, p(x) juga bisa merupakan polinomial bentuk yang lain seperti, polinomial Legendre.) •

Selisih antara p(x) dan f(x) untuk titik data tertentu: m

∆i = f(xi ) − p(xi ) = f(xi ) − ∑ ajxij

(i = 1, ..., N )

j= 0

•

Jumlah kuadrat selisih antara p(x) dan f(x) untuk semua titik data: m   2 j 2  S = ∑ ∆i = ∑ (f(xi ) − p(xi ) ) = ∑  f(xi ) − ∑ ajxi  i=1 i=1 i=1  j= 0  N

N

N

2

Fungsi p(x) ditentukan dengan mencari nilai aj (j = 0, …, m) yang membuat S bernilai minimum.

Titik Minimum g(x) g(a) merupakan titik minimum jika:

dg(x) dx x = a

a

2 dg(x) d g(x) = 0 dan >0 2 dx x = a dx x = a

x

Spesial: fungsi kuadratik g(x) = ax2 + bx + c dg(x) = 2ax + b dx 2 d g(x) = 2a dx 2

  g(x) memiliki satu titik minimun jika a > 0 atau  sebaliknya satu titik maksimum jika a < 0.  

S merupakan fungsi kuadratik dalam aj (j = 0, …, m): 2

m N  m    j 2 2 2j    S(a0 , ..., am ) = ∑  f(xi ) − ∑ ajxi  = ∑  ∑ aj xi + ... + f (xi )  j= 0 i=1  j= 0 i=1    N

(

N  m  k ∂S(a0 , ..., am ) j  = −2∑  f(xi ) − ∑ ajxi xi ∂ak i=1  j= 0 

N ∂ 2S(a0 , ..., am ) = 2∑ xi2k > 0 2 ∂ak i=1

)

(k = 0, ..., m)

(k = 0, ..., m)

S memiliki satu titik minimum pada nilai aj (j = 0, …, m) tertentu.

Mencari aj (j = 0, …, m): N  m  k ∂S(a0 , ..., am ) j  = −2∑  f(xi ) − ∑ ajxi xi = 0 ∂ak i=1  j= 0 

N  N j+k   ∑ xi aj = ∑ f(xi )xik ∑ j= 0  i=1 i=1  m

Definisikan:

N

(k = 0, ..., m)

(k = 0, ..., m)

N

ckj ≡ ∑ x

bk ≡ ∑ f(xi )xik

j +k i

i=1

i=1 m

maka diperoleh sebuah sistem persamaan linear:

∑c

a = bk

kj j

(k = 0, ..., m)

j= 0

dalam bentuk matrik:

C

A = B

atau

CA = B

Jadi, aj (j = 0, …, m) diperoleh sebagai solusi persamaan linear CA = B.

Contoh: Terdapat tiga data f(x) yaitu, f(1) = 30, f(2) = 70 dan f(3) = 120. Cari fungsi p(x) yang dapat melukiskan data itu. Dari data itu jelas p(x) bukan fungsi linear. Jadi, dicoba fungsi kuadratik:

f(x)

p(x)

120

p(x) = a0 + a1x + a2x2 Sistem persamaan linier untuk mencari aj :

 3 6 14  a0   220       6 14 36 a = 530   1     14 36 98  a   1390    2     3 6 14  a0   220        0 1 4  a1  =  45   0 0 1  a   5    2    Jadi,

p(x) = 5x(x + 5 )

70 30 1

2

3

 a0   0       a1  =  25  a   5   2   Cek:

p(1) = 30, p(2) = 70, p(3) = 120

! K O

x

Contoh: Kuat medan listrik E di sekitar sebuah benda berbentuk lempeng diukur pada jarak 10 cm dari pusat massanya dan arah yang bervariasi. Arah dinyatakan dalam sudut θ terhadap sumbu y yang ditetapkan sebelum pengukuran. Diperoleh data sebagai berikut: θ [derajat]

E [V/cm]

10 15 20 25 30 35 40 45 50 55

0.01794775 0.03808997 0.05516225 0.05598281 0.04795629 0.04807485 0.06273566 0.07853982 0.07395442 0.04201338

Cari fungsi p(x) yang dapat melukiskan data itu.

y

θ

E

Dicoba beberapa polinomial dengan orde berbeda, diperoleh: a0 = 8.983713484853211E - 03 a1 = 1.324478388111303E - 03

m = 3: a2 = 3.487808787880805E - 05 a3 = - 8.085809790211842E - 07 S = 1.0339E - 03

a0 = - 3.557800654975570E - 02 a1 = 1.061996221844471E - 03 a2 = 8.802185976358352E - 04

m = 5: a3 = - 5.862332690401015E - 05 a4 = 1.362046192596346E - 06 a5 = - 1.063951754163944E - 08 S = 8.1573E - 05

a0 = - 1.757260839248139E - 02

m = 9:

a1 = 1.596300085173997E - 02 a2 = - 3.402768734407800E - 03

a0 = 1.864754537649403E - 01

a3 = 3.358961098305538E - 04

a2 = 4.007658091692495E - 03

a4 = - 1.368895999268855E - 05 a5 = 1.132254508386570E - 07 a6 = 8.262829873458547E - 09 a7 = - 2.741786330789355E - 10 a8 = 3.317446724324134E - 12 a9 = - 1.459511835946927E - 14 S = 1.7528E - 11

a1 = - 4.631839872868015E - 02 a3 = - 8.985715636865594E - 05

m = 7: a4 = - 3.230489224228010E - 06 a5 = 1.912806006890119E - 07 a6 = - 3.252863805243949E - 09 a7 = 1.876184315740421E - 11 S = 3.1629E - 07

Interpolasi Keterangan:

f(x) p(x)

x

•

f(xi ) mewakili data; i = 1, …, N; N = jumlah data

•

p(x) merupakan fungsi interpolasi berdasarkan data f(xi )

Sifat interpolasi: p(xi ) = f(xi ) untuk semua xi .

Interpolasi Lagrange Digunakan p(x), suatu polinomial berorde m = N – 1, dengan N = jumlah data:

p(x) = a0 + a1x + a2x2 + ... + aN-1x N-1 ≅ f(x) Nilai ai (i = 0, …, N-1) ditetukan dengan menetapkan bahwa untuk semua titik data: (i = 1, ..., N ) p(xi ) = f(xi ) Jadi, diperoleh persamaan linear:

p(x1 ) = a0 + a1 x1 + a2 x12 + ... + aN-1 x1N-1 = f(x1 ) p(x2 ) = a0 + a1x2 + a2 x22 + ... + aN-1x2N-1 = f(x2 ) p(x3 ) = a0 + a1x3 + a2x32 + ... + aN-1x3N-1 = f(x3 ) ... p(xN ) = a0 + a1xN + a2xN2 + ... + aN-1 xNN-1 = f(xN ) dan ai (i = 0, …, N-1) diperoleh sebagai solusi dari persamaan linear itu.

p(x1 ) = a0 + a1x1 = f(x1 )

N = 2:

p(x2 ) = a0 + a1 x2 = f(x2 ) a0 = −

x2f(x1 ) − x1f(x2 ) x1 − x2

a1 =

f(x1 ) − f(x2 ) x1 − x2

 x − x2   x − x1  f(x1 ) +  f(x2 ) p(x) =   x1 − x2   x2 − x1 

p(x1 ) = a0 + a1 x1 + a2 x12 = f(x1 )

N = 3:

p(x2 ) = a0 + a1x2 + a2 x22 = f(x2 ) p(x3 ) = a0 + a1x3 + a2x32 = f(x3 ) a0 =

(x2 − x3 )x2 x3f(x1 ) + (x3 − x1 )x3x1f(x2 ) + (x1 − x2 )x1x2f(x3 ) (x2 − x3 )x12 + (x3 − x1 )x22 + (x1 − x2 )x32

(x22 − x32 )f(x1 ) + (x32 − x12 )f(x2 ) + (x12 − x22 )f(x3 ) a1 = − (x2 − x3 )x12 + (x3 − x1 )x22 + (x1 − x2 )x32 a2 =

(x2 − x3 )f(x1 ) + (x3 − x1 )f(x2 ) + (x1 − x2 )f(x3 ) (x2 − x3 )x12 + (x3 − x1 )x22 + (x1 − x2 )x32

 x − x1  x − x2   x − x2  x − x3   x − x1  x − x3  f(x1 ) +  f(x2 ) +   f(x3 )   p(x) =  x − x x − x x − x x − x x − x x − x 2  1 3  3  2   1  2 1  2  3 1  3

N

Secara umum, untuk N data rumus interpolasi Lagrange:

p(x) = ∑ l(x, xi )f(xi ) i=1

 x − xj   l(x, xi ) = ∏    j≠i  xi − xj 

Untuk x = xk (k = 1, …, N):

  xi − xj   = 1, (i = k) ∏   xk − xj   j≠i  xi − xj  = l(xk , xi ) = ∏    j≠i  xi − xj  ...  xk − xk  ... = 0, (i ≠ k)   x −x  j    i l(xk , xi ) = δik

p(xk ) = f(xk )

Perlukah memakai semua N data yang ada? Pada bagian sebelum ini interpolasi menggunakan seluruh N data f( xi ) yang tersedia, yang berarti menggunakan polinomial p(x) berorde N-1. Kini, misal N = 4 dan x berada di sekitar x4 , maka diperoleh:

 x − x2  x − x3  x − x4     l(x, x1 ) =   x1 − x2  x1 − x3  x1 − x4 

 x − x1  x − x3  x − x4     l(x, x2 ) =   x2 − x1  x2 − x3  x2 − x4 

 x − x1  x − x2  x − x4     l(x, x3 ) =   x3 − x1  x3 − x2  x3 − x4 

 x − x1  x − x2  x − x3     l(x, x4 ) =   x4 − x1  x4 − x2  x4 − x3 

Dapat dilihat bahwa, l(x, x1 ) < l(x, x2 ) < l(x, x3 ) < l(x, x4 ) . Ini berarti, semakin jauh dari x pengaruh data f( xi ) semakin kecil dalam menentukan nilai p(x). Data yang penting yaitu yang berada di sekitar titik x. Karena itu, cukup data-data di sekitar titik x yang digunakan. Dengan kata lain, untuk interpolasi cukup digunakan polinomial p(x) berorde rendah, contoh berorde 3 (fungsi kubik).

Interpolasi Lagrange Kubik Interpolasi Lagrange Kubik menggunakan polinomial p(x) berorde 3 sebagai fungsi interpolasi:

p(x) = a0 + a1x + a2x2 + a3x3 ≅ f(x) Untuk mencari nilai aj (j = 0, 1, 2, 3) diperlukan 4 data f( xi ) di sekitar x:

f(x0 ), f(x1 ), f(x2 ), f(x3 )

(xi ≤ x ≤ xi+1 ;

x0 = xi-1 , x1 = xi , x2 = xi+1 , x3 = xi+2 )

untuk membentuk sistem persamaan linear:

a0 + a1xj + a2xj2 + a3xj3 = f(xj )

(j = 0, 1, 2, 3)

Langkah pertama dengan begitu, menentukan xj (j = 0, 1, 2, 3) dengan melihat posisi x di antara titik data xi (i = 1, …, N). Diperoleh 3

p(x) = ∑ l(x, xj )f(xj ) j= 0

 x − xk l(x, xj ) = ∏   k ≠ j  xj − xk

   

Catatan: Karena fungsi interpolasi p(x) dicocokkan dengan data f(x0 = xi-1 ), ..., f(x3 = xi+2 ) maka p(x) berlaku hanya untuk daerah xi-1 ≤ x ≤ xi+2 . Untuk daerah x yang lain berlaku fungsi interpolasi p(x) yang lain. Pada batas antara dua daerah yang bersebelahan, masing-masing fungsi interpolasi p(x) dari kedua daerah berbeda itu menunjukan nilai yang sama, karena dalam menentukan p(x) selalu dibuat agar p(x) cocok dengan setiap titik data dalam daerah itu.

pI (x2 ) = f(xi+2 )

II

xi-1

xi+2 I

xi+5 pII (x0 ) = f(xi+2 )

Dengan kata lain, p(x) bersifat kontinyu. Tetapi, tidak begitu dengan turunannya: p’(x) bersifat diskontinyu pada batas dua daerah yang bersebelahan.

Interpolasi Hermite Kubik Dengan menggunakan polinomial p(x) berorde 3 (kubik), interpolasi dilakukan di antara dua titik data yang berurutan, yaitu dalam interval xi ≤ x ≤ xi+1:

p(x) = a0 + a1x + a2x2 + a3x3 ≅ f(x) Jadi, yang pertama dilakukan yaitu, menentukan posisi x di antara titik data xi (i = 1, …, N). Untuk mencari aj (j = 0, 1, 2, 3) diperlukan 4 persamaan, yang ditetapkan sebagai berikut:

p(x1 ) = a0 + a1 x1 + a2x12 + a3x13 = f(x1 ) p(x2 ) = a0 + a1x2 + a2x22 + a3x23 = f(x2 ) 2 3 1

p'(x1 ) = a1 + 2a2x1 + 3a x = f'(x1 )

(xi ≤ x ≤ xi+1 ; x1 = xi , x2 = xi+1 )

p'(x2 ) = a1 + 2a2x2 + 3a3x22 = f'(x2 ) Jadi, pada interpolasi Hermite diperlukan sebagai data bukan saja f(x) namun juga turunannya f’(x).

Diperoleh aj (j = 0, 1, 2, 3) sebagai berikut:

a0

=

a1

=

a2

=

a3

=

 x12 (3x2 − x1 )f(x2 ) + x22 (x2 − 3x1 )f(x1 )    3 − (x x ) 2 1    x f'(x2 ) + x2 f'(x1 )   − x2 x1  1 2 (x2 − x1 )    f(x2 ) − f(x1 )   − 6x2 x1  3  (x2 − x1 )   x1 (2x2 + x1 )f'(x2 ) + x2 (x2 + 2x1 )f'(x1 )   +  2 (x2 − x1 )    f(x2 ) − f(x1 )   3(x2 + x1 ) 3 − (x x ) 2 1    (x + 2x1 )f'(x2 ) + (2x2 + x1 )f'(x1 )   −  2 2 (x2 − x1 )    f(x2 ) − f(x1 )   f'(x2 ) + f'(x1 )   +   − 2 3 2  (x2 − x1 )   (x2 − x1 ) 

Dengan aj (j = 0, 1, 2, 3) yang sudah diperoleh, didapat fungsi interpolasi p(x) sebagai berikut: 2

(

p(x) = ∑ h1 (x, xj )f(xj ) + h2 (x, xj )f'(xj )

)

j=1

 (x − x1 )  x − x2    h1 (x, x1 ) =  1 − 2 (x1 − x2 )  x1 − x2  

2

 (x − x2 )  x − x1    h1 (x, x2 ) =  1 − 2 (x2 − x1 )  x2 − x1    x − x2   h2 (x, x1 ) = (x − x1 )  x1 − x2 

2

2

 x − x1   h2 (x, x2 ) = (x − x2 )  x2 − x1 

2

Pada interpolasi Hermite bukan saja p(x) yang dicocokkan dengan data f(x) namun juga turunannya p’(x) dicocokkan dengan data f’(x). Karena itu, baik p(x) maupun p’(x) bersifat kontinyu. Ini berbeda dari yang ditemui pada interpolasi Lagrange.

Interpolasi Hermite Orde Lebih Tinggi Interpolasi Hermite tidak terbatas hanya menggunakan polinomial p(x) berorde 3 (kubik), namun dapat juga yang berorde lebih tinggi. Untuk itu diperlukan lebih banyak data, bukan hanya f(x) dan f’(x) pada titik xi dan xi+1 . Secara umum fungsi interpolasi Hermite p(x) berupa polinomial berorde (2n - 1) memerlukan n data f(x) dan n data f’(x): n

(

p(x) = ∑ h1 (x, xj )f(xj ) + h2 (x, xj )f'(xj ) j=1

dengan:

(

)

h1 (x, xj ) = 1 − 2(x − xj )l'(xj ) l2 (x, xj ) h2 (x, xj ) = (x − xj )l2 (x, xj )

 x − xk   l(x, xj ) = ∏    k ≠ j  xj − xk  1 l'(xj ) = ∑ k ≠ j (xj − xk )

)

Interpolasi Hermite Kubik tanpa Data f’(x) Interpolasi Hermite memerlukan sebagai data selain f(x) juga f’(x). Pada beberapa kasus bisa saja data f’(x) tidak tersedia, melainkan hanya data f(x). Pada kasus ini sebenarnya interpolasi Hermite tidak bisa dipakai. Tetapi, jika f’(x) bisa diperoleh melalui pendekatan (approximation) maka, interpolasi Hermite bisa dipakai. f’( xi ) dapat dihitung sebagai turunan sebuah fungsi kuadratik g(x), yang dicocokkan dengan data f(x) pada titik-titik xi-1 , xi , xi+1 :

g(x) = ax 2 + bx + c

g(xi-1 ) = f(xi-1 )   g(xi ) = f(xi )  (a, b, c) g(xi+1 ) = f(xi+1 ) 

f'(xi ) ≅ g'(xi ) = 2axi + b

Dapat dilihat bahwa, proses pencarian f’(x) ini berdiri sendiri, berada di luar atau bukan bagian dari proses interpolasi Hermite. Dengan begitu, sifat kontinyu fungsi interpolasi Hermite p(x) dan turunannya p’(x) tidak berubah.

Dari sistem persamaan linear:

axi2-1 + bxi-1 + c = f(xi-1 ) axi2 + bxi + c = f(xi ) axi2+1 + bxi+1 + c = f(xi+1 ) diperoleh:

a=

f(xi−1 ) f(xi ) f(xi+1 ) + + (xi−1 − xi )(xi−1 − xi+1 ) (xi − xi−1 )(xi − xi+1 ) (xi+1 − xi−1 )(xi+1 − xi )

b=−

(xi + xi+1 )f(xi−1 ) (xi-1 + xi+1 )f(xi ) (xi-1 + xi )f(xi+1 ) − − (xi−1 − xi )(xi−1 − xi+1 ) (xi − xi−1 )(xi − xi+1 ) (xi+1 − xi−1 )(xi+1 − xi )

sehingga:

 xi − xi+1  f(xi−1 )  1  xi − xi-1  f(xi+1 ) 1       f'(xi ) ≅  + + f(xi ) +     xi−1 − xi+1  (xi−1 − xi )  xi − xi+1 xi − xi−1   xi+1 − xi−1  (xi+1 − xi )

Jika diaplikasikan pada interpolasi Hermite kubik: 2

(

p(x) = ∑ h1 (x, xj )f(xj ) + h2 (x, xj )f'(xj )

)

j=1

maka diperoleh fungsi interpolasi Hermite kubik p(x) sebagai berikut:

3

p(x) = ∑ h(x, xj )f(xj )

(xi ≤ x ≤ xi+1 ; x0 = xi-1 , x1 = xi , x2 = xi+1 , x3 = xi+2 )

j= 0

 x1 − x2  1  h(x, x0 ) = h2 (x, x1 )  x0 − x2  (x0 − x1 )  1  x − x3  1  1  + h2 (x, x2 ) 2  h(x, x1 ) = h1 (x, x1 ) + h2 (x, x1 ) +  x1 − x2 x1 − x0   x1 − x3  (x1 − x2 )  1  x − x0  1  1  + h2 (x, x1 ) 1  h(x, x2 ) = h1 (x, x2 ) + h2 (x, x2 ) +  x2 − x1 x2 − x3   x2 − x0  (x2 − x1 )  x − x1  1  h(x, x3 ) = h2 (x, x2 ) 2  x3 − x1  (x3 − x2 )

Interpolasi Spline Kubik Seperti interpolasi Lagrange, interpolasi Spline kubik juga memerlukan hanya f(x) sebagai data. Namun, turunan fungsi interpolasi Spline kubik p’(x) dibuat bersifat kontinyu. Interpolasi Spline kubik menggunakan polinomial p(x) orde 3, untuk xi ≤ x ≤ xi+1 :

p(x) = di + ci (x − xi ) + bi (x − xi )2 + ai (x − xi )3 ≅ f(x) Turunan pertama dan kedua p(x) yaitu:

p'(x) = ci + 2bi (x − xi ) + 3ai (x − xi )2 p''(x) = 2bi + 6ai (x − xi ) Evaluasi pada titik x = xi menghasilkan:

pi ≡ p(xi ) = di = f(xi )

p''i ≡ p''(xi ) = 2bi

dan pada titik x = xi+1:

p''i+1 ≡ p''(xi+1 ) = 2bi + 6aihi pi+1 ≡ p(xi+1 ) = di + cihi + bihi2 + aihi3 = f(xi+1 ) hi ≡ xi+1 − xi

Jadi,

di = pi

bi =

p''i 2

ai =

p''i+1 −p''i 6hi

ci =

pi+1 − pi hip''i+1 +2hip''i − hi 6

sehingga diperoleh:

 p − p h p''  p'' −p''i  h p''  p'' (x − xi )3 p(x) = pi +  i+1 i − i i+1 − i i (x − xi ) + i (x − xi )2 +  i+1 6 3  2  hi  6hi 

p'(x) =

 p'' −p''i  pi+1 − pi hip''i+1 hip''i (x − xi )2 − − + p''i (x − xi ) +  i+1 hi 6 3  2hi 

p(x) telah dicocokkan dengan data f(x) di titik-titik batas interval, sehingga bersifat kontinyu. Untuk membuat p’(x) kontinyu maka dicari ekspresi p’(x) untuk daerah sebelumnya xi-1 ≤ x ≤ xi :

p'(x) =

 p'' −p''i-1  pi − pi-1 hi-1p''i hi-1p''i-1 (x − xi-1 )2 − − + p''i-1 (x − xi-1 ) +  i hi-1 6 3  2hi-1 

dan disamakan dengan p’(x) untuk daerah xi ≤ x ≤ xi+1 di titik x = xi .

Untuk N = jumlah data, diperoleh:

 pi+1 − pi pi − pi-1   hi-1p''i-1 +2(hi-1 + hi )p''i +hip''i+1 = 6 − hi-1   hi

(i = 2, ..., N - 1)

Untuk menghitung p(x) diperlukan p’’(x) di semua N titik data. (N-2) buah persamaan di atas tidak cukup untuk mendapatkan p’’(x) di semua titik data. Masih diperlukan 2 persamaan lagi, yang diperoleh dengan mengevaluasi p’(x) di titik awal x = x1 (memakai ekspresi p’(x) untuk x1 ≤ x ≤ x2 ) dan akhir x = xN (memakai ekspresi p’(x) untuk xN-1 ≤ x ≤ xN ). Didapat:

(i = 1) (i = N)

p −p  2h1p''1 +h1p''2 = 6 2 1 − p'1   h1   p −p  hN-1p''N-1 +2hN-1p''N = 6 p'N − N N-1  hN-1  

Masalah: p’(x) di titik awal x = x1 dan akhir x = xN tidak diketahui,

??

Ada dua cara. Pertama yang disebut spline alamiah yaitu, menetapkan p’’(x) di titik awal x = x1 dan akhir x = xN sama dengan nol. Kedua, menebak nilai p’(x) di titik awal x = x1 dan akhir x = xN .

Interpolasi Multidimensi Jika data bergantung pada lebih dari satu variabel, maka dilakukan interpolasi multidimensi. Metode interpolasi yang telah disampaikan bisa dipakai untuk melakukan interpolasi multidimensi. Sebagai contoh di sini ditunjukkan interpolasi 2 dimensi. Untuk dimensi lebih tinggi berlaku cara yang sama.

n

m

i=1

j=1

p(x, y) = ∑ S(x, xi )∑ S(y, yj )f(xi , yj )

Pada contoh di atas, interpolasi menggunakan (n x m) data f(x,y). Interpolasi dilakukan per dimensi: Untuk satu titik data x tertentu dilakukan interpolasi di sepanjang sumbu y, hal yang sama dilakukan untuk semua titik data x yang lain. Prinsip yang sama berlaku untuk interpolasi berdimensi lebih tinggi.

Contoh, interpolasi Lagrange kubik:

3

3

i= 0

j= 0

p(x, y) = ∑ l(x, xi )∑ l(y, yj )f(xi , yj )  x − xk   l(x, xi ) = ∏  k ≠i  xi − xk   y − ys   l(y, yj ) = ∏    s ≠ j  yj − ys 

Kembali ke contoh problem least square: Kuat medan listrik E di sekitar sebuah benda berbentuk lempeng diukur pada jarak 10 cm dari pusat massanya dan arah yang bervariasi. Arah dinyatakan dalam sudut θ terhadap sumbu y yang ditetapkan sebelum pengukuran. Diperoleh data sebagai berikut:

θ [derajat]

E [V/cm]

10 15 20 25 30 35 40 45 50 55

0.01794775 0.03808997 0.05516225 0.05598281 0.04795629 0.04807485 0.06273566 0.07853982 0.07395442 0.04201338

y

θ

Dengan interpolasi, cari nilai p(x) di sepanjang titik data.

E

Integrasi Menghitung luas daerah di bawah kurva: f(x)

f(x) analitik

numerik

b

∑ w f(x )

∫ f(x) dx

i

i

a

a

b b

N

a

i=1

I = ∫ f(x) dx ≅ ∑ wif(xi )

i

x

a

b

x

Integral numerik sering disebut juga sebagai quadrature; integrasi numerik disebut sebagai integrasi dgn menjumlah quadrature.

Meski tidak terlihat pada rumus akhir, pada integrasi numerik integrand f(x) diinterpolasi dengan suatu polinomial: b

N

a

i=1

I = ∫ f(x) dx ≅ ∑ wif(xi )

f(x) ≅ p(x)

polinomial

Dilihat dari titik-titik xi tempat integrand f(x) dihitung, ada teknik integrasi numerik yang menggunakan xi berjarak tetap dan ada yang memakai xi berjarak tidak tetap. Contoh (akan dibahas): •

quadrature trapezoid, Simpson menggunakan xi berjarak sama,

•

quadrature Gaussian menggunakan xi berjarak tidak sama.

Quadrature Trapezoid Kurva integrand f(x) diinterpolasi dengan sebuah garis lurus (f(x) diinterpolasi dengan fungsi linier / polinomial orde 1): b

b

N

I = ∫ f(x) dx ≅ ∫ p(x) dx = ∑ wip(xi ), a

a

p(x) = r + sx

i=1

f(x) Untuk menarik garis lurus diperlukan minimal 2 titik, dipilih titik f(a) dan f(b):

p(x)

p(a) = f(a), p(b) = f(b) b

∫ p(x) dx a

a

b

x

Dengan diketahui hanya p(a) dan p(b) (r dan s tidak dicari), maka integrasi numerik dikerjakan untuk N = 2: b

2

∫ p(x) dx = ∑ w p(x ) = w p(x ) + w p(x ) = w p(a) + w p(b) i

a

1

1

2

2

1

w1 , w2 = ?

2

i=1

Mencari w1 dan w2 :

p(x) = r + sx

i

b

∫ (r + sx) dx = w (r + sa) + w (r + sb) 1

2

a

1 r(b - a) + s(b2 − a2 ) = r(w1 + w2 ) + s(aw1 + bw2 ) 2 w1 + w2 = b - a 1 aw1 + bw2 = (b2 − a2 ) 2 b

Rumus quadrature trapezoid:

I = ∫ f(x) dx ≅ a

h (f(a) + f(b) ) 2

w1 = w2 =

1 (b − a) 2

(h = b − a)

luas trapezoid (lihat gambar)

Quadrature Simpson & Boole Cara yang sama seperti pada quadrature trapezoid bisa dipakai untuk polinomial p(x) orde lebih tinggi. Contoh, quadrature Simpson memakai p(x) fungsi kuadratik / polinomial orde 2 untuk menginterpolasi integrand f(x): c

c

N

I = ∫ f(x) dx ≅ ∫ p(x) dx = ∑ wip(xi ), a

a

p(x) = r + sx + tx2

i=1

f(x) Untuk membuat kurva kuadratik diperlukan minimal 3 titik, dipilih titik f(a), f(b) dan f(c):

p(x)

p(a) = f(a), p(b) = f(b),

b

∫ p(x) dx

p(c) = f(c)

a

dengan

b=

a+c 2

a

b

c

x

Integrasi numerik dikerjakan untuk N = 3: c

3

∫ p(x) dx = ∑ w p(x ) = w p(a) + w p(b) + w p(c) i

a

i

1

2

w1 , w2 , w3 = ?

3

i=1

Mencari w1 , w2 , w3:

p(x) = r + sx + tx2 c

2 2 2 (r + sx + tx ) dx = w (r + sa + ta ) + w (r + sb + tb ) 1 2 ∫ a

r(c - a) +

1 1 s(c2 − a2 ) + t(c3 − a3 ) = 2 3

w1 + w2 + w3 = c - a 1 aw1 + bw2 + cw3 = (c2 − a2 ) 2 1 a2w1 + b2w2 + c2w3 = (c3 − a3 ) 3

+ w3 (r + sc + tc2 ) r(w1 + w2 + w3 ) + s(aw1 + bw2 + cw3 ) + t(a2w1 + b2w2 + c2w3 )

w1 = w3 = w2 =

1 (c − a) 6

2 (c − a) 3

c

Diperoleh Rumus quadrature Simpson:

I = ∫ f(x) dx ≅ a

dengan h =

h (f(a) + 4f(b) + f(c) ) 3

c−a yaitu jarak antar titik xi tempat f(x) dihitung: h = b − a = c − b 2

Dengan cara yang sama, menggunakan p(x) polinomial orde 3 diperoleh rumus quadrature Simpson 3 8 : d

I = ∫ f(x) dx ≅ a

3h (f(a) + 3f(b) + 3f(c) + f(d) ) 8

d-a   = b − a = c − b = d − c h = 3  

dan dengan p(x) polinomial orde 4 rumus quadrature Boole: e

2h (7f(a) + 32f(b) + 12f(c) + 32f(d) + 7f(e) ) I = ∫ f(x) dx ≅ 45 a

e-a  h =  4     

 = b− a  = c −b = d−c  = e − d

Integrasi Komposit Polinomial orde rendah memadai untuk menginterpolasi sebuah fungsi dalam daerah yang sempit. Untuk daerah yang lebar diperlukan orde yang lebih tinggi. Alternatif lain yaitu, membagi daerah fungsi yang lebar itu dalam beberapa daerah yang sempit, lalu di tiap daerah yang sempit itu digunakan polinomial orde rendah untuk interpolasi. Quadrature trapezoid dan Simpson pada dasarnya memadai untuk daerah integrasi yang sempit, namun dengan membagi daerah integrasi dalam beberapa daerah yang sempit, maka quadrature trapezoid dan Simpson bisa dipakai juga untuk daerah integrasi yang lebar. Integral total merupakan jumlah semua integral untuk daerah yang sempit. Integrasi seperti ini disebut integrasi komposit. Bergantung pada integrand f(x), daerah integrasi yang lebar bisa dibagi dalam beberapa daerah sempit yang sama atau berbeda panjang. Juga, semua integral untuk daerah yang sempit bisa dihitung menurut rumus quadrature yang sama, misal semuanya trapezoid, atau berbeda-beda, sesuai kurva di tiap daerah sempit itu. Kasus sederhana yaitu, bila daerah integrasi dibagi sama panjang dan untuk tiap daerah digunakan rumus quadrature yang sama.

Contoh, daerah integrasi [a,b] dibagi dalam N bagian sama panjang. b

a +d

a +2d

b -d

b

a

a

a +d

b-2d

b -d

I = ∫ f(x) dx =

•

∫ f(x) dx +

∫ f(x) dx + ... +

∫ f(x) dx + ∫ f(x) dx

b−a  d =   N  

integrasi komposit menggunakan quadrature trapezoid b

I = ∫ f(x) dx ≅ h[21 (f0 + fN ) + f1 + f2 + ... + fN −1 ] a

h=

•

b−a , fi = f(a + ih), i = 0, ..., N N

integrasi komposit menggunakan quadrature Simpson b

I = ∫ f(x) dx ≅ a

2h 1 [2 (f0 + f2N ) + 2(f1 + f3 + ... + f2N−1 ) + f2 + f4 + ... + f2N−2 ] 3 h=

b−a , fi = f(a + ih), i = 0, ..., 2N 2N

Integrasi komposit trapezoid untuk daerah integrasi [a,b] yang dibagi 8 sama panjang: b

I = ∫ f(x) dx ≅ h[21 (f0 + f8 ) + f1 + f2 + f3 + f4 + f5 + f6 + f7 ] a

f(x)

a

h

b

x

Integrasi komposit yang menggunakan quadrature trapezoid dan Simpson; daerah integrasi [a,b] yang dibagi 3: b

I = ∫ f(x) dx ≅ a

h1 (fa + 2fa +h1 + fc ) + h2 (fc + 4fc+h2 + fb ) 2 3

f(x) Simpson

trapezoid

a

b

c h1

h1

2h2

x

Quadrature Gaussian Quadrature Gaussian memanfaatkan polinomial yang memiliki sifat orthogonal dan ternormalisasi sebagai berikut: b

∫ v(x)O (x)O n

m

n

(x) dx = δnm ,

On (x) = ∑ bixi i= 0

a

Contoh: 1

On =

2n +1 2 n

P , Pn = polinomial Legendre

∫ O (x)O n

m

(x) dx = δnm

-1

On = n!1 Ln , Ln =

∞

polinomial Laguerre

∫e

-x

On (x)Om (x) dx = δnm

0

Dengan quadrature Gaussian, dievaluasi integral berbentuk: b

N

∫ v(x)f(x)dx = ∑ w f(x ) i

a

i=1

i

wi , xi = ?

Mencari xi : Anggap integrand f(x) merupakan polinomial orde 2N-1 (atau katakan saja f(x) diinterpolasi dengan polinomial p(x) orde 2N-1): 2N -1

f(x) ≅ p(x) =

N -1

∑ aix = r(x) + s(x) i

r(x) = ∑ aix , s(x) =

dengan

i= 0

i

i= 0

2N -1

i a x ∑ i i =N

s(x) bisa ditulis sebagai s(x) = q(x)ON (x) dengan q(x) polinomial orde N-1: N -1

N -1

q(x) = ∑ dix = ∑ ciOi (x) i

i= 0 b

Maka:

Secara numerik:

i= 0 b

N -1

∫ v(x)s(x)dx = ∑ c ∫ v(x)O (x)O i

a

i= 0

b

N

i

N

N -1

(x)dx = ∑ ciδiN = 0 i= 0

a N

∫ v(x)s(x)dx = ∑ w s(x ) = ∑ w q(x )O i

a

i=1

i

i

i

N

(xi ) = 0

i=1

Mengingat q(x) fungsi sembarang, persamaan terakhir dipenuhi hanya jika ON (xi ) = 0 (i = 1, ..., N) . xi = akar polinomial ON (x) (i = 1, ..., N)

Mencari wi : Untuk integrand f(x) dan s(x), yang merupakan polinomial orde 2N-1 berlaku integrasi numerik: b

N

b

N

∫ v(x)f(x)dx = ∑ w f(x )

∫ v(x)s(x)dx = ∑ w s(x )

a

a

i

i

i=1

i

i

i=1

Integrasi numerik yang sama tentu berlaku juga untuk integrand polinomial orde lebih rendah, contohnya r(x), yang berorde N-1: b

N

a

i=1

∫ v(x)r(x)dx = ∑ wir(xi ),

N-1

r(x) = ∑ aixi i= 0

Dari penurunan rumus quadrature trapezoid, Simpson dll sebelum ini diketahui bahwa untuk mencari wi bisa digunakan r(x) sembarang polinomial orde N-1 (koefisien ai tidak diperlukan). Karena itu, dipilih r(x) yang memudahkan:

 x − xj  , (i, j = 1, ..., N) r(x) = l(x, xi ) = ∏    j≠i  xi − xj  b

N

a

i=1

Diperoleh: ∫ v(x)l(x, xj )dx = ∑ wil(xi , xj ) = wj

l(xk , xi ) = δik b

wj = ∫ v(x)l(x, xj )dx (j = 1, ..., N) a

Pada integrasi numerik Gaussian, diperlukan N buah titik evaluasi xi untuk integrand f(x) ≅ p(x) polinomial orde 2N-1. Pada integrasi numerik seperti quadrature trapezoid dan Simpson, diperlukan 2N buah titik xi untuk integrand f(x) ≅ p(x) polinomial orde 2N-1: trapezoid : 2N = 2 Simpson

: 2N = 3

Simpson 3 8 : 2N = 4 Boole

: 2N = 5

dst Secara umum, dengan begitu, quadrature Gaussian memerlukan titik evaluasi lebih sedikit (separuh) dari yang diperlukan integrasi numerik yang mengikuti cara seperti quadrature trapezoid dan Simpson. Bergantung pada keperluan, integrasi komposit juga bisa diterapkan menggunakan quadrature Gaussian atau campuran quadrature Gaussian dan yang lain.

Quadrature Gauss-Legendre Quadrature Gauss-Legendre menggunakan polinomial Legendre Pn : 1

On =

2n +1 2 n

P:

∫ O (x)O n

m

(x) dx = δnm

-1

Asalnya, quadrature Gauss-Legendre dipakai untuk integral berbatas [-1,1]: 1

N

∫ f(x)dx = ∑ w f(x ) i

i

i=1

-1

Namun dengan mengganti variabel integrasi, quadrature Gauss-Legendre dapat juga dipakai untuk mengevaluasi integral dengan batas bukan [-1,1]. 1

Contoh:

N

b

2 f(x)dx = w f(x ) f(y)dy = ∑ i i ∫-1 ∫ b−a a i=1

y − a x − ( −1) x + 1 = = b − a 1 − ( −1) 2 (transformasi linier)

b

N

∫ f(y)dy = ∑ u f(y ) i

a

i=1 1 2

i

yi = (xi + 1)(b − a) + a b−a  ui =  wi  2 

Contoh xi dan wi quadrature Gauss-Legendre untuk beberapa N terkecil: 1

N

-1

i=1

∫ f(x)dx = ∑ w f(x ) i

N

x

i

w

2

± 0.577350269189626

1.000000000000000

3

± 0.774596669241483

0.555555555555556

0.000000000000000

0.888888888888889

4

5

± 0.861136311594053

0.347854845137454

± 0.339981043584856

0.652145154862546

± 0.906179845938664

0.236926885056189

± 0.538469310105683

0.478628670499367

0.000000000000000

0.568888888888889

Distribusi xi pada quadrature Gauss-Legendre tidak merata seperti distribusi pada quadrature trapezoid dan Simpson. Makin dekat ke batas-batas integral distribusi makin rapat. Distribusi itu simetris terhadap garis x = 0. Ilustrasi untuk N = 11: x = -1

x=0

x=1

x

Distribusi ini lebih cocok untuk integrand f(x) yang bentuk kurvanya lebih tajam di sekitar batas integral, sementara kurang tajam di bagian tengah.

Untuk f(x) yang berkurva tajam di bagian tengah dan kurang tajam di sekitar batas integral diperlukan beberapa penanganan (mis. membagi daerah integrasi, redistribusi x dll).

Quadrature Gauss-Laguerre Quadrature Gauss-Laguerre menggunakan polinomial Laguerre Ln : ∞

On = n!1 Ln :

-x e ∫ On (x)Om (x) dx = δnm 0 ∞

dipakai untuk integral berbentuk:

∫e 0

-x

N

f(x)dx = ∑ wif(xi ) i=1

Contoh xi dan wi quadrature Gauss-Laguerre untuk beberapa N: N

x

w

2

0.585786437626905

0.853553390593274

3.414213562373095

0.146446609406726

0.322547689619392

0.603154104341634

1.745761101158347

0.357418692437800

4.536620296921128

0.038887908515005

9.395070912301133

0.000539294705561

4

Lain-Lain Mengganti Variabel Integrasi Pada topik quadrature Gauss-Legendre terdapat contoh penggantian variabel integrasi. Penggantian variabel integrasi bisa juga diperlukan pada kasus lain. Tujuannya, agar evaluasi integral menjadi lebih mudah dan hasilnya baik. ∞

Contoh:

dx I=∫ 1 + x2 0 transformasi:

batas integral sampai tak behingga, jika dievaluasi langsung memerlukan sangat banyak titik, tidak praktis dan hasilnya bisa saja buruk

x=

1+ y 2 , dx = dy 2 1− y (1 − y) 1

I=∫ -1

2dy 2

(

(1 − y) 1 + 1

( )) 1+ y 2 1− y

dy (1 − y)2 + (1 + y)2 -1

= 2∫

quadrature Gauss-Legendre

Meringkas Daerah Integrasi Beberapa fungsi bersifat genap, ini memungkinkan daerah integrasi diringkas menjadi separuhnya (mengurangi jumlah titik evaluasi 2N menjadi N). fungsi genap: f( −x) = f(x)

Contoh:

fungsi ganjil: f( −x) = −f(x)

a

2N dx wi • I=∫ = ∑ 2 2 1 + x i=1 1 + xi -a a

2N dx wi = 2∫ = 2 ∑ 2 2 1 + x 1 + x i = N + 1 i 0 1

2N dy wi • I=∫ = ∑ 2 2 2 2 (1 − y) + (1 + y) i=1 (1 − yi ) + (1 + yi ) -1 1

2N dy wi = 2∫ = 2 ∑ 2 2 2 2 (1 − y) + (1 + y) (1 − y ) + (1 + y ) i N 1 = + i i 0

Beberapa fungsi memiliki simetri, contoh fungsi trigonometri:

sin(π ± x) = ∓ sin(x) cos(π ± x) = −cos(x)

sin( −x) = −sin(x) cos( −x) = cos(x)

Dengan memanfaatkan relasi simetri di atas batas integrasi sebuah integral tertutup (loop) seperti contoh di bawah dapat diringkas menjadi seperempatnya, sehingga jumlah titik evaluasi berkurang banyak: 2π

I=

integral tertutup bisa dimulai dari titik mana saja

im(x − a) [ ] f(sin(x − a)) + f(cos(x − a)) e dx ∫ 0

2π

=

∫ [f(sin(x)) + f(cos(x))]e

imx

dx

0 π

[

= ∫ {f(sin(x)) + f(cos(x))}e

imx

+ {f( −sin(x)) + f( −cos(x))}e

im(x + π)

]dx

0 π

=

2

∫ [f(sin(x))(e

imx

)

(

+ eim(π −x) + f(cos(x)) eimx + eim(2π −x)

0

(

)

(

telah dipakai

x = x'+ π

)

x = −x'

)]

+ f( −sin(x)) eim(π + x) + eim(2π −x) + f( −cos(x)) eim(π + x) + eim(π −x) dx

Menangani Singularitas

Kadang ditemui integrand f(x) yang memiliki singularitas dalam daerah integrasi. Salah satu cara menangani singularitas yaitu subtraksi, yang dimulai dengan menambahkan integral bernilai nol pada integral yang dihitung. a

Contoh:

dx 0 (1 + x) x

• I=∫ a

singular pada x = 0 a

a

dx dx dx =∫ −∫ +∫ x 0 x 0 (1 + x) x 0

ditambah nol

a  1 1  dx dx + ∫ = ∫  − x x 0  (1 + x) x 0

subtraksi pada integral asal

a

a

x dx + 2 a 1 + x 0

= −∫

a

x2f(x)dx • I=∫ 2 (b − x2 ) 0

singular pada x = b

(0 < b2 ≤ a2 )

a

∞

x2f(x)dx b2f(b)dx =∫ 2 −∫ 2 2 2 (b x ) (b x ) − − 0 0 a

=∫

(x f(x) − b f(b))dx − 2

2

(b2 − x2 )

0 a

=∫

ditambah nol (lihat *)

∞

b2 f(b)dx ∫a (b2 − x2 )

subtraksi pada integral asal

(x f(x) − b f(b))dx − 1 bf(b)ln a − b  2

2

(b2 − x2 )

0

∞

(*)

2

∞

 a +b 

∞

∞

dx 1  1 1  1 dx 1 dx = + dx = = ∫0 (b2 − x2 ) 2b ∫0  b − x b + x  2b −∫∞b + x 2b −∫∞ x = 0

Quadrature Filon Bisa saja ditemui integrand f(x) yang sangat berosilasi; dalam jarak yang pendek f(x) berubah-ubah naik turun. Dengan macam-macam quadrature yang sudah disampaikan, integrasi menjadi sulit karena dibutuhkan banyak sekali titik evaluasi. Integral seperti ini dapat dihitung dengan menggunakan rumus quadrature Filon (M. Abramowitz & I. A. Stegun, Handbook of Mathematical Function, Dover Publications, Inc., NY, 1972). f(x)

x

Quadrature Filon (tanpa suku kesalahan, yang bisa diabaikan): b

∫ f(x)cos(tx)dx = h[α(th)(f

2n

a

sin(tb) − f0 sin(ta)) + β(th)Cgenap + γ(th)Cganjil

n

Cgenap = ∑ f2icos(tx2i ) − 21 (f2n cos(tb) + f0 cos(ta) ) i= 0 n

Cganjil = ∑ f2i−1cos(tx2i−1 )

]

b−a = xi+1 − xi 2n x0 = a

h=

fj = f(xj )

i=1 b

∫ f(x)sin(tx)dx = h[α(th)(f cos(ta) − f 0

a

2n

cos(tb)) + β(th)Sgenap + γ(th)Sganjil

]

n

Sgenap = ∑ f2isin(tx2i ) − 21 (f2n sin(tb) + f0 sin(ta) ) i= 0 n

Sganjil = ∑ f2i−1 sin(tx2i−1 ) i=1

1 sin(2x) 2sin2x − α(x) = + x 2x2 x3  1 + cos2x sin(2x)   − β(x) = 2 2 3 x x    sinx cosx  γ(x) = 4 3 − 2  x   x

untuk nilai x kecil:

2x3 2x 5 2x 7 α(x) = − + − ... 45 315 4725 2 2x2 4x 4 2x 6 β(x) = + − + − ... 3 15 105 567 4 2x2 x 4 x6 γ(x) = − + − + ... 3 15 210 11340

Integrasi Monte Carlo Mungkin saja cara-cara integrasi numerik yang sudah disampaikan sulit atau tidak bisa diterapkan untuk mengevaluasi suatu integral. Pada keadaan ini, integrasi Monte Carlo dapat dipilih. Integrasi Monte Carlo tidak menggunakan interpolasi seperti pada cara-cara integrasi numerik sebelum ini. Integral dianggap sebagai satu persegi panjang, dengan lebar daerah integrasi dan tinggi nilai rata-rata integrand f(x), yang diperoleh melalui statistik dengan memanfaatkan bilangan acak: f(x)

1 n < f(x) >= ∑ f(xi ) n i=1 xi = bilangan acak : a ≤ xi ≤ b b

1 n I = ∫ f(x)dx ≅ (b - a) ∑ f(xi ) n i=1 a

(b-a) x a

b

Persamaan Differensial Persamaan differensial (PD) yang dibahas meliputi persamaan differensial biasa dan serta persamaan differensial parsial. Beberapa persamaan differensial merupakan juga persamaan eigenvalue, contoh persamaan untuk senar gitar (gelombang berdiri). Karena itu, akan dibahas juga persamaan eigenvalue.

Persamaan Differensial Biasa

Pada bagian ini disampaikan metode numerik untuk menyelesaikan persamaan differensial biasa orde 1 dan 2. Dua masalah yang akan dibahas yaitu: •

PD dengan syarat awal

•

PD dengan syarat batas

PD dengan Syarat Awal PD Orde 1 dy = f(x, y) dx

Bentuk umum PD orde 1:

y'=

Diketahui:

y(x0 ) = y0 y

Integrasi:

y(x) = ?

x

Masalah persamaan differensial berubah menjadi masalah persamaan integral.

∫ dy = ∫ f(x, y)dx y0

x0

x

y(x) = y0 + ∫ f(x, y)dx x0

x0 +h

Dicari y(x) pada titik x = x0 + h :

y(x0 + h) = y0 +

∫ f(x, y)dx

x0

Setelah y(x0 + h) didapat, selanjutnya dicari y(x0 + 2h) . Demikian seterusnya.

Metode Euler Menurut metode Euler: f(x,y)

f(x0 , y0 )

f(x,y) dianggap konstan dan dihitung pada x = x0.

x

x0 x0 + h Diperoleh:

y(x0 + h) yg diperoleh

y(x) x0 +h

y(x0 + h) ≅ y0 + f(x0 , y0 )

∫ dx

x0

≅ y0 + hf(x0 , y0 )

y(x0 + h) sebenarnya

y0 x

x0 x0 + h

Metode Euler yang Dimodifikasi f(x,y)

f(x0 + 21 h, y(x0 + 21 h))

Modifikasi dilakukan dalam memilih nilai f(x,y) yang dianggap konstan. Dipilih f(x,y) pada titik x = x0 + 21 h :

x

f(x0 + h, y(x0 + h)) 1 2

1 2

x0 x0 + h x0 + 21 h

dengan y(x0 + 21 h) dihitung memakai metode Euler:

y(x0 + 21 h) ≅ y0 + 21 hf(x0 , y0 )

y(x0 + h) yg diperoleh

y(x)

Diperoleh:

y(x0 + h) ≅ y0 + hf(x0 + 21 h, y(x0 + 21 h)) ≅ y0 + hf(x0 + 21 h, y0 + 21 hf(x0 , y0 ))

y0 x

x0 x0 + h x0 + 21 h

y(x0 + h) sebenarnya

Metode Euler yang Lebih Baik (Improved) Kali ini dipakai nilai f(x,y) yang merupakan rata-rata dari dua nilai f(x,y), masing-masing pada titik x0 dan x0 + h : 1 2

[f(x0 , y0 ) + f(x0 + h, y(x0 + h))]

f(x,y)

Ini sama dengan menggunakan quadrature trapezoid untuk mengevaluasi integral: x0 +h

∫ f(x, y)dx ≅ h[f(x , y 1 2

0

0

) + f(x0 + h, y(x0 + h))]

x

x0 x0 + h

x0

dengan y(x0 + h) dihitung memakai metode Euler:

y(x0 + h) ≅ y0 + hf(x0 , y0 ) Diperoleh:

y(x0 + h) ≅ y0 + 21 h[f(x0 , y0 ) + f(x0 + h, y(x0 + h))] ≅ y0 + 21 h[f(x0 , y0 ) + f(x0 + h, y0 + hf(x0 , y0 ))]

Metode Runge-Kutta Metode Euler dan variasinya sebelum ini sebetulnya termasuk metode RungeKutta, yang menyatakan solusi PD y(x) dalam turunannya f(x,y), yang dihitung untuk argumen x,y yang bervariasi. Sebuah metode Runge-Kutta disebut berorde n jika memiliki suku koreksi O(hn+1 ) (diperoleh dari ekspansi Taylor):

y(x0 + h) = ydiperoleh + O(hn +1 ) Menurut hal itu, metode Euler merupakan metode Runge-Kutta orde 1 sedangkan metode Euler yang dimodifikasi dan yang lebih baik (improved) merupakan metode Runge-Kutta orde 2: Euler

2 : y(x0 + h) = y0 + hf(x0 , y0 ) + O(h )

Euler yg dimodifikasi: y(x0 + h) = y0 + hf(x0 + 21 h, y0 + 21 hf(x0 , y0 )) + O(h3 ) Euler yg lebih baik

: y(x0 + h) = y0 + 21 h[f(x0 , y0 ) + f(x0 + h, y0 + hf(x0 , y0 ))] + O(h3 )

Metode Runge-Kutta yang paling banyak digunakan orang yaitu berorde 4, yang sering diingat sebagai metode Runge-Kutta tanpa tambahan keterangan ‘orde 4’.

Untuk mendapatkan rumus metode Runge-Kutta orde 4, orang bisa memulai dengan mengevaluasi integral f(x,y) memakai quadrature Simpson: x0 +h

1 1 1 [ f(x, y)dx ≅ h f(x , y ) + 4f(x + h, y(x + h)) + f(x0 + h, y(x0 + h))] 0 0 0 0 6 2 2 ∫

x0

≅ 61 h(f0 + 2f1 + 2f2 + f3 ) dengan:

f0 = f(x0 , y0 ) f1 = f(x0 + 21 h, y(x0 + 21 h))

f2 = f(x0 + 21 h, y(x0 + 21 h)) f3 = f(x0 + h, y(x0 + h))

f1 dan f2 memiliki nilai berbeda, karena dihitung untuk nilai argumen y(x0 + 21 h) yang berbeda: menurut metode Euler, y(x0 + 21 h) dapat diperoleh melalui 2 persamaan: (1)

y(x0 + 21 h) ≅ y0 + 21 hf0

f1 = f(x0 + 21 h, y0 + 21 hf0 )

atau (2) y0 ≅ y(x0 + 21 h) − 21 hf(x0 + 21 h, y(x0 + 21 h)) ≅ y(x0 + 21 h) − 21 hf(x0 + 21 h, y0 + 21 hf0 ) y(x0 + 21 h) ≅ y0 + 21 hf1

f2 = f(x0 + 21 h, y0 + 21 hf1 )

Untuk f3 , digunakan metode Euler yang dimodifikasi untuk mencari y(x0 + h) :

y(x0 + h) ≅ y0 + hf(x0 + 21 h, y(x0 + 21 h))

≅ y0 + hf(x0 + 21 h, y0 + 21 hf1 ) ≅ y0 + hf2

f3 = f(x0 + h, y0 + hf2 )

Jadi, menurut metode Runge-Kutta orde 4:

y(x0 + h) = y0 + 61 h(f0 + 2f1 + 2f2 + f3 ) dengan:

f0 = f(x0 , y0 )

f2 = f(x0 + 21 h, y0 + 21 hf1 )

f1 = f(x0 + 21 h, y0 + 21 hf0 )

f3 = f(x0 + h, y0 + hf2 )

Berangkat dengan quadrature Simpson, orang juga bisa memperoleh rumus metode Runge-Kutta orde 3: x0 +h

∫ f(x, y)dx ≅ h(f + 4f + f ) 1 6

0

1

2

x0

dengan:

f0 = f(x0 , y0 )

f1 = f(x0 + 21 h, y(x0 + 21 h))

f2 = f(x0 + h, y(x0 + h))

y(x0 + 21 h) dicari dengan metode Euler dan y(x0 + h) dengan metode Euler yang dimodifikasi: y(x0 + 21 h) ≅ y0 + 21 hf0 y(x0 + h) ≅ y0 + hf1

f1 = f(x0 + 21 h, y0 + 21 hf0 ) f2 = f(x0 + h, y0 + hf1 )

Jadi, menurut metode Runge-Kutta orde 3:

y(x0 + h) = y0 + 61 h(f0 + 4f1 + f2 ) f0 = f(x0 , y0 ) f1 = f(x0 + 21 h, y0 + 21 hf0 ) f2 = f(x0 + h, y0 + hf1 )

PD Orde 2

Bentuk umum PD orde 2:

d2 y y''= 2 = f(x, y, y') dx

Diketahui:

y(x0 ) = y0 , y'(x0 ) = y'0

Definisikan fungsi baru u:

u = y' u0 = y'0

Masalah PD orde 2 berubah menjadi masalah PD orde 1.

y(x) = ?

y'= u(x, y) u'= f(x, y, u)

Contoh penyelesaian dengan metode Euler yang lebih baik (improved):

u'= f(x, y, u)

y'= u(x, y)

u(x0 + h) = u0 + 21 h(f0 + f1 ) f0 = f(x0 , y0 , u0 )

y(x0 + h) = y0 + 21 h(u0 + u1 ) u0 = y'0

f1 = f(x0 + h, y0 + hu0 , u1 )

u1 = u0 + hf0

Alur perhitungan:

y0 , u0

f0

u1

f1 , y(x0 + h), u(x0 + h)

x0 + h → x0 , u(x0 + h) → u0 , y(x0 + h) → y0

Contoh penyelesaian dengan metode Runge-Kutta orde 4:

u'= f(x, y, u)

y'= u(x, y)

u(x0 + h) = u0 + 61 h(f0 + 2f1 + 2f2 + f3 ) f0 = f(x0 , y0 , u0 ) f1 = f(x0 + 21 h, y0 + 21 hu0 , u1 ) f2 = f(x0 + 21 h, y0 + 21 hu1 , u2 )

y(x0 + h) = y0 + 61 h(u0 + 2u1 + 2u2 + u3 ) u0 = y'0 u1 = u0 + 21 hf0 u2 = u0 + 21 hf1

f3 = f(x0 + h, y0 + hu2 , u3 )

u3 = u0 + hf2

Alur perhitungan:

y0 , u0

f0

u1

f1

u2

f2

u3

x0 + h → x0 , u(x0 + h) → u0 , y(x0 + h) → y0

f3 , y(x0 + h), u(x0 + h)

PD dengan Syarat Batas Contoh, gelombang yang merambat di sepanjang tali bisa digambarkan dengan PD orde 2. Jika ujung-ujung tali itu diikat sehingga tidak bisa bergerak, maka kita temui kasus PD dengan syarat batas.

terikat

terikat

Bentuk umum PD orde 1 & 2 linear:

Diketahui:

(1) y'= f(x, y) = d(x) − e(x)y (2) y''= g(x, y, y') = a(x) − b(x)y − c(x)y'

x0 ≤ x ≤ xn y(x0 ) = y0

y(x) = ?

y(xn ) = yn Dicari yi = y(xi ) pada titik xi = x0 + ih (i = 1, ..., n − 1) dengan h =

xn − x0 . n

Metode Finite Differences (1) y'+e(x)y = d(x)

yi+1 − yi−1 2h y − 2yi + yi−1 y''i ≅ i+1 h2 y'i ≅

(2) y''+c(x)y'+b(x)y = a(x) (1) y'i +ei yi = di (2) y''i +ci y'i +bi yi = ai

yi+1 − yi−1 + ei yi ≅ di 2h y − 2yi + yi−1 yi+1 − yi−1 (2) i+1 + c + bi yi ≅ ai i 2 h 2h

(1)

Jadi, pada akhirnya ditemui masalah sistem persamaan linear:

(1) − yi−1 + 2eihyi + yi+1 ≅ 2dih  ch   ch  (2)  1 − i  yi−1 − 2 − bih2 yi +  1 + i  yi+1 ≅ aih2 2  2   

(

)

yang dapat diselesaikan menggunakan metode, contoh, Eliminasi Gauss. Namun, sistem persamaan linear juga dapat diselesaikan dengan cara lain yaitu, iterasi. Contoh: iterasi Jacobi dan iterasi Gauss-Siedel.

Iterasi Jacobi sistem persamaan linear:

n

∑a x ij

j

= bi (i = 1, ..., n)

j=1

n  1  solusi: xi = bi − ∑ aijxj   aii  j≠i 

Pencarian solusi dimulai dengan nilai awal xi(0) (i = 1, …, n) hasil perkiraan / tebakan. Dengan nilai tebak awal ini diperoleh nilai perkiraan berikut xi(1) melalui: (1) i

x

n  1  (0) =  bi − ∑ aijxj  (i = 1, ..., n) aii  j≠i 

Demikian seterusnya berulang-ulang, nilai perkiraan pada langkah ke k diperoleh dari nilai perkiraan pada langkah ke k-1: (k) i

x

n  1  =  bi − ∑ aijxj(k-1)  (i = 1, ..., n) aii  j≠i 

Pencarian dihentikan setelah didapat nilai xi yang konvergen yaitu, yang tidak atau sedikit berubah dari nilai yang diperoleh pada langkah sebelumnya:

xi(k-1) 1 − (k) < ε, xi

ε = bilangan kecil

Iterasi Gauss-Siedel

Rumus iterasi Jacobi dapat ditulis:

(k) i

x

 1  (k -1) (k -1) =  bi − ∑ aijxj − ∑ aijxj  aii  ji 

Jika pada tiap langkah pencarian dilakukan dengan urutan i yang makin besar, (k) (k) maka semua xji (k) sudah diperoleh ketika mencari xi . (k) (k) Karena itu, nilai xji itu bisa langsung dipakai untuk mencari xi(k), sehingga iterasi mencapai nilai konvergen menjadi lebih cepat:

 1  (k) (k -1)  b − a x − a x ∑ ∑ i ij j ij j   aii  ji 

(i = 1, 2, ..., n)

 1  (k -1) (k) =  bi − ∑ aijxj − ∑ aijxj  aii  ji 

(i = n, ..., 2, 1)

xi(k) = (k) i

x

Iterasi seperti ini disebut iterasi Gauss-Siedel.

Sebelum ini dikenal metode Eliminasi Gauss dan LU Decomposition untuk mencari solusi sebuah sistem persamaan linear. Pada metode ini terdapat substitusi mundur dan maju. Pada substitusi mundur (maju), nilai xi dihitung dari nilai xj>i ( xji ( xj
Aplikasi Iterasi Jacobi dan Gauss-Siedel pada PD dengan Syarat Batas PD orde 1: (1) y'= d(x) − e(x)y

− yi−1 + 2eihyi + yi+1 ≅ 2dih

PD orde 2: (2) y''= a(x) − b(x)y − c(x)y'

 cih   ch  1 −  yi−1 − (2 − bih2 )yi +  1 + i  yi+1 ≅ aih2 2  2   

Iterasi Jacobi:

(1) yi(k) ≅

1 2dih + yi(k-1-1) − yi(k+1-1) 2eih

(2) yi(k) ≅

 1  cih  (k-1)  cih  (k-1)  2   − a h + (2 − bih2 )  i 1 − 2 yi-1 + 1 + 2 yi+1 

(

)

Iterasi Gauss-Siedel (contoh untuk i membesar, i = 1, …, n-1):

(1) yi(k) ≅

1 (k -1) 2dih + yi(k) -1 − yi+1 2eih

(2) yi(k) ≅

 1  cih  (k)  cih  (k-1)  2  − a h +  yi+1  1 −  yi-1 +  1 + i 2  (2 − bih )  2  2    

(

Catatan, sesuai syarat batas:

y0(k) = y0

)

yn(k) = yn

Persamaan Differensial Parsial

Pada bagian ini disampaikan metode numerik untuk menyelesaikan persamaan differensial parsial 2 dimensi tipe eliptik, parabolik dan hiperbolik.

Persamaan Differensial Eliptik

Bentuk umum PD eliptik:

  ∇2ψ(r ) = −4π ρ (r )

Untuk kasus 2 dimensi:

 ∂2 ∂2   2 + 2 ψ(x, y) = −4π ρ (x, y)  ∂x ∂y  

Gunakan metode finite differences:

ψ(xi+1 , yj ) − 2ψ (xi , yj ) + ψ(xi-1 , yj ) ψi+1, j − 2ψi, j + ψi-1, j ∂2 ≈ = ψ(x, y) 2 2 2 ∂x h h x ,y i

j

ψ(xi , yj+1 ) − 2ψ (xi , yj ) + ψ(xi , yj-1 ) ψi, j+1 − 2ψi, j + ψi, j-1 ∂2 ψ(x, y) ≈ = 2 2 ∂y2 h h x ,y i

j

(h = xi+1 − xi = yj+1 − yj )

Diperoleh:

[

ψi, j = h2 π ρi, j + 41 ψi+1, j + ψi-1, j + ψi, j+1 + ψi, j-1

]

Dicari distribusi spasial ψ .

i,j+1

Langkah: 1. Buat grid pada bidang xy, dengan jarak terdekat antar titik h.

i,j

i-1,j

i,j-1

i+1,j

2. (Dianggap nilai pada batas-batas bidang xy diketahui.) Hitung dengan rumus :

y

[

ψi, j = h2 π ρi, j + 41 ψi+1, j + ψi-1, j + ψi, j+1 + ψi, j-1

]

secara berurutan ψi, j untuk i = 1 & j = 1, 2, 3, ..., lalu i = 2 & j = 1, 2, 3, ..., i = 3 & j = 1, 2, 3, ... dan seterusnya. 3. Jika dalam langkah 2 ditemui nilai ψi, j yang belum diketahui, gunakan nilai x tebakan.

h h

ψ0,0

4. Ulangi langkah 2 – 3 sampai dicapai kestabilan untuk nilai ψi, j di semua titik:

ψi, j

(iterasi sebelum)

− ψi, j

(iterasi berikutnya)

< ε (bilangan kecil)

Persamaan Differensial Parabolik

Bentuk umum PD parabolik:

  2 1 ∂   ∇ − ψ(r , t) = −4π ρ (r , t) Γ ∂ t  

Untuk kasus 2 dimensi:

 ∂2 1 ∂   2 − ψ(x, t) = −4π ρ (x, t) ∂ x Γ ∂ t  

Gunakan metode finite differences:

ψ(xi+1 , tj ) − 2ψ (xi , tj ) + ψ(xi-1 , tj ) ψi+1, j − 2ψi, j + ψi-1, j ∂2 ψ(x, t) ≈ = 2 2 h ∂x2 h x x x ,t i

j

ψ(xi , tj+1 ) − ψ (xi , tj ) ψi, j+1 − ψi, j ∂ ψ(x, t) ≈ = ∂t h ht xi , tj t (hx = xi+1 − xi , ht = tj+1 − tj )

Diperoleh: ψi, j+1 = 4Γ ht π ρi, j + ψi, j +

Γht ψi+1, j − 2ψi, j + ψi-1, j 2 hx

(

)

Untuk tiap posisi dicari perubahan ψ terhadap waktu.

Langkah:

i,j+1

i-1,j

i,j

1. Buat grid pada bidang xt, dengan lebar hx untuk arah x dan ht untuk arah t. 2. (Dianggap nilai awal dan nilai pada batas-batas daerah x diketahui.)

i+1,j

Hitung dengan rumus : t

ψi, j+1 = 4Γ ht π ρi, j + ψi, j +

Γht ψi+1, j − 2ψi, j + ψi-1, j 2 hx

(

)

secara berurutan ψi, j+1 untuk j = 0 & i = 1, 2, 3, ..., lalu j = 1 & i = 1, 2, 3, ..., j = 2 & i = 1, 2, 3, ... dan seterusnya. Kasus khusus:

ht

x

hx ψ0,0

Jika ρ dan nilai pada batas-batas daerah x tetap (tidak bergantung waktu), maka akan tercapai suatu waktu t, bahwa ψi, j+1 tidak berubah lagi (atau berubah hanya sedikit, sehingga dapat diabaikan):

ψi, j+2 − ψi, j+1 < ε (bilangan kecil)

Persamaan Differensial Hiperbolik

Bentuk umum PD hiperbolik:

  2 1 ∂2    ∇ − 2 2 ψ(r ) = −4π ρ (r ) c ∂t  

Untuk kasus 2 dimensi:

 ∂2 1 ∂2   2 − 2 2 ψ(x, t) = −4π ρ (x, t)  ∂x c ∂t 

Gunakan metode finite differences:

ψ(xi+1 , tj ) − 2ψ (xi , tj ) + ψ(xi-1 , tj ) ψi+1, j − 2ψi, j + ψi-1, j ∂2 ψ(x, t) ≈ = 2 2 2 ∂x h h x x x ,t i

j

ψ(xi , tj+1 ) − 2ψ (xi , tj ) + ψ(xi , tj-1 ) ψi, j+1 − 2ψi, j + ψi, j-1 ∂2 ψ(x, t) ≈ = 2 2 ∂t2 h h t t xi ,tj Untuk tiap posisi dicari (hx = xi+1 − xi , ht = tj+1 − tj ) perubahan ψ terhadap 2 2 c ht waktu. 2 2 ψ = 4c h π ρ + 2ψ − ψ + ψ − 2ψ + ψ Diperoleh: i, j+1 t i, j i, j i, j-1 i + 1, j i, j i 1, j hx2

(

)

Untuk j = 0 diperoleh ψi,1 sebagai berikut: 2 2 t

ψi,1 = 4c h π ρi,0 + 2ψi,0 − ψi,-1

c2ht2 + 2 (ψi+1,0 − 2ψi,0 + ψi-1,0 ) hx

? Anggap

∂ ∂t

ψ(x, t) pada semua x dan t = 0 diketahui:

∂ ψ(x, t) = bi ∂t xi ,t0

∂ ψ(xi , t1 ) − ψ(xi , t−1 ) ψi,1 − ψi,-1 ψ(x, t) ≈ = Maka gunakan: , shg: ψi,-1 = ψi,1 − 2biht ∂t 2h 2h xi ,t0 t t Dengan demikian:

•

•

untuk j = 0:

c2ht2 (ψi+1,0 − 2ψi,0 + ψi-1,0 ) ψi,1 = 2c h π ρi,0 + biht + ψi,0 + 2hx2

untuk j > 0: ψi, j+1

2 2 t

c2ht2 = 4c h π ρi, j + 2ψi, j − ψi, j-1 + 2 ψi+1, j − 2ψi, j + ψi-1, j hx 2 2 t

(

)

i,1

Langkah: 1. Buat grid pada bidang xt, dengan lebar hx untuk arah x dan ht untuk arah t.

i-1,0 i,0

i+1,0

i,j+1

2. (Dianggap nilai awal dan nilai pada batas-batas daerah x diketahui.) Hitung dengan rumus :

i-1,j t

i,j

c2ht2 (ψi+1,0 − 2ψi,0 + ψi-1,0 ) ψi,1 = 2c h π ρi,0 + biht + ψi,0 + 2 2hx 2 2 t

i+1,j

i,j-1

ψi, j+1

c2ht2 = 4c h π ρi, j + 2ψi, j − ψi, j-1 + 2 ψi+1, j − 2ψi, j + ψi-1, j hx 2 2 t

(

secara berurutan ψi, j+1 untuk j = 0 & i = 1, 2, 3, ..., lalu j = 1 & i = 1, 2, 3, ..., j = 2 & i = 1, 2, 3, ... dan seterusnya.

ht

ψ0,0

x

hx

)

Persamaan Eigenvalue Contoh, lagi, gelombang pada tali yang kedua ujungnya diikat. Pada suatu waktu simpangan di sepanjang tali y(x) memenuhi PD orde 2:

d2 f(x) 2 y(x) = ky(x) dx dengan k berhubungan dengan frekwensi, yang nilainya tidak sembarang, yang menunjukkan modus gelombang. Untuk tiap-tiap modus/frekwensi/k yang mungkin, berlaku simpangan y(x) tertentu. Dengan kata lain, k merupakan eigenvalue untuk eigenfunction y(x). Persamaan di atas disebut persamaan eigenvalue. f Dengan metode Finite Differences, PD di atas menjadi: i2 (yi+1 − 2yi + yi−1 ) = kyi h yang membentuk persamaan matriks:

⋱ ⋱  ⋱ ⋱  fi  2 h    

⋱ − 2fi h2 ⋱

 ⋮   ⋮       yi−1  y  i−1   fi      yi  = k yi  2 h   yi+1  ⋱ ⋱ yi+1     ⋮    ⋱ ⋱ ⋮ 

k=? y=?

Metode Pangkat (Power Method) Persamaan eigenvalue:

A ui = λi ui,

ui = eigenfunct ion, λi = eigenvalue

Jika A matriks n x n, maka i = 1, …, n. Sebagai eigenfunction (atau disebut juga eigenvector), ui bersifat orthogonal dan juga komplit yaitu, sembarang fungsi (vector) x dapat dtulis sebagai kombinasi linear ui : n

komplit: x = ∑ ciui

T i j

orthogonal: u u = δij

i=1

Bermula dengan sembarang vector x, dilakukan iterasi berikut:

Ax= y

Untuk kali pertama: Setelah m kali iterasi diperoleh:

y→x

y y

(1)

(m)

n

n

n

i=1

i=1

i=1

= A x = A∑ ciui = ∑ ciAui = ∑ ciλiui m

m

=A x =A

n

n

n

∑cu = ∑cA u = ∑c λ i i

i=1

m

i

i=1

i

i=1

m i i i

u

λi≠k <1 λk

Anggap λk merupakan eigenvalue terbesar: Maka, jika m besar (banyak iterasi):

y

(m)

 λim  = A x = c λ u + ∑ c λ u = λ  ck uk + ∑ ci m ui  ≅ ck λkmuk λk  i ≠k i ≠k  m

m k k k

m i i i

m k

λk diperoleh dengan jalan: T

x y

(m)

T

m

m k k

=x A x≅c λ

n

∑cu

T i i k

m k k

u ≅c λ

i=1

n

∑cδ

i ik

x T y (m) λk = T (m-1) x y

2 m k k

≅c λ

i=1

uk diperoleh dengan jalan: y

(m)

2

=y

(m)T

y

(m)

(

) u u ≅ (c λ )

m 2 T k k k k

≅ cλ

m 2 k k

uk ≅

y (m) y (m)

Jika λk(m) merupakan nilai λk yang diperoleh setelah iterasi sebanyak m kali, maka iterasi dihentikan setelah dicapai nilai yang konvergen:

λk(m-1) 1 − (m) < ε, λk

ε = bilangan kecil

Untuk mencari eigenvalue terbesar kedua, hilangkan uk dari perhitungan. Jadi, dipakai vector awal baru x’: n

n

i=1

i=1

x'= (A − λk ) x = ∑ ci (A − λk )ui = ∑ ci (λi − λk )ui = ∑ ci (λi − λk )ui = ∑ diui Iterasi:

i ≠k

A x' = y'

Setelah m kali iterasi diperoleh:

(di ≡ ci (λi − λk ))

i≠ k

y'→ x'

y'(m) = Am x'= ∑ diλi ui m

i ≠k

Anggap λl merupakan eigenvalue terbesar kedua: sehingga setelah banyak iterasi: Memperoleh λl dan ul :

λi≠k ≠l <1 λl

y'(m) ≅ dlλlmul

x'T y'(m) λl = T (m-1) x' y'

ul ≅

y'(m) y'(m)

Pola yang sama berlaku untuk mencari eigenvalue terbesar berikutnya.

Metode Pangkat Kebalikan (Inverse Power Method) Dengan metode pangkat didapat eigenvalue terbesar. Untuk mencari eigenvalue terkecil digunakan metode pangkat kebalikan.

A-1A ui = λi A-1ui

A ui = λi ui

A-1ui = λi-1ui

Bermula dengan sembarang vector x, dilakukan iterasi berikut:

A-1 x = y

Setelah m kali iterasi diperoleh:

y→x

y

(m)

n

= (A ) x = ∑ ci (λi-1 ) m ui -1 m

i=1

Jika λs eigenvalue terkecil, maka setelah banyak kali iterasi: Jadi, λs diperoleh sebagai:

xT y (m-1) λs = T (m) x y

dan us :

us ≅

y (m) y (m)

y (m) ≅ cs (λs-1 ) m us