11

Ketika seseorang ingin menghitung turunan numerik, metode yang disajikan oleh Bengt Fornberg di sini (dan dilaporkan di sini ) sangat mudah (baik tepat dan mudah diimplementasikan). Sebagai kertas asli tanggal dari 1988, saya ingin tahu apakah ada alternatif yang lebih baik hari ini (sebagai (atau hampir sama) sederhana dan lebih tepat)?

finite-difference precision Vincent
sumber

1

Sulit dikatakan tanpa mengetahui apa yang ingin Anda bedakan. Sudahkah Anda mempertimbangkan diferensiasi otomatis ?

Biswajit Banerjee

@BiswajitBanerjee: Untuk koefisien perbedaan hingga, diferensiasi otomatis tidak berlaku.

Geoff Oxberry

@ GeoffOxberry: Saya merujuk pada masalah fisik yang sebenarnya, yaitu, bagian sains dari "ilmu komputasi".

Biswajit Banerjee

@Vincent: Apakah Anda mencoba membedakan fungsi, atau tabel? Jika data tabel, apakah data tabel berisik? Apakah Anda mencoba untuk mendiskritisasi PDE?

user14717

9

Gambaran

Pertanyaan bagus. Ada sebuah makalah yang berjudul "Meningkatkan akurasi metode diferensiasi matriks untuk titik kolokasi sewenang-wenang" oleh R. Baltensperger. Ini bukan masalah besar menurut saya, tetapi ada benarnya (yang sudah diketahui sebelum penampilan pada tahun 2000): ini menekankan pentingnya representasi akurat dari fakta bahwa turunan dari fungsi konstan $f(x)=1$ harus menjadi nol (ini berlaku persis dalam arti matematika, tetapi tidak harus dalam representasi numerik).

Sangat mudah untuk melihat bahwa ini membutuhkan jumlah baris dari matriks turunan ke-n $D^{(n)}$ menjadi nol. Hal ini umum untuk menegakkan kendala ini dengan menyesuaikan entri diagonal, yaitu dengan menetapkan

\begin{matrix} (1) & D_{j j}^{(n)} := - \sum_{\begin{matrix} i = 1 \\ i \neq j \end{matrix}}^{N} D_{i j} . \end{matrix}

$\tag 1 D^{(n)}_{jj} := -\sum_{\substack{i=1\\i \neq j }}^N D_{ij}\,.$ Jelas bahwa fitur ini tidak berlaku ketika bekerja pada komputer karena kesalahan pembulatan dalam perhitungan floating point. Apa yang lebih mengejutkan adalah bahwa kesalahan ini bahkan lebih parah ketika menggunakan rumus analitik untuk matriks turunan (yang tersedia untuk banyak titik kolokasi klasik, misalnya Gauss-Lobatto).

Sekarang, makalah (dan referensi di dalamnya) menyatakan bahwa kesalahan turunan adalah dalam urutan penyimpangan dari jumlah baris dari nol. Karena itu tujuannya adalah untuk membuat angka-angka ini sekecil mungkin.

Tes numerik

Poin baiknya adalah prosedur Fornberg tampaknya cukup bagus dalam hal ini. Pada gambar di bawah ini saya telah membandingkan perilaku yang tepat, yaitu analitik, matriks turunan pertama dan yang diturunkan oleh algoritma Fornberg, untuk jumlah titik kolokasi Chebyshev-Lobatto yang bervariasi.

Sekali lagi, meyakini pernyataan dalam makalah yang dikutip, ini menyiratkan bahwa algoritma Fornberg akan menghasilkan hasil yang lebih akurat untuk turunannya.

Untuk membuktikannya, saya akan menggunakan fungsi yang sama seperti di kertas,

\begin{matrix} (2) & f (x) = \frac{1}{1 + x^{2}} . \end{matrix}

$\tag 2 f(x) = \frac{1}{1+x^2}\,.$

\begin{matrix} (3) & E_{n} = max_{i \in {0, \dots, n}} | f^{'} (x_{i}) - \sum_{j = 1}^{n} D_{i j} f (x_{j}) | . \end{matrix}

$\tag 3 E_n = \max_{i\in \{0,\ldots,n\}} \bigg| f^\prime(x_i) - \sum_{j=1}^n D_{ij} f(x_j) \bigg|\,.$

\begin{matrix} (4) & {\tilde{D}}_{j j} = D_{j j} - (\sum_{i = 1}^{n} D_{j i}), for all j . \end{matrix}

$\tag 4 \tilde D_{jj} = D_{jj} - \left(\sum_{i=1}^n D_{ji} \right), \qquad \text{for all}\ j\,.$

Kesimpulan

Kesimpulannya, metode Fornberg tampaknya cukup akurat, dalam kasus bahkan sekitar 3 urutan besarnya lebih akurat daripada hasil dari rumus analitik. Ini harus cukup akurat untuk sebagian besar aplikasi. Selain itu, ini luar biasa karena Fornberg tampaknya tidak secara eksplisit memasukkan fakta ini dalam metodenya (setidaknya tidak disebutkan dalam dua makalah Fornberg). $N=512$

Urutan besar lainnya dapat diperoleh untuk contoh ini melalui penyertaan langsung Persamaan (4). Karena ini adalah pendekatan yang cukup sederhana dan hanya diterapkan sekali untuk setiap turunan, saya tidak melihat alasan untuk tidak menggunakannya.

Metode dari makalah Baltensperger - yang menggunakan pendekatan yang lebih canggih untuk mengevaluasi jumlah dalam Persamaan (1) untuk mengurangi kesalahan pembulatan - menghasilkan sekitar urutan yang sama besarnya untuk kesalahan. Jadi, setidaknya untuk contoh ini, kira-kira setara dengan metode "Adjusted Fornberg" di atas.

Davidhigh
sumber

4

Dengan asumsi Anda mencoba untuk membedakan implementasi numerik dari fungsi kontinu, ada sejumlah besar metode:

1) Diferensiasi otomatis. Metode yang paling akurat dan umum. Nyeri untuk kode, membutuhkan operator kelebihan dan pencarian tergantung argumen. Memberi beban pada pengguna untuk memahami konsep-konsep ini. Juga berjuang dengan singularitas yang dapat dilepas, seperti membedakan sinc pada . $x= 0$

2) Transformasi Chebyshev. Proyeksikan fungsi Anda ke rentang polinomial Chebyshev dan bedakan perulangan tiga istilah. Sangat cepat, sangat akurat. Tetapi mengharuskan Anda memiliki domain minat yang didukung secara kompak; di luar domain yang dipilih , perulangan istilah tiga tidak stabil. $[a,b]$

3) Perbedaan terbatas. Underrated dalam 1D; lihat Tip dan Trik Nick Higham dalam Komputasi Numerik . Idenya adalah bahwa jika Anda menyeimbangkan kesalahan pemotongan dan kesalahan roundoff maka Anda tidak perlu memilih stepsize; itu dapat dipilih secara otomatis. Dalam Boost, ide ini digunakan untuk memulihkan (secara default) 6/7 digit yang benar untuk tipe tersebut. (Higham hanya menunjukkan ide untuk case yang lebih sederhana dari 1/2 digit yang benar, tetapi idenya mudah diperpanjang.) Koefisiennya berasal dari tabel Fornberg yang sama, tetapi ukuran langkah dipilih dengan asumsi bahwa fungsi tersebut dapat dievaluasi menjadi 1ULP ketepatan. Kerugiannya adalah itu membutuhkan 2 evaluasi fungsi untuk memulihkan setengah digit dari tipe, 4 untuk memulihkan 3/4 digit, begitu seterusnya. Dalam 1D, bukan transaksi yang buruk. Dalam dimensi yang lebih tinggi, ini adalah bencana besar.

4) Langkah turunan kompleks. Gunakan . Ambil untuk menjadi unit roundoff dan ini akan pulih hampir setiap bit dengan benar. Namun, itu agak curang, karena umumnya lebih sulit untuk mengimplementasikan suatu fungsi dalam bidang kompleks daripada memberikan kode turunan nyata. Masih ide yang keren dan bermanfaat dalam kondisi tertentu. $f'(x) \approx \Im(f(x+ih))$ $h$

pengguna14717
sumber

1

Saya tidak mengetahui ada yang memperbaiki algoritma Fornberg (lihat juga makalahnya yang sedikit lebih baru ). Sebagai tambahan, bagi saya tampaknya melihat algoritme sebagai cara untuk menghitung turunan numerik tidak benar. Semua yang dia lakukan adalah mendapatkan algoritma yang efisien untuk menghitung bobot untuk metode beda hingga. Keuntungan dari metodenya adalah memberi Anda bobot untuk semua turunan hingga turunan yang diinginkan dalam sekali jalan.

Brian Zatapatique
sumber

Saya tidak setuju dengan pernyataan "melihat algoritme sebagai cara menghitung turunan numerik tidak benar". Setelah Anda memiliki bobot dievaluasi untuk titik , Anda dapat langsung menghitung turunan numerik fungsi setiap melalui .

w_{i}

$w_i$

z

$z$

f (x)

$f(x)$

f^{'} (z) = \sum_{i} w_{i} f (x_{i})

$f^\prime(z) = \sum_i w_i f(x_i)$

davidhigh

@davidhigh: Jika Anda membaca makalah Fornberg, mereka berbicara tentang perhitungan bobot untuk perkiraan perbedaan hingga, bukan tentang menghitung perkiraan itu sendiri. Tentu saja Anda dapat menggunakan algoritme untuk menghitung turunannya juga, tetapi lebih masuk akal untuk menghitung bobot, menyimpannya, dan kemudian menggunakannya berulang kali untuk menghitung perkiraan turunannya. Kalau tidak, Anda menghitung bobot berulang kali ketika tidak perlu.

Brian Zatapatique

1

Skema yang lebih sederhana

Selain jawaban saya yang lain yang lebih lanjut tentang perluasan metode Fornberg, saya akan membahas di sini pertanyaan untuk alternatif yang lebih sederhana.

Untuk ini saya membuat sketsa skema alternatif yang menghasilkan koefisien turunan interpolasi Lagrangian lebih langsung. Implementasinya hanya membutuhkan beberapa baris kode, berfungsi untuk grid sewenang-wenang, dan menurut eksperimen pertama saya, seakurat Fornberg.

Dasar implementasi adalah turunan langkah imajiner mana adalah variabel dalam urutan ketepatan mesin. Derivatif-langkah imajiner diketahui secara stabil menghasilkan nilai-nilai turunan dan tidak mengalami ketidakstabilan numerik dari implementasi beda hingga dengan .

f^{'} (x) = \frac{1}{h} Aku (f (x + saya h)),

$f'(x) \ = \ \frac{1}{h}\,\text{Im}\big(f(x+ih)\big)\,,$

h

$h$

h \to 0

$h\rightarrow 0$

Bahan kedua adalah polinomial interpolasi Lagrange pada grid dievaluasi dengan salah satu bentuk barycentric, misalnya mana Untuk menggunakan turunan langkah-kompleks, kita harus memastikan formula ini juga berfungsi untuk kompleks argumen . Selain itu, untuk fungsi yang diberikan f (x) dan vektor koefisien , kami menunjukkan polinomial interpolasi melalui oleh $L_i(x)$ $\{x_1,\ldots, x_N\}$

{L.}_{saya} (z) = {\begin{cases} 1 & jika z = x_{saya} \\ \frac{\frac{μ_{saya}}{z - x_{k}}}{\sum_{k} \frac{μ_{k}}{z - x_{k}}} & jika tidak . \end{cases}

$L_i(z) \ = \ \begin{cases}1 & \text{if }z = x_i\\\ \frac{\frac {\mu_i}{z-x_k}}{ \sum_k\frac {\mu_k}{z-x_k} } & \text{otherwise}\,. \end{cases}$

μ_{i} = \frac{1}{\prod_{k \neq i} (x_{i} - x_{k})}

$\mu_i = \frac{1}{\prod_{k\neq i} (x_i - x_k)}$

z

$z$

f_{i} = f (x_{i})

$f_i = f(x_i)$

(x_{1}, f_{i})

$(x_1,f_i)$

P (x; f) = \sum_{saya = 1}^{N} f_{saya} {L.}_{saya} (x) .

$P(x;f) = \sum_{i=1}^N f_i L_i(x)\,.$

Algoritma

Algoritma ini dibuat sketsa sebagai berikut. Ini memiliki input dan output parameter yang sama dengan Fornberg, tetapi jauh lebih pintar.

Memasukkan:

$x$ : kisi dengan N titik jaringan yang berbeda
$ord$ : pesanan turunan
z: titik di mana turunan akan dievaluasi
Kemungkinan: suatu fungsi atau nilai-nilai funcitonnya pada titik-titik (hanya diperlukan untuk varian keluaran 2.) $f(x)$ $f_i$

Inisialisasi

Inisialisasi polinomial Lagrange melalui interpolasi barycentric. $L$
Inisialisasi larik -tingkat , untuk $N\times N$ $D^{(o)}$ $o=0, \ldots, ord$
Set , yaitu matriks satuan dimensi $D^{(0)} := I_N$ $N\times N$
Set $o := 0$

Algoritma

Sementara : $o < ord$

Hitung oleh langkah turunan kompleks untuk semua dan . Di sini, menunjukkan baris -th . $D^{(o+1)}_{ik} = CSD \bigg( P\big(x_k; D^{(o)}_{i}\big) \bigg)$ $CSD$ $i$ $k$
$D^{(o)}_{i}$ $i$ $D^{(o)}$
Set o = o + 1;

Putuskan apa yang akan Output :

Vektor dari koefisien Perbedaan Hingga pada titik , di mana . Inilah yang dilakukan Fornberg. $d^{(ord)}$ $z$ $d_i = P(z;D^{(ord)}_i)$
Fungsi interpolasi ke turunan dari pesanan . Untuk ini, Anda harus memasukkan fungsi resp. nilai fungsi di ke algoritma. $p^{(ord)}(x) = \sum_{i=1}^N f(x_i) P(x;D^{(ord)}_i)$ $f^{(ord)}(x)$ $ord$ $f_i$ $x_i$
Sebuah fungsi-meta yang mengembalikan fungsi interpolasi varian 2., tetapi untuk fungsi sewenang-wenang yang akan diinterpolasi pada titik-titik grid. $\text{diff}(\text{int } ord, \text{function } g)$ $g$

Secara pribadi, saya suka varian 3. yang paling.

Analisis algoritma

Seperti Fornberg, algoritma ini adalah . Saya akan memposting hasil yang lebih empiris mengenai akurasi, stabilitas, dll. Jika saya menemukan waktu. $\mathcal O(ord \cdot N^2)$

Davidhigh
sumber

0

Untuk meningkatkan ketepatan diferensiasi numerik, lakukan hal berikut:

1) Pilih metode "standar" presisi tinggi favorit Anda berdasarkan beberapa ukuran langkah h .

2) Hitung nilai derivatif dengan metode yang dipilih dalam 1) berkali-kali dengan ukuran langkah yang berbeda tetapi masuk akal h . Setiap kali Anda dapat memilih h sebagai angka acak dari interval (0,5 * H / 10, 1,5 * H / 10) di mana H adalah ukuran langkah yang tepat untuk metode yang Anda gunakan.

3) Rata-rata hasilnya.

Hasil Anda mungkin mendapatkan 2-3 urutan besarnya dalam kesalahan kesalahan mutlak. hasil tidak rata-rata.

https://arxiv.org/abs/1706.10219

F. Jatpil
sumber

3

Selamat datang di SciComp.SE! Jawaban ini akan lebih baik jika Anda meringkas metode ini secara singkat.

Christian Clason

2

Nama pengguna Anda juga menyarankan agar Anda menulis makalah itu. Harap baca panduan kami tentang promosi diri dan edit posting Anda sesuai dengan itu.

Wrzlprmft

Jujur saya pikir pemungutan suara negatif tidak adil jika jawaban saya memang menunjuk ke jawaban yang valid.

F. Jatpil

1

Saya juga ... jadi ambil +1 saya ... :-)