Metode untuk integrasi numerik dari osilasi integral yang sulit

25

Saya perlu mengevaluasi secara integral angka-angka di bawah ini:

0sinc(xr)rE(r)dr

di mana , dan . Di sini adalah fungsi Bessel yang dimodifikasi dari jenis kedua. Dalam kasus khusus saya, saya memiliki , dan .E(r)=r4(λκ2+r2)ν5/2Kν5/2(λκ2+r2)xR+λ,κ,ν>0Kλ=0.00313κ=0.00825ν=0.33

Saya menggunakan MATLAB, dan saya telah mencoba fungsi bawaan integraldan 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 .kxπ(k+1)xπ

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?xkψk
  • 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

torbonde
sumber
Apa dalam integral Anda? x
Pedro
Setiap bilangan real positif. Saya baru saja memperbarui posting saya.
torbonde
Jika Anda dapat menunjukkan beberapa kode dan kesalahan, mungkin tidak terlalu sulit untuk menyelesaikannya. Tentu saja silakan coba membaca kesalahan dengan hati-hati terlebih dahulu dan lihat apakah Anda dapat membuatnya sendiri.
Dennis Jaheruddin
Saya akan membuat komentar hari ini dengan beberapa kode dan kesalahan. Atau besok.
torbonde
Oke, jadi saya lupa. Tapi sekarang saya memperbarui posting saya dengan contoh (saya telah membagi integral menjadi dua dengan menghitung secara eksplisit). sinc
torbonde

Jawaban:

12

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:

>> lambda = 0.00313; kappa = 0.00825; nu = 0.33;
>> x = 10;
>> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* besselk(-nu-5/2,lambda*sqrt(kappa^2 + r.^2));
>> sincp = @(x) cos(x)./x - sin(x)./x.^2;
>> f = @(r) sincp(x*r) .* r .* sqrt( E(r) );

Fungsi fini sekarang adalah integrasi Anda. Perhatikan bahwa saya baru saja menetapkan nilai lama apa pun ke x.

Untuk mengintegrasikan pada domain yang tidak terbatas, saya menerapkan substitusi variabel:

>> g = @(x) f ( tan ( pi / 2 * x ) ) .* ( 1 + tan ( pi * x / 2 ).^2 ) * pi / 2;

yaitu mengintegrasikan gdari 0 ke 1 harus sama dengan mengintegrasikan fdari 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 menangani NaNs di kedua ujungnya:

>> [ int , err , npoints ] = quadcc( g , 0 , 1 , 1e-6 )
int =
  -1.9552e+06
err =
   1.6933e+07
npoints =
       20761

Perhatikan bahwa perkiraan kesalahan sangat besar, yaitu quadcctidak 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.

Pedro
sumber
Terima kasih banyak. Saya akan melihat metode yang berbeda. Untuk tujuan saya, kesalahan tidak harus sekecil standar dalam persamaan integral(1e-10 saya pikir), tetapi 1.7e + 07 masih benar-benar besar. Mungkin transformasi lain akan baik, seperti yang Anda sebutkan.
torbonde
@ cimrg.joe: Perhatikan bahwa perkiraan kesalahan adalah perkiraan berdasarkan kesalahan absolut, antara lain, pada nilai absolut maksimum dari integrand. Dalam beberapa kasus ekstrim, nilai yang dikembalikan mungkin sebenarnya cukup baik. Jika Anda mencari akurasi sepuluh digit, maka saya sangat merekomendasikan menggunakan metode tipe-Levin yang saya sebutkan di akhir posting saya.
Pedro
Saya mungkin tidak membutuhkan akurasi sepuluh digit, tetapi saya pikir saya perlu setidaknya lima digit. Bisakah metode Anda menghasilkan itu?
torbonde
Metode ini tidak dapat menjamin ketepatan semacam itu untuk integral Anda karena nilai-nilai di ujung kanan interval adalah beberapa urutan besarnya lebih besar dari integral itu sendiri.
Pedro
11

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:

In[1]:= e[r_] := 
 r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]]

In[2]:= {l, k, v} = {0.00313, 0.00825, 0.33};

In[3]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3]]

Out[3]= -112494.

Berikut adalah plot pada rentang nilai x:

In[4]:= ListLinePlot[
 Table[NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
   PrecisionGoal -> 3], {x, .5, 10, 0.1}]]

Plot dari x = 0,5 hingga x = 10

Anda juga dapat secara manual menentukan metode tipe-Levin tertentu untuk diterapkan, yang dalam hal ini dapat menghasilkan sedikit peningkatan kinerja:

In[5]:= method = {"LevinRule", "Kernel" -> {Cos[r x], Sin[r x]}, 
   "DifferentialMatrix" -> {{0, -x}, {x, 0}}, 
   "Amplitude" -> {(
     3497.878840962873` Sqrt[(
      r^4 BesselK[-2.17`, 
        0.00313` Sqrt[
         0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
        r^2)^1.415`])/
     x, -((3497.878840962873` Sqrt[(
       r^4 BesselK[-2.17`, 
         0.00313` Sqrt[
          0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
         r^2)^1.415`])/(r x^2))}, "AdditiveTerm" -> 0};

In[6]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3, Method -> method]]

Out[6]= -112495.

Lihat dokumentasi untuk rincian metode tipe-Levin di Mathematica .

Andrew Moylan
sumber
Sayangnya saya tidak memiliki akses ke Mathematica - hanya MATLAB. Saya hanya akan memperbarui pertanyaan saya dengan beberapa pertanyaan tambahan, tentang makalah @Pedro yang tertaut.
torbonde
OK, seperti yang Anda katakan Anda harus puas dengan Matlab. Saya akan menambahkan jawaban lain tentang itu.
Andrew Moylan
5

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.

Andrew Moylan
sumber
Saya sudah berpikir tentang menerapkan metode Levin sendiri, tetapi saya belum yakin apakah saya siap untuk tantangan ini. Saya pikir saya perlu memahami metode ini sedikit lebih baik. Mungkin saya bisa berbicara dengan penasihat saya tentang hal itu. Lagi pula, alasan saya bertanya tentang metode Filon, adalah bahwa mereka tampak lebih mudah diimplementasikan. Dan karena saya tidak membutuhkan akurasi yang sangat tinggi, tetapi ini adalah bagian dari tesis master saya, kesulitan menimbang.
torbonde
Saya telah melihat perpustakaan chebfun (yang mengesankan) dan contoh integrasi Levin. Tapi saya tidak bisa menjalankannya. Saya sebenarnya telah memposting pertanyaan mengenai hal itu di sini .
torbonde
0

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 )

Joel Vroom
sumber