Verifikasi terhadap kebenaran

Program yang berhasil dikompilasi belum tentu benar, sehingga proyek numerik menuntut verifikasi terhadap solusi yang telah diketahui melalui estimasi error. Topik ini memverifikasi poisson1d dengan tiga alat, yaitu solusi terbuat (manufactured solution), dua norma error, dan studi konvergensi pada enam jaring yang makin halus. Modul poisson_mod.f90 dari topik pemilihan proyek disalin utuh ke direktori ini, sehingga yang diverifikasi adalah kode yang sama dengan yang sudah diuji secara kasar. Orde konvergensi yang diukur dari keluaran menjadi bukti kuantitatif bahwa diskretisasi orde dua diterapkan dengan benar.

Solusi terbuat#

Metode solusi terbuat membalik arah perancangan persamaan: solusi yang diinginkan dipilih lebih dahulu, kemudian ruas kanan yang menghasilkannya dihitung dari turunannya. Cara ini dipakai karena contoh Poisson yang dapat diselesaikan secara analitis terbatas pada bentuk $f$ tertentu, sedangkan metode ini menghadirkan pasangan $u$ dan $f$ yang konsisten secara eksak pada setiap titik jaring. Untuk proyek ini dipilih solusi eksak berikut.

Fungsi ini memenuhi kedua syarat batas karena $\sin(0) = \sin(\pi) = 0$. Turunan keduanya adalah $u''(x) = -\pi^2 \sin(\pi x)$, sehingga persamaan $-u'' = f$ terpenuhi apabila ruas kanannya berupa

Dalam kode, ruas kanan ini diisi sebagai satu operasi array utuh tanpa perulangan, yaitu f = pi**2*sin(pi*g%x), dengan konstanta pi yang dideklarasikan memakai kind dp dari modul. Nama solusi terbuat menegaskan sifat artifisial pasangan ini: solusi dan ruas kanan sama-sama dirakit untuk keperluan uji, sehingga tidak ada ketidakpastian data masukan yang menutupi error diskretisasi. Akibatnya, setiap selisih yang terukur pada topik ini murni berasal dari diskretisasi dan dapat dibandingkan dengan nilai eksaknya pada titik yang sama.

Norma error#

Selisih antara solusi numerik dan solusi eksak berupa vektor berukuran $n + 2$, sehingga perlu diringkas menjadi satu bilangan. Dua norma yang dipakai adalah norma maksimum dan norma L2 diskrit.

Norma maksimum menyatakan error terburuk pada satu titik, sedangkan norma L2 diskrit menaksir error pada keseluruhan domain dengan pembobotan lebar sel $h$. Keduanya saling melengkapi: satu lonjakan lokal langsung terlihat pada norma maksimum, sedangkan norma L2 mengikuti kecenderungan umum dan tidak bergeser karena satu titik saja. Pemilihan bentuk kuadrat berbobot $h$ pada norma L2 memiliki alasan teknis. Kuantitas di bawah akar, yaitu $h \sum_i e_i^2$ dengan $e_i = u_i - u(x_i)$, merupakan jumlah Riemann dari $\int_0^1 e(x)^2\,dx$, sehingga $e_2$ menghampiri norma kontinu $\left(\int_0^1 e(x)^2\,dx\right)^{1/2}$ ketika jaring dihaluskan. Berkat bobot itu pula nilai $e_2$ dari jaring yang berbeda tetap sebanding, sesuatu yang tidak berlaku bagi jumlah kuadrat error tanpa bobot karena jumlah itu membengkak seiring bertambahnya titik.

Jumlah pada norma L2 mencakup seluruh $n + 2$ titik termasuk kedua ujung, tetapi kontribusi ujung itu nol karena solusi numerik maupun eksaknya memenuhi syarat batas yang sama. Keduanya dinyatakan sebagai fungsi norma_inf dan norma_l2 di dalam modul, dan keduanya menuju nol apabila skema diskretisasinya konvergen. Nilai kedua norma inilah yang dibandingkan antarjaring pada studi konvergensi.

Studi konvergensi#

Program studi_konvergensi.f90 menyelesaikan masalah yang sama pada enam jaring yang makin halus, yaitu $n = 7, 15, 31, 63, 127, 255$. Daftar ini dipilih dengan sengaja: setiap kenaikan mengubah $n$ menjadi $2n + 1$, sehingga $h$ selalu tepat terbagi dua antara dua kasus berurutan. Orde konvergensi terukur dihitung dari dua error berurutan dengan rumus

Dalam kode, err_inf_lama menyimpan error kasus sebelumnya sehingga orde setiap baris terhitung di perulangan yang sama, dan pembagi $\ln(h_{i-1}/h_i)$ disederhanakan menjadi $\log 2$ karena daftar $n$ menjamin $h$ tepat terbagi dua.

program studi_konvergensi
   use poisson_mod, only: dp, pi, jaring, buat_jaring, selesaikan, norma_inf, norma_l2
   implicit none
   integer, parameter :: n_kasus = 6
   integer :: daftar_n(n_kasus) = [7, 15, 31, 63, 127, 255]
   type(jaring) :: g
   real(dp), allocatable :: f(:), u(:), u_eksak(:)
   real(dp) :: err_inf, err_l2, err_inf_lama
   integer :: kasus

   err_inf_lama = 0.0_dp
   print '(A)', '    n        h      error L-inf     error L2      orde'
   do kasus = 1, n_kasus
      g = buat_jaring(daftar_n(kasus))
      f = pi**2*sin(pi*g%x)
      call selesaikan(g, f, u)
      u_eksak = sin(pi*g%x)
      err_inf = norma_inf(u, u_eksak)
      err_l2 = norma_l2(u, u_eksak, g%h)
      if (kasus == 1) then
         print '(I6, F10.7, 2ES15.6, A)', g%n, g%h, err_inf, err_l2, '       -'
      else
         print '(I6, F10.7, 2ES15.6, F8.3)', g%n, g%h, err_inf, err_l2, &
              log(err_inf_lama/err_inf)/log(2.0_dp)
      end if
      err_inf_lama = err_inf
   end do
end program studi_konvergensi

# kompilasi program studi konvergensi bersama modulnya
gfortran -Wall -Wextra poisson_mod.f90 studi_konvergensi.f90 -o studi_konvergensi
# jalankan studi konvergensi
./studi_konvergensi

    n        h      error L-inf     error L2      orde
     7 0.1250000   1.295075E-02   9.157561E-03       -
    15 0.0625000   3.218964E-03   2.276152E-03   2.008
    31 0.0312500   8.035777E-04   5.682152E-04   2.002
    63 0.0156250   2.008218E-04   1.420025E-04   2.001
   127 0.0078125   5.020092E-05   3.549741E-05   2.000
   255 0.0039062   1.254995E-05   8.874152E-06   2.000

Beda hingga pusat orde dua memiliki sisa (truncation error) berorde $O(h^2)$, sehingga orde terukur harus mendekati 2,0 ketika jaring makin halus. Orde konvergensi yang terbaca dari kolom terakhir juga memiliki makna praktis: setiap kali $h$ terbagi dua, error diperkirakan menyusut sekitar empat kali, dan untuk menambah satu angka penting akurasi, $h$ cukup diperkecil sekitar $10^{1/2} \approx 3{,}16$ kali. Keluarannya memenuhi tuntutan teori itu: orde dimulai dari 2,008 pada jaring paling kasar, lalu menurun menuju 2,000 karena pengaruh suku orde lebih tinggi pada sisa diskretisasi semakin kecil. Error L-inf sendiri menyusut dari 1.295075E-02 pada $n = 7$ menjadi 1.254995E-05 pada $n = 255$.

Tabel ini mudah direproduksi pembaca. Cukup letakkan poisson_mod.f90 dan studi_konvergensi.f90 dalam satu direktori, jalankan perintah kompilasi dan eksekusi di atas, dan keenam baris tercetak seketika karena setiap kasus diselesaikan algoritma Thomas berbiaya $O(n)$. Baris pertama tidak memiliki orde karena belum ada kasus sebelumnya sebagai pembanding, sehingga kolom terakhirnya diisi tanda hubung. Salah satu baris juga dapat diperiksa dengan tangan: untuk $n = 15$, rasio error terhadap kasus sebelumnya adalah $1{,}295075 \times 10^{-2} / 3{,}218964 \times 10^{-3} \approx 4{,}023$, dan $\ln 4{,}023 / \ln 2 \approx 2{,}008$, persis seperti yang tercetak.

Penyusutan total sepanjang tabel dapat diperiksa dari kedua baris ujungnya. Error L-inf pada $n = 255$ menyisakan $1{,}295075 \times 10^{-2} / 1{,}254995 \times 10^{-5} \approx 1032$ bagian dari error awal, jadi error menyusut sekitar seribu kali lipat sepanjang studi. Teori memberikan dugaan untuk angka itu dari sisi $h$: jumlah selang bertambah dari 8 menjadi 256, sehingga $h$ terbagi dua tepat 32 kali, dan karena error berbanding lurus dengan $h^2$, error diperkirakan menyusut $32^2 = 1024$ kali. Rasio terukur 1032 bertepatan dengan dugaan tersebut; penyimpangannya yang hanya sekitar satu persen berasal dari suku orde lebih tinggi pada sisa diskretisasi yang masih berpengaruh pada jaring kasar. Ketika orde terukur menyimpang jauh dari nilai teoretisnya, dugaan pertama jatuh pada kesalahan kode atau pada cara mengukurnya, dan kedua kemungkinan itu dibahas pada bagian berikut.

Kesalahan umum#

• Error membesar saat jaring makin halus. Gejala ini biasanya berarti tanda atau skala sistem keliru, misalnya ruas kanan lupa dikalikan $h^2$ atau diagonal utama tidak bernilai 2. Perbaikannya dengan menyelesaikan kasus kecil $n = 1$ dan $n = 2$ secara manual, lalu mencocokkan hasilnya dengan sistem yang dirakit program.
• Orde dihitung dari jaring yang tidak bersarang. Rumus $p_i$ mengasumsikan $h_{i-1} = 2h_i$. Daftar $n$ yang asal, misalnya 10, 20, 40, membuat pembagi $\ln(h_{i-1}/h_i)$ tidak lagi sama dengan $\ln 2$, sehingga orde terukur tidak menuju 2 meskipun kodenya benar. Perbaikannya dengan mempertahankan pola $n$ dan $2n + 1$.
• Toleransi yang lebih ketat daripada round-off. Pada $h$ yang sangat kecil, error pembulatan aritmetika double precision menumpuk dan error total berhenti menyusut. Ambang kelulusan uji diset jauh di atas lantai tersebut, misalnya 1.0e-3 ketika error terukur pada $n = 127$ masih 5.020092E-05.
• Menilai benar dari satu titik. Satu nilai $u$ yang kebetulan cocok tidak membuktikan apa pun, karena kesalahan perakitan sering membiarkan sebagian titik tetap tepat. Verifikasi menilai seluruh jaring melalui norma, bukan satu entri.

Ringkasan#

• Solusi terbuat $u = \sin(\pi x)$ menghasilkan ruas kanan $f = \pi^2 \sin(\pi x)$ dan memberi pembanding eksak pada setiap titik jaring.
• Error diringkas dengan norma maksimum dan norma L2 diskrit, keduanya tersedia sebagai fungsi modul.
• Studi konvergensi pada enam jaring bersarang menghasilkan orde 2,008 hingga 2,000, sesuai teori beda hingga orde dua.
• Orde konvergensi hanya bermakna ketika jaring berurutan saling bersarang, yaitu $h$ tepat terbagi dua.

Referensi#

• Learn: tutorial resmi Fortran — fortran-lang — Rujukan fitur bahasa yang dipakai pada program studi konvergensi ini.
• The GNU Fortran Compiler — GNU Project — Manual gfortran untuk format keluaran ES dan intrinsik log yang dipakai.
• Finite difference method — Wikipedia — Uraian beda hingga pusat beserta orde sisa diskretisasinya.