Docstoc

METODE NUMERIK

Document Sample
METODE NUMERIK Powered By Docstoc
					               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)


                    1
    Contoh:    x-     =0            x2 = 1                   x = 1 dan -1
                    x

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


Pada dua contoh di atas akar fungsi dapat
                                                         x = 2/3 dan -3
dicari secara analitik.
Secara umum, tidak selalu begitu keadaannya.
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
                       ω                                       cos (t) = t
           celah



             sensor
                           plat
                              v                                      ?
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                 akar fungsi
    dengan cara memplot fungsi).
 2. Tentukan batas awal yang
    mengurung akar fungsi.
                                             a      df e c            b
 3. Belah dua daerah berisi akar fungsi
    itu.
 4. Tentukan daerah yang berisi akar
    fungsi.                                                    Batas e, f atau
                                                                   nilai di
 5. Ulangi langkah 3 dan 4 sampai                              tengahnya bisa
    dianggap cukup.                                            dipilih sebagai
 6. Tentukan akar fungsi.                                       akar fungsi.
•   Menentukan daerah yang berisi akar fungsi:
                                                             a         c   b
       Jika z merupakan akar fungsi, maka
       f(x < z) dan f(x > z) saling berbeda           f(x)
       tanda.
       f(a)*f(c) negatif, berarti di antara a & c                              x
                                                                   z
       ada akar fungsi.
       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.


                                    perkiraan sebelum - perkiraan berikut
         kesalahan relatif semu =
                                               perkiraan berikut
                            Kesalahan

               kesalahan mutlak = | perkiraan – nilai sebenarnya |

                                     perkiraan – nilai sebenarnya
               kesalahan relatif =
                                           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 |

                                 perkiraan sebelum - perkiraan berikut
      kesalahan relatif semu =
                                          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



                                                 garis lurus sebagai
                                                   pengganti f(x)
      akar fungsi yang diperoleh
         f(x)
                                          b

                                                        x
                             c
                   a


                                                p(x)



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

                                        af(b) − bf(a)
                p(c) = 0           c=
                                         f(b) − f(a)
Langkah:
 1. Perkirakan akar fungsi (bisa                                akar fungsi
                                     f(x)
    dengan cara memplot fungsi).
                                                                         b
 2. Tentukan batas awal yang
    mengurung akar fungsi.                                                    x
                                                            c
 3. Tarik garis lurus penghubung
                                                a
    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.        f(x)
    Ulangi langkah 3.
                                                                         b
 5. Ulangi langkah 4 sampai
    dianggap cukup.                                                           x
                                                                c
 6. Titik potong garis nol dan                              a
    garis lurus yang terakhir               af(b) − bf(a)
    dinyatakan sebagai akar          c=
    fungsi.                                  f(b) − f(a)
                                   f(x)
                                                       af(b) − bf(a)
                                                  c=
                                                        f(b) − f(a)

                                              b
                                          c
                                                                       x
Metode false position juga
menggunakan dua batas seperti        a
metode bisection. Namun, berbeda
dari metode bisection, pada
metoda false position hanya satu
batas yang berubah.                f(x)
Pada contoh sebelum ini, batas a
berubah sementara batas b tetap.
Pada contoh berikut terjadi                                  b
sebaliknya.                                   c
                                                                       x
                                     a
Menghitung akar fungsi dengan metode false position,
menggunakan a dan b sebagai batas awal:
 •   jika batas a tetap, batas b berubah:

                  af(xi ) − xif(a)
         xi+1 =                           (i = 0, 1, 2, ...;   x0 = b)
                   f(xi ) − f(a)

 •   jika batas b tetap, batas a berubah:

                  bf(xi ) − xif(b)
         xi+1 =                           (i = 0, 1, 2, ...;   x0 = a )
                   f(xi ) − f(b)

 •   kesalahan relatif semu:


                                       xi − xi+1
                             ∆ rel =
                                          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.




 f(x)                                 akar fungsi yang diperoleh


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


                            c
                                                 x

              a
                                          p(x)




Diperoleh:   p(x) = f(a) + (x − a)f'(a)

             (f’(a) turunan pertama f(x) pada x = a)


                                           f(a)
             p(c) = 0               c=a−
                                           f'(a)
Langkah:                             f(x)
 1. Perkirakan akar fungsi.
                                                              c
 2. Buat garis singgung pada titik                                x
    sesuai akar fungsi yang                    a
    diperkirakan itu, lalu cari
    titik potongnya dengan garis                akar fungsi
    nol.                                        sebenarnya
 3. Titik potong itu merupakan
    perkiraan akar fungsi baru.
 4. Ulangi langkah 2 dan 3 sampai
    dianggap cukup.                  f(x)
 5. Titik potong garis nol dan                                a
    garis singgung kurva yang
                                                                  x
    terakhir dinyatakan sebagai                     c
    akar fungsi.                            f(a)
                                     c=a−
                                            f'(a)
                                     1
                          f(x)
Contoh perkiraan akar
fungsi awal yang baik
  perkiraan akar fungsi                          x
  makin mendekati akar
  fungsi sebenarnya.


                                         2

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



                                             2
Menghitung akar fungsi dengan metode Newton-Raphson:

                        f(xi )
          xi+1 = xi −                  (i = 0, 1, 2, ...;   x0 = a )
                        f'(xi )



kesalahan relatif semu:


                                       xi − xi+1
                             ∆ rel =
                                          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:

       bf(x0 ) − x0 f(b)          bf(x1 ) − x1f(b)          bf(x2 ) − x2f(b)
x1 =                       x2 =                      x3 =
        f(x0 ) − f(b)              f(x1 ) − f(b)             f(x2 ) − f(b)


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

       bf(x0 ) − x0 f(b)          bf(x1 ) − x1f(b)          x1f(x2 ) − x2f(x1 )
x1 =                       x2 =                      x3 =
        f(x0 ) − f(b)              f(x1 ) − f(b)              f(x2 ) − f(x1 )


                                                            xi-2f(xi-1 ) − xi-1f(xi-2 )
Jadi, mulai dari i = 3, akar fungsi dihitung dengan: xi =
                                                                f(xi-1 ) − f(xi-2 )
f(x)                                f(x)
                            b                                   b
                                x                                   x
       a                                    x1

            x1                                   x2

   I                                  II


           false position                   secant
                                                           x3
f(x)                                f(x)
                            b
                                x                                   x
                   x2                       x1        x2

                   x3

  III                                 III
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:
                                                                                  -1
       x f(x ) − xi-1f(xi-2 )                               f(xi-1 ) − f(xi-2 ) 
   xi = i-2 i-1                                xi = xi-1 − 
                                                            x −x                 f(xi-1 )
                                                                                 
           f(xi-1 ) − f(xi-2 )                                   i-1    i-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 )
                      ≅          = f'(xi-1 )                    xi ≅ xi-1 −
      xi-1 − xi-2       dxi-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
                                         x2
                               x3
                                    x4

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

       ε1 = x2 − x1 , ε2 = x3 − x2 , ε3 = x4 − x3
                   ε2 = 2 ε1 , ε3 = 2 ε2
                        1           1




Kecepatan konvergensi bersifat linear:        εi+1 = 2 εi
                                                     1
Pada metode False Position, Newton-Raphson dan Secant akar fungsi dicari
dengan rumus yang bentuknya serupa:

                                                     -1
                                      f(xi ) − f(a) 
   False Position:     xi+1           x − a  f(xi )
                              = xi −                            (atau a diganti b)
                                           i        

                                     f(xi )
   Newton-Raphson:     xi+1 = xi −
                                     f'(xi )
                                                          -1
                                    f(xi ) − f(xi-1 ) 
   Secant:             xi+1 = xi − 
                                    x −x               f(xi )
                                                       
                                         i    i-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.
              f(xi )                               f(xi )                   f(xi+1 )
xi+1 = xi −                     εi ≡ xi − xi+1 =                   εi+1 =             =?
              f'(xi )                              f'(xi )                  f'(xi+1 )

ekspansi deret Taylor:

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



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

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

Kecepatan konvergensi pada metode Newton-                                         f''(xi ) 2
Raphson (kira-kira demikian juga False Position                         εi+1 ≅             εi
                                                                                  2f'(xi )
dan Secant) bersifat kurang lebih kuadratik:

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:
                                                   n buah
                                                 persamaan
a11x1   + a12x2   + a13x3 ⋯ + a1n xn     = b1   dengan n buah
                                                  unknown
a21x1   + a22x2   + a23x3 ⋯ + a2n xn     = b2          xj
a31x1   + a32x2   + a33x3 ⋯ + a3n xn     = b3
  ⋮        ⋮        ⋮      ⋱    ⋮         ⋮       aij dan bi
                                                  diketahui
an1x1   + an2x2   + an3x3 ⋯ + ann xn     = bn
                                                i, j = 1, 2, …, n

                                xj = ?
Soal:     2x − 3y + 2z = −6      (1)
                                               3 persamaan dan
          − x + 2y − 3z = 2      (2)           3 unknown
            x+ y− z =0           (3)

Jawab:   2x −     3y + 2z = −6   (1)
                                               eliminasi x:
                0.5y − 2z = −1   (2)           pers. (2) + 0.5 pers. (1)
                2.5y − 2z = 3    (3)           pers. (3) – 0.5 pers. (1)


         2x −    3y + 2z = −6    (1)
                                               eliminasi y:
                0.5y − 2z = −1   (2)           pers. (3) – 5 pers. (2)
                      8z = 8     (3)

                       z =1
                                               substitusi mundur:
                          − 1 + 2z             pers. (3)   mencari z
                        y=         =2
                            0.5                pers. (2)   mencari y
                          − 6 + 3y − 2z        pers. (1)   mencari x
                       x=               = −1
                                 2
Dalam bentuk matriks:

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

          Jawab:        2 −3     2  x   − 6 
                                      
                         0 0.5 − 2  y  =  − 1 
                         0 2.5 − 2  z   3 
                                      
                        2 −3     2  x   − 6 
                                      
                         0 0.5 − 2  y  =  − 1 
                        0    0   8  z   8 
                                      
                          z =1
                             − 1 + 2z
                          y=           =2
                                0.5
                              − 6 + 3y − 2z
                          x=                = −1
                                    2
                                Eliminasi Gauss
Metode Eliminasi Gauss mencari solusi sebuah sistem persamaan linear
dengan cara seperti ditunjukkan pada contoh sebelum ini:

    a11     a12    a13 ⋯ a1n  x1   b1                (0)
                                                         aij = aij , bi(0) = bi
                                
    a21     a22    a23 ⋯ a2n  x2   b2                                     (k
                                                                              aik -1) (k-1)
   a        a32    a33 ⋯ a3n  x3  =  b3            a (k)
                                                          ij     =a (k -1)
                                                                   ij        − (k-1) akj
    31                                                                  akk
    ⋮        ⋮      ⋮ ⋱ ⋮  ⋮   ⋮ 
   a
                                                                               (k
                                                                              aik -1) (k −1)
                    an3 ⋯ ann  xn   bn 
                                                           (k)      (k −1)
    n1      an2                                    bi      =bi         − (k-1) bk
                                                                              akk
                                                         (k = 1, ..., n − 1;
                                                         i = k + 1, ..., n; j = k, ..., n)
 a(0)
           a(0)
                   a(0)
                          ⋯ a  x1   b 
                                 (0)             (0)
  11       12
            (1)
                    13
                    (1)
                                 1n
                                     
                                 (1)
                                                 1
                                                  (1)
                                                     
  0       a       a      ⋯ a  x2   b 
 
            22      23           2n
                                          
                                                 2
                                                             aij , bi(m) ≡ aij , bi
                                                                (m)


  0        0      a(2)
                          ⋯ a  x3  =  b 
                                 (2)             (2)
                    33           3n            3           pada langkah ke m
  ⋮         ⋮       ⋮    ⋱   ⋮  ⋮   ⋮ 
                                                
  0        0       0        (n
                          ⋯ ann-1)  xn   bn(n-1)                  halaman berikut
Substitusi mundur:


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



              bn(n-1)
          xn = (n-1)
              ann
                                        n
                     (n- j-1)
                     b
                     n-j        −     ∑a
                                    k =n- j+1
                                                  (n- j-1)
                                                        x
                                                  n - j,k  k

          xn − j =                  (n- j-1)
                                                               (j = 1, ..., n − 1)
                                a   n - j,n - j
                 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.



                                                      U
       A       X = B                                       X = B
                                      L
Langkah:

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


                                    U                                  U
          A     =                                                             X = B
                     L                              L


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                              Y = B
                                        L


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

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

 l11 0 0         ⋯ 0         1 u12 u13     ⋯ u1n   a11         a12     a13    ⋯        a1n 
                                                                                           
l21 l22 0        ⋯ 0         0 1 u23       ⋯ u2n   a21         a22     a23    ⋯        a2n 
l l      l       ⋯ 0        0 0     1      ⋯ u3n  =  a31       a32     a33    ⋯        a3n 
 31 32 33                                                                                  
⋮ ⋮ ⋮            ⋱ ⋮        ⋮ ⋮     ⋮      ⋱ ⋮   ⋮              ⋮       ⋮     ⋱         ⋮ 
                                                      
 ln1 ln2 ln3     ⋯ lnn      0 0 0
                                             ⋯ 1   an1
                                                                   an2     an3    ⋯        ann 
                                                                                                

Diperoleh:

      ai1 = li1                        → li1 = ai1                        (i = 1, ..., n)
                                                a1j
      a1j = l11u1j                     → u1j =                            (j = 2, ..., n)
                                                 l11
      ai2 = li1u12 + li2               → li2 = ai2 − li1u12               (i = 2, ..., n)
                                                 a2j − l21u1j
      a2j = l21u1j + l22u2j            → u2j =                            (j = 3, ..., n)
                                                     l22
      ai3 = li1u13 + li2u23 + li3      → li3 = ai3 − li1u13 − li2u23      (i = 3, ..., n)
                                                 a3j − l31u1j − l32u2j
      a3j = l31u1j + l32u2j + l33u3j   → u3j =                            (j = 4, ..., n)
                                                          l33
      …
Jadi, elemen matriks L dan U dicari menurut:

       li1 = ai1               (i = 1, ..., n)
               a1j
       u1j =                   (j = 2, ..., n)
               l11
                     j−1
       lij = aij − ∑ likukj    (i = j, ..., n; j = 2, ..., n)
                     k =1
                       i−1
               aij − ∑likukj
                      k =1
       uij =                   (j = i + 1, ..., n; i = 2, ..., n - 1)
                       lii

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       ⋯ 0  y1   b1                     b1
                                            y1 =
   l21 l22 0      ⋯ 0  y2   b2                     l11
   l l      l     ⋯ 0  y3  =  b3                          i−1
    31 32 33                                       bi − ∑ lij yj
   ⋮ ⋮ ⋮          ⋱ ⋮  ⋮   ⋮                              j=1
                                            yi =                               (i = 2, ..., n)
    ln1 ln2 ln3   ⋯ lnn  yn   bn                           lii



Substitusi mundur untuk menghitung x:


   1 u12 u13      ⋯ u1n  x1   y1 
                           
   0 1 u23        ⋯ u2n  x2   y2      xn = yn
  0 0     1       ⋯ u3n  x3  =  y3                        n
  
  ⋮ ⋮     ⋮
                            
                   ⋱ ⋮  ⋮   ⋮ 
                                            xn-i = yn −i −     ∑u
                                                             j=n −i+1
                                                                        n −i, j   xj (i = 1, ..., n − 1)
  0 0 0                     
                  ⋯ 1  xn   yn 
                         
                                     2 −3    2        − 6
                                                      
Kembali ke soal AX = B , dengan A =  − 1 2 − 3 , B =  2  .
                                     1   1 − 1        0
                                                      


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

                                        − 3
                                        
          Y = UX, LY = B           Y =  − 2
                                        1
                                        

                            −1
                            
          UX = Y       X =  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        X    =    B     atau AX = B



Keterangan: •    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   a12   a13     ⋯ a1n  x11   b11       a11     a12   a13   ⋯ a1n  x12   b12 
                                                                             
 a21   a22   a23     ⋯ a2n  x21   b21       a21     a22   a23   ⋯ a2n  x22   b22 
a      a32   a33     ⋯ a3n  x31  =  b31    a        a32   a33   ⋯ a3n  x32  =  b32 
 31                                         31                               
 ⋮      ⋮     ⋮      ⋱ ⋮  ⋮   ⋮             ⋮        ⋮     ⋮    ⋱ ⋮  ⋮   ⋮ 
a      an2   an3     ⋯ ann  xn1   bn1      a        an2   an3   ⋯ ann  xn2   bn2 
 n1                                         n1                               




                     a11   a12   a13 ⋯ a1n  x11       x12   b11 b12 
                                                                      
                     a21   a22   a23 ⋯ a2n  x21       x22   b21 b22 
                    a      a32   a33 ⋯ a3n  x31       x32  =  b31 b32 
                     31                                               
                     ⋮      ⋮     ⋮ ⋱ ⋮  ⋮             ⋮   ⋮       ⋮ 
                    a      an2   an3 ⋯ ann  xn1       xn2   bn1 bn2 
                     n1                                               
Metode Eliminasi Gauss:

•     rumus triangulasi:
      (0)         (0)
    aij = aij , bir = bir                   (i, j = 1, ..., n; r = 1, ..., m)

                                                                                                        (m)   (m)
                       a  (k -1)                                                                      aij , bir ≡ aij , bir
      (k)   (k
    aij = aij -1) −      ik
                          (k -1)
                                    (k
                                   akj -1) (k = 1, ..., n − 1;i = k + 1, ..., n;                    pada langkah ke m
                       a kk

      (k)
                         (k
                       aik -1) (k −1)
    bir         (k
            = bir −1) − (k-1) bkr               j = k, ..., n;r = 1, ..., m)
                       akk

•     rumus substitusi mundur:
                           (n
                          bnr -1)
                     xnr = (n-1)                                                 (r = 1, ..., m)
                          ann
                                                        n
                                    b(n- j-1)
                                     n- j,r     −     ∑a
                                                    k =n - j+1
                                                                   (n- j-1)
                                                                   n- j,kx  kr

                     xn − j,r =                      (n- j-1)
                                                                                 (j = 1, ..., n − 1; r = 1, ..., m)
                                                 a   n - j,n - j
Metode LU Decomposition:

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


                           b1r
                   y1r =                                 (r = 1, ..., m)
                           l11
                                   i−1
                           bir − ∑ lij yjr
                                   j=1
                   yir =                                 (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
                                  j=n −i+1
                                             n −i, j   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:     2 −4   1  3  x1   2 
                                                        baris 2 ditukar
          −1  2 3 − 2  x2   2 
          3 −4         x  =  2                          dengan baris 3
                  1  2 3
                         
          1 − 3 −1  5  x4   2 
                         
Jawab:
         2 − 4    1     3  x1   2           2 − 4    1     3  x1   2 
                                                                  
         0   0    3.5 − 0.5  x2   3         0   2 − 0.5 − 2.5  x2   − 1
         0   2 − 0.5 − 2.5 3 x  =  − 1      0                   x  =  3 
                                                 0    3.5 − 0.5 3
          0 − 1 − 1 .5                                                 
                        3.5  x4   1
                                              0 − 1 − 1 .5
                                                                 3.5  x4   1
                                                                         
   2 − 4   1     3  x1   2           2 − 4   1     3  x1   2 
                                                                 
   0   2 − 0.5 − 2.5  x2   − 1       0   2 − 0.5 − 2.5  x2   − 1 
                             =                                       =
   0   0   3.5 − 0.5  x3   3         0   0   3.5 − 0.5  x3   3 
                                                                 
   0             2  x4   2           0   0 − 1.75        x   0.5 
       0   0                                       2.25  4       

                   x4 = 1                       − 1 + 0.5x3 + 2.5x4
                                           x2 =                     =1
                                                         2
                          3 + 0.5x4
                   x3 =             =1          2 + 4x2 − x3 − 3x4
                             3.5           x1 =                    =1
                                                         2
   2 −4   1  3       2        2 0 0 0                 1 − 2 0.5 1.5 
                                                                    
   −1  2 3 − 2       2        − 1 l22 0 0            0    1 u23 u24 
A=              B =  2     L=
                                    3 l32 l33 0 
                                                         U=
  
     3 −4  1  2
                                                        0   0 1   u34 
   1 − 3 −1          2        1 l                                    
             5                     42 l43 l44 
                                                          0
                                                                0 0 1    
    2 −4   1  3          2 0 0 0                 2 0 0 0
                                                              
    3 −4   1  2          3 2 0 0                −1 0 0 0
 A=                    L=                       L=
     −1  2 3 − 2           − 1 0 l33 0              3 2 l33 0 
                                                              
    1 − 3 −1              1 −1 l    l44           1 −1 l   l44 
              5
                                 43                     43     

         2        1 −2   0.5   1.5                    2 0     0      0
                                                                       
         2       0   1 − 0.25 − 1.25                  3 2     0      0
      B=       U=                                   L=
          2          0  0   1        u34                  − 1 0 3.5      0
                                                                       
         2       0   0   0      1                     1 − 1 − 1.75 l 
                                                                    44 


               2 0     0      0         1 −2   0.5    1.5 
         ¨                                                 
       2         3 2    0      0        0   1 − 0.25 − 1.25 
   r is 3 L =                         U=
 ba ris        − 1 0 3.5      0          0  0   1    − 1/7 
   ba         
               1 − 1 − 1.75
                                                            
                              2
                                
                                         0
                                             0   0      1    
                                                              
       Data Fitting dengan
       Metode Least Square

f(x)                  Keterangan:

           p(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 .
                  x
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:
                                                                             2
                   N         N                      
                                                    N           m      
              S = ∑ ∆i = ∑ (f(xi ) − p(xi ) ) = ∑   f(xi ) − ∑ ajxij 
                      2                      2

                                                i=1 
                                                                       
                  i=1    i=1                                   j= 0    


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)
                                            dg(x)             d2g(x)
                                dx x = a              = 0 dan     2
                                                                        >0
                                             dx x = a          dx x = a


                                 x
                a


Spesial: fungsi kuadratik g(x) = ax2 + bx + c

     dg(x)             
           = 2ax + b
      dx                g(x) memiliki satu titik minimun jika a > 0 atau
      2
     d g(x)             sebaliknya satu titik maksimum jika a < 0.
            = 2a       
      dx 2             
S merupakan fungsi kuadratik dalam aj (j = 0, …, m):

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

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


  ∂ 2S(a0 , ..., am )     N

             2
                      = 2∑ xi2k > 0       (k = 0, ..., m)
       ∂ak               i=1




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

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

                  ∑  ∑ xi aj = ∑ f(xi )xik
                  j= 0  i=1    i=1
                                                              (k = 0, ..., m)

                         N                      N
Definisikan:      ckj ≡ ∑ x    i
                                j +k
                                           bk ≡ ∑ f(xi )xik
                         i=1                    i=1


                                                                       m
maka diperoleh sebuah sistem persamaan linear:                        ∑c
                                                                      j= 0
                                                                              a = bk
                                                                             kj j      (k = 0, ..., m)


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.                    f(x)           p(x)
Jadi, dicoba fungsi kuadratik:
                                                           120
             p(x) = a0 + a1x + a2x2

Sistem persamaan linier untuk mencari aj :                  70

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


Jadi,     p(x) = 5x(x + 5 )            Cek:        p(1) = 30, p(2) = 70,
                                                                                 OK!
                                                         p(3) = 120
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]
                                                           y

                  10         0.01794775
                  15         0.03808997
                  20         0.05516225                        θ   E
                  25         0.05598281
                  30         0.04795629
                  35         0.04807485
                  40         0.06273566
                  45         0.07853982
                  50         0.07395442
                  55         0.04201338


        Cari fungsi p(x) yang dapat melukiskan data itu.
Dicoba beberapa polinomial dengan orde berbeda, diperoleh:

                                                  a0 = - 3.557800654975570E - 02
         a0 = 8.983713484853211E - 03             a1 = 1.061996221844471E - 03
         a1 = 1.324478388111303E - 03             a2 = 8.802185976358352E - 04
m = 3: a2 = 3.487808787880805E - 05       m = 5: a3 = - 5.862332690401015E - 05
         a3 = - 8.085809790211842E - 07           a4 = 1.362046192596346E - 06
         S = 1.0339E - 03                         a5 = - 1.063951754163944E - 08

                                                  S = 8.1573E - 05

         a0 = - 1.757260839248139E - 02
         a1 = 1.596300085173997E - 02             a0 = 1.864754537649403E - 01
         a2 = - 3.402768734407800E - 03           a1 = - 4.631839872868015E - 02
         a3 = 3.358961098305538E - 04             a2 = 4.007658091692495E - 03
         a4 = - 1.368895999268855E - 05           a3 = - 8.985715636865594E - 05
m = 9:   a5 = 1.132254508386570E - 07     m = 7: a4 = - 3.230489224228010E - 06
         a6 = 8.262829873458547E - 09             a5 = 1.912806006890119E - 07
         a7 = - 2.741786330789355E - 10           a6 = - 3.252863805243949E - 09
         a8 = 3.317446724324134E - 12             a7 = 1.876184315740421E - 11
         a9 = - 1.459511835946927E - 14
                                                  S = 3.1629E - 07
         S = 1.7528E - 11
       Interpolasi

f(x)              Keterangan:

       p(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 .
              x
                        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:
                            p(xi ) = f(xi ) (i = 1, ..., N )
Jadi, diperoleh persamaan linear:

                 p(x1 ) = a0 + a1 x1 + a2 x12 + ... + aN-1 x1N-1 = f(x1 )
                                          2              N
                 p(x2 ) = a0 + a1x2 + a2 x2 + ... + aN-1x2 -1 = f(x2 )
                                         2              N
                 p(x3 ) = a0 + a1x3 + a2x3 + ... + aN-1x3 -1 = f(x3 )
                 ...
                                         2               N
                 p(xN ) = a0 + a1xN + a2xN + ... + aN-1 xN -1 = f(xN )

dan ai (i = 0, …, N-1) diperoleh sebagai solusi dari persamaan linear itu.
N = 2:                             p(x1 ) = a0 + a1x1 = f(x1 )
                                   p(x2 ) = a0 + a1 x2 = f(x2 )
                               x2f(x1 ) − x1f(x2 )              f(x1 ) − f(x2 )
                      a0 = −                             a1 =
                                    x1 − x2                        x1 − x2
                                  x − x2           x − x1 
                          p(x) = 
                                 x −x              x − x f(x2 )
                                          f(x1 ) +          
                                  1    2           2     1 



N = 3:                          p(x1 ) = a0 + a1 x1 + a2 x12 = f(x1 )
                                                         2
                                p(x2 ) = a0 + a1x2 + a2 x2 = f(x2 )
                                                        2
                                p(x3 ) = a0 + a1x3 + a2x3 = f(x3 )

                   (x2 − x3 )x2 x3f(x1 ) + (x3 − x1 )x3x1f(x2 ) + (x1 − x2 )x1x2f(x3 )
            a0 =                                          2              2
                               (x2 − x3 )x12 + (x3 − x1 )x2 + (x1 − x2 )x3
                     2    2             2                         2
                   (x2 − x3 )f(x1 ) + (x3 − x12 )f(x2 ) + (x12 − x2 )f(x3 )
            a1 = −                                  2
                        (x2 − x3 )x12 + (x3 − x1 )x2 + (x1 − x2 )x3 2


                   (x2 − x3 )f(x1 ) + (x3 − x1 )f(x2 ) + (x1 − x2 )f(x3 )
            a2 =                                   2
                       (x2 − x3 )x12 + (x3 − x1 )x2 + (x1 − x2 )x3 2



        x − x2  x − x3           x − x1  x − x3            x − x1  x − x2 
p(x) =         
        x − x  x − x  f(x1 ) +         
                                     x − x  x − x              x − x  x − x f(x3 )
                                                        f(x2 ) +                  
        1    2  1    3           2 1  2        3           3 1  3        2 
                                                  N
Secara umum, untuk N data                p(x) = ∑ l(x, xi )f(xi )
                                                 i=1
rumus interpolasi Lagrange:

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


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 
                                      i        j 


                   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                    x − x1  x − x3  x − x4 
             x − x  x − x  x − x 
l(x, x1 ) =                                            x − x  x − x  x − x 
                                                l(x, x2 ) =                          
             1    2  1    3  1    4                    2     1  2    3  2    4 


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


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:
                                  2      3
                   a0 + a1xj + a2xj + a3xj = 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                                      x − xk    
                  p(x) = ∑ l(x, xj )f(xj )      l(x, xj ) = ∏                
                                                                  
                                                            k ≠ j  xj − xk
                                                                              
                           j= 0                                               
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           xi+5

                                          I          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 )
                           2      3
   p(x2 ) = a0 + a1x2 + a2x2 + a3x2 = f(x2 )
                               2
                                                  (xi ≤ x ≤ xi+1 ; x1 = xi , x2 = xi+1 )
   p'(x1 ) = a1 + 2a2x1 + 3a x = f'(x1 )
                             3 1
                              2
   p'(x2 ) = a1 + 2a2x2 + 3a3x2 = 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:


                    x12 (3x2 − x1 )f(x2 ) + x2 (x2 − 3x1 )f(x1 ) 
                                                 2
          a0   =   
                                                                 
                                                                  
                                     (x2 − x1 )3                 
                            x f'(x2 ) + x2 f'(x1 ) 
                   − x2 x1  1
                                             2
                                                       
                                                       
                                  (x2 − x1 )          
                              f(x2 ) − f(x1 ) 
          a1   =   − 6x2 x1  (x − x )3       
                                   2    1      
                      x1 (2x2 + x1 )f'(x2 ) + x2 (x2 + 2x1 )f'(x1 ) 
                   +                                2
                                                                     
                                                                     
                                         (x2 − x1 )                 
                                 f(x2 ) − f(x1 ) 
          a2   =   3(x2 + x1 )  (x − x )3        
                                     2     1       
                      (x + 2x1 )f'(x2 ) + (2x2 + x1 )f'(x1 ) 
                   − 2
                                                  2
                                                                
                                                                
                                     (x2 − x1 )                
                        f(x2 ) − f(x1 )   f'(x2 ) + f'(x1 ) 
          a3   =   − 2 (x − x )3  +  (x − x )2 
                                                             
                             2     1               2  1     
         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

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




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(xi-1 ) = f(xi-1 ) 
                                             
g(x) = ax 2 + bx + c      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:

                                 axi21 + bxi-1 + c = f(xi-1 )
                                    -

                                axi2 + bxi + c = f(xi )
                                axi2 1 + bxi+1 + c = f(xi+1 )
                                   +


diperoleh:

                     f(xi−1 )                   f(xi )                   f(xi+1 )
        a=                             +                        +
             (xi−1 − xi )(xi−1 − xi+1 ) (xi − xi−1 )(xi − xi+1 ) (xi+1 − xi−1 )(xi+1 − xi )
                 (xi + xi+1 )f(xi−1 )      (xi-1 + xi+1 )f(xi )      (xi-1 + xi )f(xi+1 )
        b=−                              −                        −
               (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            1               x − x  f(xi+1 )
            
  f'(xi ) ≅                           +          +          f(xi ) +  i i-1 
             xi−1 − xi+1  (xi−1 − xi )  xi − xi+1 xi − xi−1 
                                                             
                                                                          x − x  (x − x )
                                                                          i +1 i−1  i +1 i
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 )         
                          x − x  (x − x )
                          0     2  0   1

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

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

                          x − x1    1
  h(x, x3 ) = h2 (x, x2 ) 2      
                          x − x  (x − x )
                          3 1 3        2
                         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'' ≡ p''(xi ) = 2bi
                                                           i


dan pada titik x = xi+1:

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

sehingga diperoleh:


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

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

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 :

                  pi − pi-1 hi-1p'' hi-1p''-1                        p'' −p''-1 
        p'(x) =            −      i
                                    −      i
                                              + p''-1 (x − xi-1 ) +  i
                                                   i                 2h
                                                                               i
                                                                                 (x − xi-1 )2
                                                                                 
                     hi-1      6        3                                 i-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''-1 +2(hi-1 + hi )p'' +hip''+1 = 6
              i                  i       i      h − h               
                                                                        (i = 2, ..., N - 1)
                                                     i        i-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:

                                                    p −p          
                    (i = 1)      2h1p'' +h1p'' = 6 2 1 − p' 
                                       1       2     h           1
                                                        1         
                                                                  p −p 
                    (i = N)       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
                        p(x, y) = ∑ S(x, xi )∑ S(y, yj )f(xi , yj )
                                  i=1         j=1




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
     p(x, y) = ∑ l(x, xi )∑ l(y, yj )f(xi , yj )
                 i= 0         j= 0

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

                        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]
                                                         y

                  10         0.01794775
                  15         0.03808997
                  20         0.05516225                      θ       E
                  25         0.05598281
                  30         0.04795629
                  35         0.04807485
                  40         0.06273566
                  45         0.07853982
                  50         0.07395442
                  55         0.04201338


      Dengan interpolasi, cari nilai p(x) di sepanjang titik data.
                              Integrasi
Menghitung luas daerah di bawah kurva:


f(x)                                      f(x)
                analitik                                     numerik




                b

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


                                     x                                         x
       a                       b                 a                         b

           b        N              Integral numerik sering disebut juga sebagai
 I = ∫ f(x) dx ≅ ∑ wif(xi )        quadrature; integrasi numerik disebut sebagai
           a        i=1            integrasi dgn menjumlah quadrature.
Meski tidak terlihat pada rumus akhir, pada integrasi numerik integrand
f(x) diinterpolasi dengan suatu polinomial:

                         b           N
                     I = ∫ f(x) dx ≅ ∑ wif(xi )
                         a           i=1




                             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 ),   p(x) = r + sx
                   a            a          i=1


                                           f(x)

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


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

                                                          ∫ p(x) dx
                                                          a


                                                                              x
                                                  a                       b
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)
  a             i=1
                      i   i    1   1      2       2   1         2                 w1 , w2 = ?


Mencari w1 dan w2 :
                                   b

                                   ∫ (r + sx) dx = w (r + sa) + w (r + sb)
                                                      1             2
p(x) = r + sx                       a
                                    1
                          r(b - a) + s(b2 − a2 ) = r(w1 + w2 ) + s(aw1 + bw2 )
                                    2

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

                                              b
                                                          h
Rumus quadrature trapezoid:             I = ∫ f(x) dx ≅     (f(a) + f(b) )    (h = b − a)
                                              a
                                                          2


                                                          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 ),             p(x) = r + sx + tx2
                 a           a           i=1

                                               f(x)

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

     p(a) = f(a), p(b) = f(b),
                                                                     b
             p(c) = f(c)                                             ∫ p(x) dx
                                                                     a

                     a+c
    dengan     b=                                                                          x
                      2                               a                   b            c
Integrasi numerik dikerjakan untuk N = 3:
       c              3

       ∫ p(x) dx = ∑ w p(x ) = w p(a) + w p(b) + w p(c)
       a             i=1
                               i   i   1       2         3                  w1 , w2 , w3 = ?

Mencari w1 , w2 , w3:

p(x) = r + sx + tx2
                           c

                           ∫
                           a
                             (r + sx + tx2 ) dx = w1 (r + sa + ta2 ) + w2 (r + sb + tb2 )

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


        w1 + w2 + w3 = c - a                                                    1
                                                                   w1 = w3 =      (c − a)
                         1                                                      6
       aw1 + bw2 + cw3 = (c2 − a2 )
                         2                                                 2
                         1                                         w2 =      (c − a)
     a2w1 + b2w2 + c2w3 = (c3 − a3 )                                       3
                         3
                                                        c
                                                                      h
Diperoleh Rumus quadrature Simpson:                 I = ∫ f(x) dx ≅     (f(a) + 4f(b) + f(c) )
                                                        a
                                                                      3

             c−a
dengan h =       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
                    3h                                          d-a                        
  I = ∫ f(x) dx ≅      (f(a) + 3f(b) + 3f(c) + f(d) )       h =     = b − a = c − b = d − c
        a
                     8                                           3                         

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

    e
                2h                                                          e-a            
I = ∫ f(x) dx ≅    (7f(a) + 32f(b) + 12f(c) + 32f(d) + 7f(e) )           h=         = b− a
                45                                                           4             
    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
                                                                                                b−a
 I = ∫ f(x) dx =    ∫ f(x) dx +      ∫ f(x) dx + ... +    ∫ f(x) dx + ∫ f(x) dx              d=    
     a              a               a +d                 b-2d            b -d                    N 



 •   integrasi komposit menggunakan quadrature trapezoid
                              b
                          I = ∫ f(x) dx ≅ h[2 (f + f ) + f + f + ... + f −1 ]
                                            1
                                                0   N     1   2         N
                              a
                                   b−a
                              h=       , f = f(a + ih), i = 0, ..., N
                                          i
                                    N

 •   integrasi komposit menggunakan quadrature Simpson
            b
                             2h 1
         I = ∫ f(x) dx ≅        [2 (f0 + f2N ) + 2(f1 + f3 + ... + f2N−1 ) + f2 + f4 + ... + f2N−2 ]
             a
                              3
                                  b−a
                             h=       , f = f(a + ih), i = 0, ..., 2N
                                         i
                                  2N
  Integrasi komposit trapezoid untuk daerah integrasi [a,b] yang dibagi 8
  sama panjang:

                 b
             I = ∫ f(x) dx ≅ h[2 (f0 + f ) + f + f + f + f4 + f + f + f ]
                               1
                                        8     1   2   3        5   6   7
                 a


f(x)




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

                b
                              h1
            I = ∫ f(x) dx ≅      (fa + 2fa +h1 + fc ) + h2 (fc + 4fc+h2 + fb )
                a
                              2                          3
f(x)
                                                                                 Simpson
                     trapezoid




                                                                                       x
       a                                                 c                       b
                h1                        h1                       2h2
                            Quadrature Gaussian

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

                       b                                              n

                       ∫ v(x)O (x)O
                               n     m   (x) dx = δnm ,   On (x) = ∑ bixi
                       a                                             i= 0



Contoh:
                                                                1

   On =   2n +1
           2 nP , P = polinomial Legendre
                   n                                            ∫ O (x)O
                                                                -1
                                                                      n         m   (x) dx = δnm

                                                                ∞
   On = n! Ln , Ln =
        1
                           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 )
                  a                i=1
                                          i   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                                                                    N -1              2N -1
  f(x) ≅ p(x) =   ∑ aix = r(x) + s(x)
                  i= 0
                              i
                                                        dengan                r(x) = ∑ aix , s(x) =
                                                                                           i= 0
                                                                                                         i
                                                                                                             ∑ aixi
                                                                                                             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               i= 0

                          b             N -1            b                                         N -1
Maka:                     ∫ v(x)s(x)dx = ∑ c ∫ v(x)O (x)O
                          a             i= 0
                                                i
                                                        a
                                                                  i           N   (x)dx = ∑ ciδiN = 0
                                                                                                  i= 0

                          b              N                            N
Secara numerik:           ∫ v(x)s(x)dx = ∑ w s(x ) = ∑ w q(x )O
                          a             i=1
                                                    i       i
                                                                  i=1
                                                                          i        i   N   (xi ) = 0


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 )
                 a                   i=1
                                           i         i     ∫ v(x)s(x)dx = ∑ w s(x )
                                                           a                      i=1
                                                                                        i    i



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

                        ∫ v(x)r(x)dx = ∑ wir(xi ),
                        a                      i=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 
        r(x) = l(x, xi ) = ∏            , (i, j = 1, ..., N)                      l(xk , xi ) = δik
                                        
                           j≠i  xi − xj 


             b                      N                                 b
Diperoleh: ∫ v(x)l(x, xj )dx = ∑ wil(xi , xj ) = wj              wj = ∫ v(x)l(x, xj )dx (j = 1, ..., N)
             a                     i=1                                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 P :
                                                           n
                                        1

                      On =   2n +1
                              2 nP:     ∫ O (x)O
                                        -1
                                             n             m   (x) dx = δnm


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

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



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

           1          N                b
                                     2
Contoh:    ∫1 f(x)dx = ∑ wif(xi ) =     ∫ f(y)dy
                                                                      b            N

           -           i=1          b−a a                             ∫ f(y)dy = ∑ u f(y )
                                                                      a           i=1
                                                                                        i   i


                  y − a x − ( −1) x + 1                                       yi = (xi + 1)(b − a) + a
                                                                                  1
                                                                                  2
                       =          =
                  b − a 1 − ( −1)   2                                              b−a 
                                                                              ui =      wi
                  (transformasi linier)                                             2 
Contoh xi dan wi quadrature Gauss-Legendre untuk beberapa N terkecil:
                           1          N

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




  N                    x                                  w


  2         ± 0.577350269189626                    1.000000000000000


  3         ± 0.774596669241483                    0.555555555555556
              0.000000000000000                    0.888888888888889


  4         ± 0.861136311594053                    0.347854845137454
            ± 0.339981043584856                    0.652145154862546


  5         ± 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        Untuk f(x) yang berkurva tajam di
integrand f(x) yang bentuk kurvanya     bagian tengah dan kurang tajam di
lebih tajam di sekitar batas            sekitar batas integral diperlukan
integral, sementara kurang tajam di     beberapa penanganan (mis. membagi
bagian tengah.                          daerah integrasi, redistribusi x dll).
                 Quadrature Gauss-Laguerre

Quadrature Gauss-Laguerre menggunakan polinomial Laguerre Ln :
                               ∞
                     1
                On = n! Ln :   ∫
                               0
                                 e -xOn (x)Om (x) dx = δnm
                                      ∞                   N
dipakai untuk integral berbentuk:     ∫e
                                           -x
                                                f(x)dx = ∑ wif(xi )
                                      0                  i=1


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

N                       x                                       w

2           0.585786437626905                        0.853553390593274
             3.414213562373095                       0.146446609406726

4           0.322547689619392                        0.603154104341634
             1.745761101158347                       0.357418692437800
             4.536620296921128                       0.038887908515005
             9.395070912301133                       0.000539294705561
                               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.

               ∞
                 dx            batas integral sampai tak behingga, jika
Contoh:    I=∫                 dievaluasi langsung memerlukan sangat banyak
               1 + x2
             0                 titik, tidak praktis dan hasilnya bisa saja buruk

                                   1+ y           2
           transformasi:      x=        , dx =         2
                                                         dy
                                   1− y        (1 − y)
                                   1
                                               2dy
                              I=∫
                                   -1
                                               2
                                        (1 − y) 1 +(   ( ))
                                                       1+ y 2
                                                       1− y
                                        1
                                                   dy            quadrature
                               = 2∫
                                           (1 − y)2 + (1 + y)2   Gauss-Legendre
                                        -1
                             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)                 fungsi ganjil: f( −x) = −f(x)


                     a     2N
                      dx         wi
Contoh:      • I=∫        =∑
                 -a
                    1 + x2 i=1 1 + xi2
                         a        2N
                         dx              wi
                  = 2∫        =2 ∑
                     0
                       1 + x2               2
                                i=N +1 1 + xi


                    1                   2N
                            dy                        wi
             • I=∫                     =∑
                 -1
                    (1 − y)2 + (1 + y)2 i=1 (1 − yi )2 + (1 + yi )2
                         1                     2N
                               dy                              wi
                  = 2∫                     =2 ∑
                     0
                       (1 − y)2 + (1 + y)2                   2
                                             i=N +1 (1 − yi ) + (1 + yi )
                                                                          2
Beberapa fungsi memiliki simetri, contoh fungsi trigonometri:

                              sin( −x) = −sin(x)               sin(π ± x) = ∓ sin(x)
                              cos( −x) = cos(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π
       [f(sin(x − a)) + f(cos(x − a))]eim(x −a)dx                             integral tertutup bisa
I=   ∫
     0
                                                                              dimulai dari titik mana saja
     2π

     ∫ [f(sin(x)) + f(cos(x))]e
                                        imx
 =                                            dx
     0
     π                                                                                               telah
         [
 = ∫ {f(sin(x)) + f(cos(x))}e               imx
                                                  + {f( −sin(x)) + f( −cos(x))}e   im(x + π)
                                                                                               ]dx   dipakai
     0
     π
         2
                                                                                                     x = x'+ π
     ∫ [f(sin(x))(e                     )                 (               )
                        imx
 =                            + eim(π −x) + f(cos(x)) eimx + eim(2π −x)                              x = −x'
     0

                        (                          )              (
         + 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
                       dx
Contoh:   • I=∫                                            singular pada x = 0
                 0 (1 + x) x
                 a             a      a
                      dx       dx    dx
               =∫           −∫    +∫                       ditambah nol
                0 (1 + x) x  0  x 0 x
                 a
                         1      1 
                                         a
                                           dx              subtraksi pada
               = ∫           −   dx + ∫
                 0  (1 + x) x    x
                                        0  x              integral asal
                     a
                         x
               = −∫        dx + 2 a
                     0
                       1+x
     a
      x2f(x)dx
• I=∫ 2                          (0 < b2 ≤ a2 )                   singular pada x = b
    0
      (b − x2 )
     a                       ∞
       x2f(x)dx     b2f(b)dx
    =∫ 2
       (b − x2 ) ∫ (b2 − x2 )
                −                                                 ditambah nol (lihat *)
     0            0
     a
    =∫
         (x f(x) − b f(b))dx −
           2             2               ∞
                                           b2 f(b)dx              subtraksi pada
     0
                   (b2 − x2 )            ∫ (b2 − x2 )
                                         a
                                                                  integral asal
     a
    =∫
         (x f(x) − b f(b))dx − 1 bf(b)ln a − b 
           2             2
                                               
     0
                   (b2 − x2 )            2          a +b 



               ∞                     ∞                        ∞         ∞
                    dx       1  1         1        1  dx    1 dx
     (*)       ∫ (b2 − x2 ) 2b ∫  b − x b + x  2b −∫∞b + x 2b −∫∞ x = 0
               0
                           =
                               0
                                        +
                                               
                                                dx =        =
                      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
a
                                 2n   sin(tb) − f sin(ta)) + β(th)Cgenap + γ(th)Cganjil
                                                 0                                        ]
                       n
                                                                                      b−a
             Cgenap = ∑ f cos(tx2i ) − 2 (f cos(tb) + f cos(ta) )
                         2i
                                       1
                                           2n          0                          h=      = xi+1 − xi
                      i= 0                                                             2n
                       n
                                                                                  x0 = a
             Cganjil = ∑ f −1cos(tx2i−1 )
                          2i
                      i=1                                                         fj = f(xj )
b

∫ f(x)sin(tx)dx = h[α(th)(f cos(ta) − f
a
                                0               2n   cos(tb)) + β(th)Sgenap + γ(th)Sganjil    ]
                       n
             Sgenap = ∑ f sin(tx2i ) − 2 (f sin(tb) + f sin(ta) )
                         2i
                                       1
                                           2n          0
                      i= 0
                       n
             Sganjil = ∑ f −1 sin(tx2i−1 )
                          2i
                      i=1


       1 sin(2x) 2sin2x                                                 2x3 2x 5 2x 7
α(x) = +            −                                            α(x) =    −    +     − ...
       x      2x2      x3                                               45 315 4725
         1 + cos2x sin(2x)                   untuk nilai              2 2x2 4x 4 2x 6
β(x) = 2
               2
                    −    3
                            
                                                                β(x) = +    −    +       − ...
             x        x                      x kecil:                 3 15 105 567
         sinx cosx                                                    4 2x2 x 4     x6
γ(x) = 4 3 − 2                                                 γ(x) = −     +   −          + ...
         x        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

                                <f(x)>
                                                b
                                                                  1 n
                                            I = ∫ f(x)dx ≅ (b - a) ∑ f(xi )
            (b-a)<f(x)>                         a
                                                                  n i=1
                                   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
Bentuk umum PD orde 1:        y'=      = f(x, y)
                                    dx

Diketahui:                    y(x0 ) = y0                y(x) = ?
                              y       x
                                                                      Masalah persamaan
Integrasi:                    ∫ dy = ∫ f(x, y)dx
                             y0       x0
                                                                      differensial
                                                                      berubah menjadi
                                                                      masalah persamaan
                                              x                       integral.
                              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)
                                                  y(x)
                                                                       yg diperoleh
                               x0 +h

y(x0 + h) ≅ y0 + f(x0 , y0 )    ∫ dx
                                x0
                                                                           y(x0 + h)
          ≅ y0 + hf(x0 , y0 )                    y0
                                                                          sebenarnya
                                                                      x
                                                          x0 x0 + h
                          Metode Euler yang Dimodifikasi

                                                     f(x,y)

Modifikasi dilakukan dalam memilih                                       f(x0 + 2 h, y(x0 + 2 h))
                                                                                1           1

nilai f(x,y) yang dianggap konstan.
Dipilih f(x,y) pada titik x = x0 + 2 h :
                                   1

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

dengan y(x0 + 2 h) dihitung memakai
              1                                               x0 + 2 h
                                                                   1


metode Euler:
                                                                              y(x0 + h)
     y(x0 + 2 h) ≅ y0 + 2 hf(x0 , y0 )
            1           1                            y(x)
                                                                             yg diperoleh

Diperoleh:

y(x0 + h) ≅ y0 + hf(x0 + 2 h, y(x0 + 2 h))
                         1           1
                                                                                 y(x0 + h)
                                                     y0
          ≅ y0 + hf(x0 + 2 h, y0 + 2 hf(x0 , y0 ))
                         1         1
                                                                                sebenarnya
                                                                            x
                                                              x0 x0 + h
                                                              x0 + 2 h
                                                                   1
                     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 x0 + h
dengan y(x0 + h) dihitung memakai metode Euler:
                                    y(x0 + h) ≅ y0 + hf(x0 , y0 )

Diperoleh:

                     y(x0 + h) ≅ y0 + 2 h[f(x0 , y0 ) + f(x0 + h, y(x0 + h))]
                                      1


                               ≅ y0 + 2 h[f(x0 , y0 ) + f(x0 + h, y0 + hf(x0 , y0 ))]
                                      1
                              Metode Runge-Kutta

Metode Euler dan variasinya sebelum ini sebetulnya termasuk metode Runge-
Kutta, 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:

                                                             2
Euler                   : y(x0 + h) = y0 + hf(x0 , y0 ) + O(h )
Euler yg dimodifikasi: y(x0 + h) = y0 + hf(x0 + 2 h, y0 + 2 hf(x0 , y0 )) + O(h3 )
                                                1         1


Euler yg lebih baik     : y(x0 + h) = y0 + 2 h[f(x0 , y0 ) + f(x0 + h, y0 + hf(x0 , y0 ))] + O(h3 )
                                           1



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

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

      x0

                      ≅ 6 h(f + 2f + 2f + f )
                        1
                             0    1    2   3



dengan:            f = f(x0 , y0 )
                    0                                f = f(x0 + 2 h, y(x0 + 2 h))
                                                      2
                                                                1           1


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


f dan f memiliki nilai berbeda, karena dihitung untuk nilai argumen y(x0 + 2 h)
 1     2
                                                                           1

yang berbeda: menurut metode Euler, y(x0 + 2 h) dapat diperoleh melalui 2
                                              1

persamaan:

             (1)    y(x0 + 2 h) ≅ y0 + 2 hf0
                           1           1
                                                           f = f(x0 + 2 h, y0 + 2 hf0 )
                                                            1
                                                                      1         1


             atau
             (2) y0 ≅ y(x0 + 2 h) − 2 hf(x0 + 2 h, y(x0 + 2 h))
                             1      1         1           1


                        ≅ y(x0 + 2 h) − 2 hf(x0 + 2 h, y0 + 2 hf )
                                 1      1         1         1
                                                                0

                          y(x0 + 2 h) ≅ y0 + 2 hf
                                 1           1
                                                 1          f = f(x0 + 2 h, y0 + 2 hf )
                                                             2
                                                                       1         1
                                                                                     1
Untuk f , digunakan metode Euler yang dimodifikasi untuk mencari y(x0 + h) :
       3


          y(x0 + h) ≅ y0 + hf(x0 + 2 h, y(x0 + 2 h))
                                   1           1


                    ≅ y0 + hf(x0 + 2 h, y0 + 2 hf )
                                   1         1
                                                 1

                    ≅ y0 + hf2                         f = f(x0 + h, y0 + hf )
                                                        3                   2




Jadi, menurut metode Runge-Kutta orde 4:

                         y(x0 + h) = y0 + 6 h(f + 2f + 2f + f )
                                          1
                                               0    1    2   3



dengan:


               f = f(x0 , y0 )
                0                              f = f(x0 + 2 h, y0 + 2 hf )
                                                2
                                                          1         1
                                                                        1

               f = f(x0 + 2 h, y0 + 2 hf0 )
                1
                          1         1
                                               f = f(x0 + h, y0 + hf )
                                                3                   2
Berangkat dengan quadrature Simpson, orang juga bisa memperoleh rumus
metode Runge-Kutta orde 3:
                             x0 +h

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



dengan:   f = f(x0 , y0 )
           0                  f = f(x0 + 2 h, y(x0 + 2 h))
                               1
                                         1           1
                                                                   f = f(x0 + h, y(x0 + h))
                                                                    2

y(x0 + 2 h) dicari dengan metode Euler dan y(x0 + h) dengan metode Euler
       1

yang dimodifikasi:

            y(x0 + 2 h) ≅ y0 + 2 hf
                   1           1
                                   0                  f = f(x0 + 2 h, y0 + 2 hf )
                                                       1
                                                                 1         1
                                                                               0

              y(x0 + h) ≅ y0 + hf1                    f = f(x0 + h, y0 + hf )
                                                       2                   1


Jadi, menurut metode Runge-Kutta orde 3:

                            y(x0 + h) = y0 + 6 h(f + 4f + f )
                                             1
                                                  0    1   2

                                     f = f(x0 , y0 )
                                      0

                                     f = f(x0 + 2 h, y0 + 2 hf )
                                      1
                                                 1        1
                                                              0

                                     f = f(x0 + h, y0 + hf )
                                      2                   1
                             PD Orde 2


                                  d2 y
Bentuk umum PD orde 2:        y''= 2 = f(x, y, y')
                                  dx

Diketahui:                    y(x0 ) = y0 , y'(x0 ) = y'0          y(x) = ?


Definisikan fungsi baru u:     u = y'                       y'= u(x, y)
                              u0 = y'0                      u'= f(x, y, u)




                                     Masalah PD orde 2
                                     berubah menjadi
                                     masalah PD orde 1.
Contoh penyelesaian dengan metode Euler yang lebih baik (improved):


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

          u(x0 + h) = u0 + 2 h(f + f )
                           1
                                 0    1           y(x0 + h) = y0 + 2 h(u0 + u1 )
                                                                   1


                 f = f(x0 , y0 , u0 )
                  0                                      u0 = y'0
                  f = f(x0 + h, y0 + hu0 , u1 )
                   1                                     u1 = u0 + hf0




Alur perhitungan:


     y0 , u0             f0               u1            f , y(x0 + h), u(x0 + h)
                                                         1




               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 + 6 h(f + 2f + 2f + f )
                     1
                           0    1    2      3    y(x0 + h) = y0 + 6 h(u0 + 2u1 + 2u2 + u3 )
                                                                  1


           f = f(x0 , y0 , u0 )
            0                                           u0 = y'0
            f = f(x0 + 2 h, y0 + 2 hu0 , u1 )
             1
                        1         1
                                                        u1 = u0 + 2 hf
                                                                  1
                                                                      0

           f = f(x0 + 2 h, y0 + 2 hu1 , u2 )
            2
                        1         1
                                                        u2 = u0 + 2 hf
                                                                  1
                                                                      1

            f = f(x0 + h, y0 + hu2 , u3 )
             3                                           u3 = u0 + hf2




Alur perhitungan:

y0 , u0       f0        u1        f1
                                            u2      f2
                                                             u3          f , y(x0 + h), u(x0 + h)
                                                                          3




                   x0 + h → x0 , u(x0 + h) → u0 , y(x0 + h) → y0
                        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:          (1) y'= f(x, y) = d(x) − e(x)y
                                           (2) y''= g(x, y, y') = a(x) − b(x)y − c(x)y'

Diketahui:                             x0 ≤ x ≤ xn
                                       y(x0 ) = y0                   y(x) = ?
                                       y(xn ) = yn
                                                                            xn − x0
Dicari yi = y(xi ) pada titik xi = x0 + ih (i = 1, ..., n − 1) dengan h =           .
                                                                               n
                            Metode Finite Differences
(1) y'+e(x)y = d(x)
                                                                        yi+1 − yi−1
(2) y''+c(x)y'+b(x)y = a(x)                                         y' ≅
                                                                     i
                                                                            2h
                                                                        y − 2yi + yi−1
                      (1) y' +ei yi = di
                           i                                       y'' ≅ i+1
                                                                     i
                                                                               h2
                      (2) y'' +ci y' +bi yi = ai
                            i      i

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

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

                                   n
                                                                                     1       n      
sistem persamaan linear:           ∑a x  ij   j   = bi (i = 1, ..., n)   solusi: xi =  bi − ∑ aijxj 
                                                                                     aii            
                                   j=1                                                      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)    1          n        
                                         bi − ∑ aijxj  (i = 1, ..., n)
                                                     (0)
                         x  i      =                    
                                    aii       j≠i       
Demikian seterusnya berulang-ulang, nilai perkiraan pada langkah ke k diperoleh
dari nilai perkiraan pada langkah ke k-1:

                                    1       n          
                        x (k)
                                   =  bi − ∑ aijxj -1)  (i = 1, ..., n)
                                                  (k

                                    aii                
                         i
                                           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) < ε,                    ε = bilangan kecil
                            xi
                                 Iterasi Gauss-Siedel


                                                   (k)    1                                   
                                                               bi − ∑ aijxj -1) − ∑ aijxj -1) 
                                                                           (k            (k
Rumus iterasi Jacobi dapat ditulis:              xi      =                                    
                                                          aii       j<i           j>i         

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

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

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 ( xj<i ), sehingga kesalahan
(ketidakakuratan) pada xj>i ( xj<i ) terakumulasi pada xi . Dengan
kata lain, terjadi perambatan kesalahan.


Pada metode iterasi tidak terdapat perambatan kesalahan seperti
itu. Semua elemen x dilihat secara sama. Pada tiap langkah
dilakukan pemeriksaan konvergensi untuk semua elemen x. Jadi,
untuk tiap elemen x terdapat kesempatan yang sama untuk
mencapai keakuratan yang diinginkan.


Namun, pada metode iterasi ada keharusan menentukan nilai awal,
yang bisa saja sulit dilakukan atau menimbulkan masalah, misalnya
membuat iterasi terlalu lama mencapai konvergensi.
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
                                                           cih                           ch 
PD orde 2: (2) y''= a(x) − b(x)y − c(x)y'                 1 −    yi−1 − (2 − bih2 )yi +  1 + i  yi+1 ≅ aih2
                                                              2                              2 
Iterasi Jacobi:

                                 1
               (1) yi(k) ≅
                                2eih
                                      (
                                     2dih + yi(k-1) − yi(k1-1)
                                              -1        +          )
                                    1                  ch                ch              
               (2) yi(k) ≅                   − aih2 +  1 − i  yi(k-1) +  1 + i  yi(k1-1) 
                                (2 − bih2 ) 
                                                           2 
                                                                   -1
                                                                                2 
                                                                                       +      
                                                                                              
Iterasi Gauss-Siedel (contoh untuk i membesar, i = 1, …, n-1):

                                 1
                  (1) yi(k) ≅
                                2eih
                                       (
                                     2dih + yi(k) − yi(k1-1)
                                              -1      +        )
                                    1                  ch              ch              
                  (2) yi(k) ≅                − aih2 +  1 − i  yi(k) +  1 + i  yi(k1-1) 
                                (2 − bih2 ) 
                                                           2 
                                                                   -1
                                                                              2 
                                                                                     +      
                                                                                            
                                               (k)                  (k)
Catatan, sesuai syarat batas:                 y0 = y0              yn = 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 )

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

Gunakan metode finite differences:

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


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


                               (h = xi+1 − xi = yj+1 − yj )
                                                                                            Dicari
                                                                                          distribusi
Diperoleh:                        1
                                     [
             ψi, j = h2 π ρi, j + 4 ψi+1, j + ψi-1, j + ψi, j+1 + ψi, j-1   ]             spasial ψ .
                   i,j+1             Langkah:
                                       1. Buat grid pada bidang xy, dengan jarak
                                          terdekat antar titik h.
    i-1,j             i,j   i+1,j      2. (Dianggap nilai pada batas-batas bidang
                                          xy diketahui.)
                   i,j-1                  Hitung dengan rumus :
    y
                                                               1
                                                                         [
                                          ψi, j = h2 π ρi, j + 4 ψ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
h                                        yang belum diketahui, gunakan nilai
                                    x    tebakan.
        h                              4. Ulangi langkah 2 – 3 sampai dicapai
                                          kestabilan untuk nilai ψi, j di semua titik:
                                                   (iterasi sebelum)             (iterasi berikutnya)
            ψ0,0                           ψi, j                       − ψi, j                          < ε (bilangan kecil)
                               Persamaan Differensial Parabolik


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

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

Gunakan metode finite differences:

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


             ∂                    ψ(xi , tj+1 ) − ψ (xi , tj ) ψi, j+1 − ψi, j
                ψ(x, t)         ≈                             =
             ∂t         xi , tj               ht                     ht
                                                                                              Untuk tiap
                                   (hx = xi+1 − xi , ht = tj+1 − tj )                        posisi dicari
                                                                                              perubahan ψ
                                                 Γht                                           terhadap
Diperoleh: ψi, j+1 = 4Γ ht π ρi, j + ψi, j +      2
                                                 hx
                                                       (
                                                     ψi+1, j − 2ψi, j + ψi-1, j   )             waktu.
                                      Langkah:
                  i,j+1
                                       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
         i-1,j    i,j     i+1,j           daerah x diketahui.)
                                          Hitung dengan rumus :

     t                                                                        Γht
                                          ψi, j+1 = 4Γ ht π ρi, j + ψi, j +    2
                                                                              hx
                                                                                   (
                                                                                  ψi+1, j − 2ψi, j + ψi-1, j   )
                                          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                                        Jika ρ dan nilai pada batas-batas daerah x
                                  x       tetap (tidak bergantung waktu), maka akan
                                          tercapai suatu waktu t, bahwa ψi, j+1 tidak
          hx                              berubah lagi (atau berubah hanya sedikit,
                                          sehingga dapat diabaikan):

           ψ0,0                            ψi, j+2 − ψi, j+1 < ε (bilangan kecil)
                            Persamaan Differensial Hiperbolik


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

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

Gunakan metode finite differences:

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


   ∂2                  ψ(xi , tj+1 ) − 2ψ (xi , tj ) + ψ(xi , tj-1 ) ψi, j+1 − 2ψi, j + ψi, j-1
      ψ(x, t)        ≈                                              =
  ∂t2         xi ,tj
                                          ht2                                    2
                                                                                ht               Untuk tiap
                       (hx = xi+1 − xi , ht = tj+1 − tj )                                       posisi dicari
                                                                                                 perubahan ψ
                                                                                                  terhadap
                                                             2 2
                                                          ch                                       waktu.
Diperoleh:
                              2

                                                            hx
                                                                   (
             ψi, j+1 = 4c2ht π ρi, j + 2ψi, j − ψi, j-1 + 2t ψi+1, j − 2ψi, j + ψi-1, j    )
Untuk j = 0 diperoleh ψi,1 sebagai berikut:
                                                            c2ht2
                             2 2
                    ψi,1 = 4c h π ρi,0 + 2ψi,0 − ψi,-1
                               t                           + 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
Maka gunakan:    ψ(x, t)        ≈                           =             , shg: ψi,-1 = ψi,1 − 2biht
              ∂t         xi ,t0              2ht                 2ht


Dengan demikian:

                                                                     2
                                                                  c2ht
       •    untuk j = 0:
                                        2 2
                               ψi,1 = 2c h π ρi,0 + biht + ψi,0 +
                                          t                          2
                                                                       (ψi+1,0 − 2ψi,0 + ψi-1,0 )
                                                                  2hx

                                                                       c2ht2

       •    untuk j > 0: ψi, j+1
                                         2 2

                                                                        hx
                                                                             (
                                    = 4c h π ρi, j + 2ψi, j − ψi, j-1 + 2 ψi+1, j − 2ψi, j + ψi-1, j
                                           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    2. (Dianggap nilai awal dan nilai pada batas-batas
                   i,j+1               daerah x diketahui.)
                                        Hitung dengan rumus :
                                                                              2
                                                                           c2ht
                                                     2 2
                                        ψi,1 = 2c h π ρi,0 + biht + ψi,0 +
                                                       t                      2
                                                                                (ψi+1,0 − 2ψi,0 + ψi-1,0 )
           i-1,j   i,j     i+1,j                                           2hx

                                                                                     c2ht2
       t           i,j-1                ψi, j+1        2 2
                                                  = 4c h π ρi, j + 2ψi, j − ψi, j-1 + 2 ψi+1, j − 2ψi, j + ψi-1, j
                                                         t
                                                                                      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.


ht
                                    x

ψ0,0        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:
         ⋱ ⋱                    ⋮       ⋮ 
                                             
         ⋱ ⋱       ⋱            yi−1    yi−1 
           f      − 2f     f                               k=?
                                 yi  = k yi 
             i          i    i
         
           h2      h 2
                            h 2
                                                                y=?
                   ⋱       ⋱ ⋱ yi+1     yi+1 
                                         ⋮ 
                           ⋱ ⋱ ⋮ 
                                                
                    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
                                   T
                  orthogonal: u u = δij
                                  i j                       komplit: x = ∑ ciui
                                                                                        i=1


Bermula dengan sembarang vector x, dilakukan iterasi berikut:

                           Ax= y                             y→x


                                                                n                 n                       n
Untuk kali pertama:                       y   (1)
                                                    = A x = A∑ ciui = ∑ ciAui = ∑ ciλiui
                                                                i=1               i=1                 i=1

                                                                      n                 n                     n
Setelah m kali iterasi diperoleh:         y   (m)       m
                                                    =A x =A     m
                                                                    ∑cu = ∑cA u = ∑c λ
                                                                      i=1
                                                                            i i
                                                                                        i=1
                                                                                              i
                                                                                                  m
                                                                                                      i
                                                                                                              i=1
                                                                                                                       m
                                                                                                                       u
                                                                                                                    i i i
                                                                                  λi≠k
Anggap λk merupakan eigenvalue terbesar:                                               <1
                                                                                   λk
Maka, jika m besar (banyak iterasi):

                                                                  λim 
              y   (m)
                        = A x = c λ u + ∑ c λ u = λ  ck uk + ∑ ci m ui  ≅ ck λk uk
                             m            m
                                                    
                                                              m          m      m

                                                                   λk 
                                        k k k              i i i         k
                                        i ≠k                 i ≠k      

λk diperoleh dengan jalan:
                                             n                     n
                                                                                                        x T y (m)
  T
 x y   (m)        T
             =x A x≅c λ m           m
                                  k k       ∑cu
                                            i=1
                                                    T
                                                  i i ku ≅c λ m
                                                            k k    ∑cδ
                                                                   i=1
                                                                         i ik
                                                                                  2 m
                                                                                ≅c λ
                                                                                  k k              λk = T (m-1)
                                                                                                       x y

uk diperoleh dengan jalan:
                        2                                                                          y (m)
              y   (m)
                            =y   (m)T
                                        y   (m)
                                                   (
                                                  ≅ cλ   m 2 T
                                                       k k) u u ≅ (c λ )
                                                             k k
                                                                           m 2
                                                                         k k                uk ≅
                                                                                                   y (m)

      (m)
Jika λk merupakan nilai λk yang diperoleh setelah iterasi sebanyak m kali,
maka iterasi dihentikan setelah dicapai nilai yang konvergen:
                                             (m-
                                           λk 1)
                                        1 − (m) < ε,               ε = bilangan kecil
                                            λk
Untuk mencari eigenvalue terbesar kedua, hilangkan uk dari perhitungan. Jadi,
dipakai vector awal baru x’:
                    n                  n
x'= (A − λk ) x = ∑ ci (A − λk )ui = ∑ ci (λi − λk )ui = ∑ ci (λi − λk )ui = ∑ diui     (di ≡ ci (λi − λk ))
                   i=1                i=1                i ≠k                    i≠ k


Iterasi:                         A x' = y'                      y'→ x'



                                              y'(m) = Am x'= ∑ diλi ui
                                                                           m
Setelah m kali iterasi diperoleh:
                                                                  i ≠k


                                                                         λi≠k ≠l
Anggap λl merupakan eigenvalue terbesar kedua:                                   <1
                                                                           λl

sehingga setelah banyak iterasi:              y'(m) ≅ dlλlmul

                                      x'T y'(m)                    y'(m)
Memperoleh λl dan ul :            λl = T (m-1)            ul ≅
                                      x' y'                        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 ui = λi ui         A-1A ui = λi A-1ui                  A-1ui = λi-1ui

Bermula dengan sembarang vector x, dilakukan iterasi berikut:

                              A-1 x = y                         y→x


                                                                   n
Setelah m kali iterasi diperoleh:         y   (m)
                                                    = (A ) x = ∑ ci (λi-1 ) m ui
                                                         -1 m

                                                                  i=1


Jika λs eigenvalue terkecil, maka setelah banyak kali iterasi:                                       -1
                                                                                        y (m) ≅ cs (λs ) m us

                                     xT y (m-1)                                 y (m)
Jadi, λs diperoleh sebagai:      λs = T (m)               dan us :       us ≅
                                     x y                                        y (m)

				
DOCUMENT INFO