Saya perlu mengevaluasi secara integral angka-angka di bawah ini:
di mana , dan . Di sini adalah fungsi Bessel yang dimodifikasi dari jenis kedua. Dalam kasus khusus saya, saya memiliki , dan .
Saya menggunakan MATLAB, dan saya telah mencoba fungsi bawaan integral
dan quadgk
, yang memberi saya banyak kesalahan (lihat di bawah). Saya secara alami telah mencoba banyak hal lain juga, seperti mengintegrasikan dengan bagian-bagian, dan menjumlahkan integral dari ke .
Jadi, apakah Anda punya saran mengenai metode mana yang harus saya coba selanjutnya?
UPDATE (menambahkan pertanyaan)
Saya membaca makalah @Pedro yang ditautkan, dan saya pikir itu tidak terlalu sulit untuk dipahami. Namun, saya punya beberapa pertanyaan:
- Apakah boleh menggunakan sebagai basis-elemen , dalam metode Levin univariat yang dijelaskan?
- Bisakah saya menggunakan metode Filon saja, karena frekuensi osilasi diperbaiki?
Kode contoh
>> integral(@(r) sin(x*r).*sqrt(E(r)),0,Inf)
Warning: Reached the limit on the maximum number of intervals in use. Approximate
bound on error is 1.6e+07. The integral may not exist, or it may be difficult to
approximate numerically to the requested accuracy.
> In funfun\private\integralCalc>iterateScalarValued at 372
In funfun\private\integralCalc>vadapt at 133
In funfun\private\integralCalc at 84
In integral at 89
ans =
3.3197e+06
sumber
Jawaban:
Saya telah menulis integrator saya sendiri
quadcc
, yang mengatasi secara substansial lebih baik daripada integrator Matlab dengan singularitas, dan memberikan perkiraan kesalahan yang lebih andal.Untuk menggunakannya untuk masalah Anda, saya melakukan hal berikut:
Fungsi
f
ini sekarang adalah integrasi Anda. Perhatikan bahwa saya baru saja menetapkan nilai lama apa pun kex
.Untuk mengintegrasikan pada domain yang tidak terbatas, saya menerapkan substitusi variabel:
yaitu mengintegrasikan∞
g
dari 0 ke 1 harus sama dengan mengintegrasikanf
dari 0 hingga . Transformasi yang berbeda dapat menghasilkan hasil kualitas yang berbeda: Secara matematis semua transformasi harus memberikan hasil yang sama, tetapi transformasi yang berbeda dapat menghasilkan hasil yang lebih halus, atau lebih mudah diintegrasikan .g
Saya kemudian memanggil integrator saya sendiri
quadcc
,, yang dapat menanganiNaN
s di kedua ujungnya:Perhatikan bahwa perkiraan kesalahan sangat besar, yaitu
quadcc
tidak memiliki banyak kepercayaan pada hasil. Melihat fungsi, meskipun, ini tidak mengherankan karena berosilasi pada nilai tiga orde besarnya di atas integral yang sebenarnya. Sekali lagi, menggunakan transformasi interval yang berbeda dapat menghasilkan hasil yang lebih baik.Anda mungkin juga ingin melihat metode yang lebih spesifik seperti ini . Ini sedikit lebih terlibat, tetapi jelas metode yang tepat untuk jenis masalah ini.
sumber
integral
(1e-10 saya pikir), tetapi 1.7e + 07 masih benar-benar besar. Mungkin transformasi lain akan baik, seperti yang Anda sebutkan.Seperti yang ditunjukkan Pedro, metode tipe-Levin adalah metode terbaik untuk masalah-masalah seperti ini.
Apakah Anda memiliki akses ke Mathematica? Untuk masalah ini, Mathematica akan mendeteksi dan menggunakannya secara default:
Berikut adalah plot pada rentang nilai x:
Anda juga dapat secara manual menentukan metode tipe-Levin tertentu untuk diterapkan, yang dalam hal ini dapat menghasilkan sedikit peningkatan kinerja:
Lihat dokumentasi untuk rincian metode tipe-Levin di Mathematica .
sumber
Jika Anda tidak memiliki akses ke Mathematica, Anda bisa menulis metode tipe-Levin (atau osilasi khusus lainnya) di Matlab seperti yang disarankan Pedro.
Apakah Anda menggunakan perpustakaan chebfun untuk Matlab? Saya baru belajar itu berisi implementasi metode tipe-Levin dasar di sini . Implementasinya ditulis oleh Olver (salah satu ahli di bidang quadrature osilasi). Itu tidak berurusan dengan singularitas, subdivisi adaptif dll, tetapi mungkin hanya apa yang Anda butuhkan untuk memulai.
sumber
Transformasi yang direkomendasikan oleh Pedro adalah ide bagus. Sudahkah Anda mencoba untuk bermain-main dengan parameter dalam fungsi "quadgk" Matlab? Misalnya, menggunakan transformasi Pedro, Anda dapat melakukan hal berikut:
quadgk(f, 0.0+eps, 1.0-eps, 'AbsTol', eps, 'MaxIntervalCount', 100000)
Menggunakan ini memberi saya solusi:
-2184689.50220729
dan hanya membutuhkan 0,8 detik (menggunakan nilai yang disebutkan di atas: x = 10)
Walter Gander dan Walter Gautschi memiliki makalah tentang quadrature adaptif dengan Matlab kode yang dapat Anda gunakan juga (tautan di sini )
sumber