Mencari urutan Runge-Kutta 8 di C / C ++

10

Saya ingin menggunakan metode urutan ke-8 Runge-Kutta (89) dalam aplikasi mekanika / astrodinamika surgawi, yang ditulis dalam C ++, menggunakan mesin Windows. Karena itu saya bertanya-tanya apakah ada yang tahu perpustakaan / implementasi yang baik yang didokumentasikan dan bebas untuk digunakan? Tidak apa-apa jika ditulis dalam C, selama tidak ada masalah kompilasi yang diharapkan.

Sejauh ini saya telah menemukan perpustakaan ini (mymathlib) . Tampaknya kode itu ok, tetapi saya belum menemukan informasi tentang lisensi.

Bisakah Anda membantu saya dengan mengungkapkan beberapa alternatif yang mungkin Anda ketahui dan akan sesuai dengan masalah saya?

EDIT:
Saya melihat bahwa sebenarnya tidak ada banyak kode sumber C / C ++ yang tersedia seperti yang saya harapkan. Oleh karena itu versi Matlab / Oktaf akan ok juga (masih harus bebas untuk digunakan).

c++ reference-request libraries c runge-kutta James C
sumber

Jawaban:

8

Kedua Perpustakaan GNU Ilmiah (GSL) (C) dan Meningkatkan Odeint (C ++) fitur-8 agar metode Runge-Kutta.

Keduanya adalah opensource, dan di bawah linux dan mac mereka harus langsung tersedia dari manajer paket. Di bawah windows, mungkin Anda akan lebih mudah menggunakan Boost daripada GSL.

GSL diterbitkan di bawah lisensi GPL, dan Boost Odeint di bawah lisensi Boost.

Sunting: Ok, Boost Odeint TIDAK memiliki metode Runge-Kutta 89, hanya 78, tetapi itu memang memberikan resep untuk membuat stepper Runge-Kutta yang sewenang-wenang.

Metode urutan ke-8 cukup tinggi, dan kemungkinan besar terlalu banyak untuk masalah Anda.

Prince-Dormand mengacu pada jenis tertentu Runge-Kutta, dan tidak terkait langsung dengan pesanan tetapi yang paling umum adalah 45. Matlabs ode45, yang merupakan algoritma ODE yang direkomendasikan adalah implementasi Prince-Dormand 45. Ini adalah algoritma yang sama seperti yang diterapkan dalam Boost Odeint Runge_Kutta_Dopri5 .

LKlevin
sumber
1
Terimakasih untuk jawaban. OK itu memalukan sekarang, saya telah melihat Boost Odeint bahkan sebelum bertanya di sini, dan hanya menemukan "runge_kutta_fehlberg78". Apakah ini hal yang benar? Sebenarnya saya tidak tahu perbedaan antara methonds ketika digunakan dalam praktek, tetapi saya sedang mencari RK89 (disebut juga Dormand-Prince ketika saya mencari di internet). Bisakah Anda berkomentar atau memperluas jawaban Anda mengenai masalah ini? Terima kasih.
James C
Pos yang diperbarui untuk menjawab pertanyaan Anda. Prince-Dormand 45 kemungkinan besar akan menyelesaikan masalah Anda dengan baik.
LKlevin
15

Jika Anda melakukan mekanika langit dalam skala waktu yang lama, menggunakan integrator Runge-Kutta klasik tidak akan menghemat energi. Dalam hal ini, menggunakan integrator symplectic mungkin akan lebih baik. Boost.odeint juga mengimplementasikan skema Runge-Kutta symplectic tingkat 4 yang akan bekerja lebih baik untuk interval waktu yang lama. GSL tidak menerapkan metode symplectic, sejauh yang saya tahu.

Geoff Oxberry
sumber
Terimakasih untuk jawaban. Akankah Runge-Kutta symplectic urutan ke-4 memberikan hasil yang lebih baik daripada RKF78, jika digunakan dengan satelit Bumi (orbit rendah dan juga orbit ruang angkasa lebih dalam), mungkin selama periode 1-3 orbit?
James C
@ James C Ya. Dalam jangka waktu yang lama, metode symplectic jauh lebih baik.
eccstartup
@eccstartup - Apa yang akan Anda pertimbangkan untuk waktu yang lama di sini? Karena bisa sepanjang satu orbit planet di sekitar Matahari, atau beberapa orbit satelit cuaca di sekitar Bumi dll.
James C
@ James C Saya belum melihat masalah besar itu. Tetapi untuk masalah model saya, dengan banyak orbit yang dihitung, metode symplectic memberikan orbit yang sangat sempurna.
eccstartup
Jadi, ini adalah saran untuk memprogram Anda memiliki versi metode Runge-Kutta implisit, yang mencakup banyak metode symplectic dengan urutan setinggi yang Anda inginkan.
eccstartup
4

merangkum beberapa poin:

  1. Jika itu adalah integrasi jangka panjang dari model non-dissapative, integrator symplectic adalah apa yang Anda cari.
  2. Kalau tidak, karena ini adalah persamaan gerak, metode Runge-Kutta Nystrom akan lebih efisien daripada transformasi ke sistem orde pertama. Ada metode RKN pesanan tinggi karena DP. Ada beberapa implementasi, seperti di sini di Julia mereka didokumentasikan dan inilah MATLAB .
  3. Metode Runge-Kutta orde tinggi hanya diperlukan jika Anda menginginkan solusi dengan akurasi tinggi. Jika toleransi lebih rendah maka RK pesanan ke-5 kemungkinan akan lebih cepat (untuk kesalahan yang sama). Hal terbaik untuk dilakukan jika Anda harus sering menyelesaikan ini adalah dengan menguji berbagai metode yang berbeda. Dalam serangkaian tolok ukur pada masalah 3-tubuh ini, kita melihat bahwa (untuk kesalahan yang sama) metode RK tingkat tinggi hanya merupakan peningkatan marjinal dalam kecepatan, meskipun sebagai kesalahan -> 0 Anda dapat melihat bahwa peningkatan sudah berjalan ke> 5x terhadap Dormand -Prince 45 ( DP5) ketika Anda melihat 4 digit akurasi (toleransi jauh lebih rendah untuk ini. Toleransi hanya rata-rata dalam masalah apa pun). Ketika Anda menarik toleransi bahkan semakin rendah peningkatan dari metode RK tingkat tinggi tumbuh, tetapi Anda mungkin perlu mulai menggunakan angka presisi yang lebih tinggi.
  4. Algoritma Dormand-Prince order 7/8 memiliki tableau urutan ke-8 yang berbeda dari metode DP853 dari Hairer's dop853dan DifferentialEquations.jl DP8(yang sama). Metode 853 yang terakhir tidak dapat diimplementasikan dalam versi tablo standar dari metode Runge-Kutta karena estimator kesalahannya adalah non-standar. Tetapi metode ini jauh lebih efisien dan saya tidak akan merekomendasikan bahkan menggunakan metode Fehlberg 7/8 atau DP 7/8 yang lebih lama.
  5. Untuk metode RK pesanan tinggi, metode Verner "Efficient" adalah standar utama. Itu muncul di tolok ukur yang saya tautkan. Anda bisa mengkodekannya ke dalam Boost diri Anda, atau menggunakan salah satu dari 2 paket yang mengimplementasikannya jika Anda menginginkannya lebih mudah (Mathematica atau DifferentialEquations.jl).
Chris Rackauckas
sumber
2

Saya ingin menambahkan bahwa sementara apa yang disarankan Geoff Oxberry untuk integrasi jangka panjang (menggunakan integrator symplectic) adalah benar, dalam beberapa kasus itu tidak akan berfungsi. Lebih khusus lagi, jika Anda memiliki kekuatan disipatif, sistem Anda tidak menyimpan energi lagi, dan karena itu Anda tidak dapat menggunakan integrator symplectic dalam kasus itu. Orang yang mengajukan pertanyaan itu berbicara tentang orbit Bumi yang rendah, dan orbit seperti itu menunjukkan sejumlah besar gaya hambat atmosfer, yang merupakan kekuatan disipatif yang menghalangi penggunaan integrator symplectic.

Dalam kasus tertentu (dan untuk kasus di mana Anda tidak dapat menggunakan / tidak memiliki akses ke / tidak ingin menggunakan integrator symplectic), saya akan merekomendasikan penggunaan integrator Bulirsch-Stoer jika Anda membutuhkan presisi dan efisiensi dalam jangka waktu yang lama. Ini bekerja dengan baik berdasarkan pengalaman, dan juga direkomendasikan oleh Numerical Recipes (Press et al., 2007).

viiv
sumber
Tidak, jangan rekomendasikan Numerical Recipes. Terutama, dalam banyak kasus Burlirsch-Stoer tidak direkomendasikan. Ini adalah masalah yang terkenal dengan buku itu. Lihat sekelompok bantahan dari para peneliti top di bidang ini di sini: uwyo.edu/buerkle/misc/wnotnr.html . Jika Anda menginginkan tolok ukur dalam hal ini, lihat buku pertama Hairer di mana Anda akan melihat BS hampir tidak pernah berhasil. Urutan yang lebih tinggi hanya lebih efisien ketika kesalahan cukup rendah, dan kami (dan yang lain) telah melakukan pembandingan untuk menunjukkan secara cukup konsisten bahwa itu hanya efisien untuk ketelitian sub-floating point.
Chris Rackauckas
Saya tidak dapat berbicara untuk NR terlalu banyak karena saya menggunakannya sebagian besar untuk ODE, tetapi bagi saya kelihatannya keluhan pada halaman yang Anda tautkan sudah tua dan telah dikuatkan oleh penulis NR dalam tanggapan mereka (akhir halaman), tapi ini di luar topik. Mengenai integrasi jangka panjang orbit dengan presisi tinggi (katakanlah, 13-14 digit) yang merupakan apa yang saya sebutkan dalam jawaban saya, terbukti sejak lama bahwa metode ekstrapolasi bekerja dengan baik (lihat bab Montenbruck & Gill tentang Numerical Integration). Makalah yang lebih baru menggunakannya juga, dan telah membuktikan kepada saya dan orang lain metode yang dapat diandalkan dan efisien.
viiv
M&G hanya mengujinya terhadap dop853 dan metode RK tingkat tinggi yang lebih modern, seperti yang disebabkan oleh Verner, jauh lebih efisien. M&G juga tampaknya hanya mengukur menggunakan evaluasi fungsi, yang merupakan indikator waktu yang lemah. Ini juga tidak waktu terhadap metode Runge-Kutta Nystrom yang khusus untuk ODE urutan ke-2 dan lebih efisien daripada metode RK orde pertama yang diterapkan ke urutan kedua dengan sedikit. Pada 13-14 digit BS mungkin kompetitif pada sebagian besar masalah, tetapi jauh dari pilihan yang jelas dan saya belum melihat diagram presisi-kerja dengan metode terbaru tidak setuju dengan itu.
Chris Rackauckas
M&G melakukan uji RKN terhadap RK, dan BS dan lainnya terhadap RKN (halaman 123-132 dan 151-154) dan mengatakan mereka adalah yang paling efisien dari metode RK (tidak termasuk Verner meskipun mereka mengutipnya). BS terbukti efisien pada 13-14 digit, yang merupakan klaim saya, saya telah melihatnya diuji terhadap dop853, ABM (12), Taylor, dan RK8 standar dan berkinerja baik. Saya harus mengakui bahwa saya belum melihatnya diuji terhadap RKN tetapi dari apa yang saya lihat dari M&G tidak jauh dari FILG11 misalnya. Saya benar-benar tertarik dengan Verner RK dan akan melihat tautan Anda di atas. Apakah Anda memiliki kertas yang menguji semuanya untuk dilihat?
viiv
Saya kembali dan menjalankan kembali tolok ukur di DiffEqBenchmarks.jl dan odextidak cenderung adil. Jadi setidaknya untuk ODE urutan pertama dan untuk toleransi >=1e-13, ekstrapolasi tampaknya tidak berfungsi dengan baik dan biasanya tidak mendekati. Ini sesuai dengan klaim di atas.
Chris Rackauckas