Pemilihan proyek

Topik ini menyeleksi proyek akhir dari tiga kandidat. Ketiganya adalah solver PDE, Monte Carlo, dan penghitung nilai eigen. Pemilihan dilakukan dengan kriteria yang ditetapkan lebih dahulu, sehingga keputusan pemilihan dapat dinilai ulang dengan dasar yang sama. Hasil pemilihan ini adalah poisson1d, yaitu penyelesai persamaan Poisson satu dimensi, dan kandidat tersebut lalu dirancang hingga siap ditulis. Topik ini kemudian menetapkan tiga hal bagi poisson1d: spesifikasi matematisnya, arsitektur berkasnya, dan modul inti poisson_mod.f90. Terakhir, modul poisson_mod.f90 langsung diuji dengan program kecil yang dikompilasi dan dijalankan sungguhan.

Tiga kandidat proyek akhir#

Ketiga kandidat mewakili kelas pekerjaan numerik yang berbeda. Solver PDE mengubah persamaan diferensial menjadi sistem persamaan linear yang diselesaikan dengan algoritma eksak, Monte Carlo mengestimasi besaran melalui sampel acak, sedangkan penghitung nilai eigen mengandalkan sifat aljabar matriks simetris. Tabel berikut membandingkan ketiganya menurut cara memverifikasi hasil tiap kandidat, cakupan materi jalur yang tersentuh oleh kandidat tersebut, dan risiko pengerjaannya.

Keterverifikasian menjadi pembeda yang paling menentukan. Monte Carlo menghasilkan error yang menyusut lambat seiring banyaknya sampel, dengan laju $O(N^{-1/2})$, sehingga klaim akurasi hanya dapat dibuat setelah banyak eksperimen dilakukan. Penghitung nilai eigen mengharuskan penulis memahami urutan Sturm pada matriks tridiagonal, yang logikanya tidak dibahas pada jalur ini dan mudah salah tanpa alat pembanding. Solver PDE dengan metode beda hingga justru memiliki jalur verifikasi yang lurus, yaitu memilih persamaan yang solusi eksaknya diketahui, lalu mengukur error dan laju penyusutannya.

Kriteria pemilihan#

Kriteria yang dipakai ditetapkan lebih dahulu, sebelum satu kandidat pun dinilai.

• Keterverifikasian: program harus dapat dibandingkan dengan solusi yang diketahui secara eksak, karena kebenaran program numerik tidak dapat dinilai dari keberhasilan kompilasi semata.
• Cakupan materi: proyek harus memanfaatkan array, tipe turunan, modul, prosedur, keluaran berkas, dan argumen baris perintah.
• Kemandirian dependensi: proyek dibangun hanya dengan gfortran dan fpm tanpa pustaka tambahan, sehingga build-nya dapat direproduksi di mana pun.
• Ukuran: seluruh kode sumber berada di bawah tiga ratus baris, sehingga setiap bagiannya dapat dibaca tuntas.

Solver Poisson satu dimensi memenuhi keempat kriteria tersebut. Monte Carlo memenuhi kemandirian dependensi tetapi lemah pada cakupan materi, sedangkan penghitung nilai eigen melampaui ukuran yang wajar untuk proyek penutup karena teorinya sendiri perlu dibangun terlebih dahulu. Dengan demikian, proyek akhir yang dipilih adalah poisson1d.

Desain solver Poisson 1D#

Masalah yang diselesaikan proyek ini adalah persamaan Poisson satu dimensi dengan syarat batas Dirichlet homogen.

Domain $[0, 1]$ dibagi seragam menjadi $n + 2$ titik, yaitu $x_i = ih$ dengan $h = 1/(n+1)$. Nilai $u_0$ dan $u_{n+1}$ diketahui dari syarat batas, sehingga yang dicari adalah $n$ nilai interior $u_1, \dots, u_n$. Turunan kedua dihampiri beda hingga pusat orde dua.

Substitusi ke persamaan menghasilkan sistem tridiagonal berukuran $n \times n$.

Matriksnya bernilai 2 pada diagonal utama dan -1 pada dua diagonal di sampingnya, sedangkan ruas kanannya adalah $d_i = h^2 f_i$. Sistem tridiagonal tidak perlu difaktorkan dengan dekomposisi LU penuh yang biayanya $O(n^3)$, karena algoritma Thomas menyelesaikannya dengan dua sapuan linier berbiaya $O(n)$. Keputusan desain ini membentuk tata letak kode: jaring disimpan sebagai tipe turunan, sistem dirakit sebagai tiga array diagonal, dan penyelesaian sistemnya dipisahkan menjadi prosedur tersendiri.

Modul inti poisson_mod#

Proyek disusun sebagai satu modul pustaka bernama poisson_mod ditambah program-program yang memanggilnya. Modul memuat parameter kind dp, konstanta pi, tipe turunan jaring, dan empat prosedur publik.

module poisson_mod
   use, intrinsic :: iso_fortran_env, only: dp => real64
   implicit none
   private
   public :: dp, pi, jaring, buat_jaring, selesaikan, norma_inf, norma_l2

   real(dp), parameter :: pi = acos(-1.0_dp)

   ! tipe turunan penyimpan jaring (grid) seragam satu dimensi
   type :: jaring
      integer :: n                    ! jumlah titik interior
      real(dp) :: h                   ! jarak antartitik
      real(dp), allocatable :: x(:)   ! koordinat titik, ukuran n + 2
   end type jaring

contains

   ! membentuk jaring seragam [0, 1] dengan n titik interior
   function buat_jaring(n) result(g)
      integer, intent(in) :: n
      type(jaring) :: g
      integer :: i

      g%n = n
      g%h = 1.0_dp/real(n + 1, dp)
      allocate (g%x(0:n + 1))
      do i = 0, n + 1
         g%x(i) = real(i, dp)*g%h
      end do
   end function buat_jaring

   ! menyelesaikan -u'' = f pada [0, 1] dengan syarat batas u(0) = u(1) = 0;
   ! f memuat nilai ruas kanan pada setiap titik jaring (ukuran n + 2)
   subroutine selesaikan(g, f, u)
      type(jaring), intent(in) :: g
      real(dp), intent(in) :: f(0:)
      real(dp), allocatable, intent(out) :: u(:)
      real(dp), allocatable :: a(:), b(:), c(:), d(:)
      integer :: i, n

      n = g%n
      allocate (u(0:n + 1))
      u(0) = 0.0_dp
      u(n + 1) = 0.0_dp
      allocate (a(n), b(n), c(n), d(n))
      ! sistem diskret: -u_{i-1} + 2u_i - u_{i+1} = h^2 f(x_i)
      a = -1.0_dp
      b = 2.0_dp
      c = -1.0_dp
      do i = 1, n
         d(i) = g%h**2*f(i)
      end do
      call thomas(a, b, c, d, u(1:n))
   end subroutine selesaikan

   ! algoritma Thomas untuk sistem tridiagonal, biaya O(n)
   subroutine thomas(a, b, c, d, x)
      real(dp), intent(in) :: a(:), b(:), c(:), d(:)
      real(dp), intent(out) :: x(:)
      real(dp), allocatable :: cp(:), dq(:)
      real(dp) :: m
      integer :: i, n

      n = size(b)
      allocate (cp(n), dq(n))
      cp(1) = c(1)/b(1)
      dq(1) = d(1)/b(1)
      do i = 2, n
         m = b(i) - a(i)*cp(i - 1)
         cp(i) = c(i)/m
         dq(i) = (d(i) - a(i)*dq(i - 1))/m
      end do
      x(n) = dq(n)
      do i = n - 1, 1, -1
         x(i) = dq(i) - cp(i)*x(i + 1)
      end do
   end subroutine thomas

   ! norma maksimum (L-inf) dari selisih dua vektor
   function norma_inf(u, v) result(e)
      real(dp), intent(in) :: u(:), v(:)
      real(dp) :: e

      e = maxval(abs(u - v))
   end function norma_inf

   ! norma L2 diskrit dari selisih dua vektor pada jaring berjarak h
   function norma_l2(u, v, h) result(e)
      real(dp), intent(in) :: u(:), v(:), h
      real(dp) :: e

      e = sqrt(h*sum((u - v)**2))
   end function norma_l2

end module poisson_mod

Tipe jaring menyimpan jumlah titik interior $n$, jarak antartitik $h$, dan koordinat $x$ berukuran $n + 2$. Fungsi buat_jaring membentuk jaring seragam tersebut dan mengembalikannya sebagai satu objek, sehingga pemanggil tidak menyalin tiga variabel terpisah. Prosedur selesaikan menerima nilai ruas kanan $f$ pada setiap titik jaring, merakit tiga diagonal beserta ruas kanannya dengan faktor $h^2$, lalu memanggil thomas. Dua fungsi terakhir dipakai pada verifikasi: norma_inf menghitung selisih maksimum antartitik, sedangkan norma_l2 menghitung norma L2 diskrit, yaitu akar kuadrat dari jumlah kuadrat selisih yang ditimbang jarak jaring $h$.

Uji coba pertama#

Program uji berikut menyelesaikan masalah dengan $n = 4$ dan mencetak hasil numeriknya berdampingan dengan solusi eksaknya. Ruas kanannya diisi sebagai operasi array utuh, yaitu f = pi**2*sin(pi*g%x), tanpa perulangan eksplisit.

program uji_modul
   use poisson_mod, only: dp, pi, jaring, buat_jaring, selesaikan
   implicit none
   type(jaring) :: g
   real(dp), allocatable :: f(:), u(:)
   integer :: i

   g = buat_jaring(4)
   f = pi**2*sin(pi*g%x)
   call selesaikan(g, f, u)
   print '(A)', '    x    u numerik     u eksak'
   do i = 0, g%n + 1
      print '(F6.3, 2F13.6)', g%x(i), u(i), sin(pi*g%x(i))
   end do
end program uji_modul

# kompilasi program uji bersama modulnya
gfortran -Wall -Wextra poisson_mod.f90 uji_modul.f90 -o uji_modul
# jalankan program uji
./uji_modul

    x    u numerik     u eksak
 0.000     0.000000     0.000000
 0.200     0.607510     0.587785
 0.400     0.982972     0.951057
 0.600     0.982972     0.951057
 0.800     0.607510     0.587785
 1.000     0.000000     0.000000

Kompilasi selesai tanpa satu pun peringatan meskipun opsi -Wall -Wextra diaktifkan. Keluarannya menunjukkan perilaku yang diharapkan dari solusi kasar dengan $n = 4$: kedua ujung bernilai tepat nol sesuai syarat batas, solusi numeriknya simetris terhadap titik tengah domain, dan nilai $u$ numeriknya berada dekat tetapi tidak sama dengan solusi eksak. Selisih pada $x = 0{,}2$ adalah sekitar 0,02, dan selisih sebesar itu wajar karena jaringnya hanya memuat empat titik interior. Penyusutan selisih tersebut diukur secara sistematis pada verifikasi kebenaran.

Ringkasan#

• Dari tiga kandidat, solver Poisson satu dimensi dipilih karena paling dapat diverifikasi, paling padat cakupan materi, bebas dependensi, dan berukuran kecil.
• Masalah diskret menjadi sistem tridiagonal dengan diagonal 2 dan -1, yang diselesaikan algoritma Thomas dengan biaya $O(n)$.
• Modul poisson_mod memuat tipe jaring, perakit sistem, algoritma Thomas, dan dua norma error, sedangkan program terpisah memanggilnya.
• Uji pertama dengan $n = 4$ berjalan tanpa peringatan dan menghasilkan solusi yang simetris serta memenuhi syarat batas.

Referensi#

• Learn: tutorial resmi Fortran — fortran-lang — Rujukan fitur bahasa modern yang dipakai pada seluruh kode proyek ini.
• The GNU Fortran Compiler — GNU Project — Manual gfortran yang menjadi rujukan opsi peringatan -Wall -Wextra pada kompilasi proyek.
• Tridiagonal matrix algorithm — Wikipedia — Penjelasan algoritma Thomas beserta turunan sapuan maju dan mundurnya.