Evaluasi numerik integral sangat berosilasi

11

Dalam kursus lanjutan ini pada aplikasi teori fungsi kompleks pada satu titik dalam latihan integral yang sangat berosilasi

I(λ)=cos(λcosx)sinxxdx

harus didekati untuk nilai λ menggunakan metode sadel di bidang kompleks.

Karena sifatnya yang sangat berosilasi, integral ini sangat sulit untuk dievaluasi menggunakan sebagian besar metode lain. Ini adalah dua fragmen grafik integrand untuk λ=10 pada skala yang berbeda:

cos (10 cos (x)) sin (x) / x

Perkiraan asimtotik urutan terkemuka adalah

I1(λ)=cos(λ14π)2πλ

dan penyempurnaan lebih lanjut (jauh lebih kecil) menambahkan istilah

I2(λ)=18sin(λ14π)2πλ3

Grafik nilai yang diperkirakan sebagai fungsi λ terlihat sebagai berikut:

I (lambda) kira-kira

Sekarang tiba pertanyaan saya: untuk melihat secara visual seberapa bagus perkiraannya, saya ingin membandingkannya dengan "nilai sebenarnya" dari integral, atau lebih tepatnya ke perkiraan yang baik untuk integral yang sama menggunakan algoritma independen. Karena kecilnya koreksi subleading, saya berharap ini akan menjadi sangat dekat.

λtanh(sinh)

sekitar mpmath

Akhirnya saya mencoba keberuntungan saya dengan integrator Monte-Carlo menggunakan sampel penting yang saya terapkan, tetapi saya tidak berhasil mendapatkan hasil yang stabil juga.

λ>1

doetoe
sumber
Apakah fungsinya seimbang?
nicoguaro
Ya, itu bahkan
doetoe
Sudahkah Anda mencoba mengubah integral Anda menjadi ODE?
nicoguaro
1
x
1
λλxx(cos(λxcosx)sincx)

Jawaban:

12

Gunakan teorema Plancherel untuk mengevaluasi integral ini.

f,g

I=f(x)g(x)dx=F(k)G(k)dk

F,Gf,gsinx/xrect(k)cos(λcosx)λ|Jn(x)|n>|x|

πJ0(λ)[0,2π]

smh
sumber
Terima kasih, ini ide yang sangat bagus!
doetoe
7

πN+π2

Asimptotik

I(λ)2πλ[cos(λπ4)+c1sin(λπ4)λ+c2cos(λπ4)λ2+c3sin(λπ4)λ3+]
c1=18

int := NIntegrate[Cos[l*Cos[x]]*Sinc[x], {x, 0, 20.5*Pi}]; 
Plot[{l*(Sqrt[2*l/Pi]*int - Cos[l-Pi/4]), Sin[l-Pi/4]/8}, {l, Pi/4, 20}]

Sebagai output Anda mendapatkan sinus yang cukup bagus yang bertepatan dengan yang Anda dapatkan di atas.

18

Jika Anda ingin menemukan koefisien berikut, sedikit lebih canggih jika diperlukan. Gagasan kode di bawah ini adalah untuk mengambil beberapa nilai batas atas yang tinggi dan untuk "rata-rata" hasilnya.

J[l_?NumericQ] := Block[{n=500},
  f[k_] := NIntegrate[Cos[l*Cos[x]]*Sinc[x], {x,0,(n+k)*Pi+Pi/2},
  Method->{"DoubleExponential"}, AccuracyGoal->14, MaxRecursion->100];
  1/2*((f[0]+f[1])/2+(f[1]+f[2])/2)
]
t = Table[{l, l^2*(Sqrt[2*l/Pi]*J[l] - Cos[l-Pi/4] - 1/8*Sin[l-Pi/4]/l)}, 
    {l, 4*Pi+Pi/4, 12*Pi+Pi/4, Pi/36}];
Fit[t, Table[Cos[l-Pi/4+Pi/2*n]/l^n, {n, 0, 10}], l]

c2=9128,c3=751024,c4=367532768,

Penjelasan

Contoh sederhana

S(x)=0xsin(y)ydy.
S()=π2

sinus

S(x)

SN=n=1N(1)nn.
SSN+12(1)N+1N+1.
S(x)0πN+π2sinxxdx
max|S(x)|

Masalahmu

Kembali ke integral dari program Konstantin dan Yaroslav, Anda dapat melihat bahwa ia berperilaku persis sama dengan sinus - integral sebagai fungsi batas atas. Itu berarti Anda hanya perlu menghitung nilai dengan . Di bawah ini adalah plot dari beberapa nilai seperti itu dengan .

Ix0(λ)=20x0cos(λcos(x))sinc(x)dx
x0=πN+π2λ=12π

tab = Table[{x0, 2*NIntegrate[Cos[12*Pi*Cos[x]]*Sinc[x], {x, 0, x0}, 
     Method->{"DoubleExponential"}, AccuracyGoal->12, MaxRecursion->100]},
    {x0, 10*Pi+Pi/2, 30*Pi+Pi/2, Pi}];
tab1 = Table[(tab[[i]] + tab[[i+1]])/2, {i,1,Length[tab]-1}];
ListPlot[{tab, tab1}]

mnrt

Di sini Anda dapat melihat hasil dari metode percepatan lain. Saya mengatur ulang jumlah parsial dengan cara berikut dan mendapatkan urutan baru yang konvergennya jauh lebih cepat. Trik itu juga berguna jika Anda ingin mengevaluasi integral dengan presisi tinggi.

SN=12(SN+SN+1)
SN

David Saykin
sumber
Bagus! Apakah instruktur tentu saja profesor kehidupan nyata Anda? Kursus mereka fantastis, meskipun sangat tangguh dan serba cepat
doetoe
@ terlalu ya, saya mahasiswa Konstantin. Dia membagikan kepada saya tautan ke pertanyaan Anda.
David Saykin
6

Metode Ooura untuk integral sinus Fourier bekerja di sini, lihat:

Ooura, Takuya, dan Masatake Mori, Formula eksponensial ganda yang kuat untuk integral tipe Fourier. Jurnal matematika komputasi dan terapan 112.1-2 (1999): 229-241.

Saya menulis sebuah implementasi dari algoritma ini tetapi tidak pernah bekerja untuk membuatnya cepat (oleh, katakanlah caching node / bobot), tetapi tetap saja, saya mendapatkan hasil yang konsisten pada segala sesuatu di luar presisi float:

float     = 0.0154244
double    = 0.943661538060268
long dbl  = 0.943661538066058702
quad      = 0.943661538066060288748574485677942
oct       = 0.943661538066060288748574485680878906503533004997613278231689169604876
asymptotic= 0.944029734

Berikut kodenya:

#include <iostream>
#include <boost/math/quadrature/ooura_fourier_integrals.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>

template<class Real>
Real asymptotic(Real lambda) {
    using std::sin;
    using std::cos;
    using boost::math::constants::pi;
    Real I1 = cos(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/lambda);
    Real I2 = sin(lambda - pi<Real>()/4)*sqrt(2*pi<Real>()/(lambda*lambda*lambda))/8;
    return I1 + I2;
}

template<class Real>
Real osc(Real lambda) {
    using std::cos;
    using boost::math::quadrature::ooura_fourier_sin;
    auto f = [&lambda](Real x)->Real { return cos(lambda*cos(x))/x; };
    Real omega = 1;
    Real Is = ooura_fourier_sin<decltype(f), Real>(f, omega);
    return 2*Is;
}

template<class Real>
void print(Real val) {
   std::cout << std::defaultfloat;
   std::cout << std::setprecision(std::numeric_limits<Real>::digits10);
   std::cout <<  val <<  " = " << std::hexfloat << val;
}

int main() {
    using boost::multiprecision::float128;
    float128  s = 7;
    std::cout << "Asymptotic = " << std::setprecision(std::numeric_limits<float128>::digits10) << asymptotic(s) << "\n";
    std::cout << "float precision = ";
    print(osc(7.0f));
    std::cout << "\n";
    std::cout << "double precision= ";
    print(osc(7.0));
    std::cout << "\n";
    std::cout << "long double     = ";
    print(osc(7.0L));
    std::cout << "\n";
    print(osc(s));

    print(osc(boost::multiprecision::cpp_bin_float_oct(7)));
    std::cout << "\n";
}

Anda tidak dapat benar-benar melihat perbedaan antara quadrature dan asimptotik karena mereka berada tepat di atas satu sama lain kecuali sebagai : λ0masukkan deskripsi gambar di sini

pengguna14717
sumber
Terima kasih, ini sangat bagus! Saya belum membuatnya berfungsi, instalasi boost saya tidak kompatibel, tapi saya mengunduh versi terbaru sekarang.
doetoe
Hanya untuk memastikan: di 23 Anda memiliki cos (lambda * cos (x)) / x, tanpa faktor sin (x) dari integrand. Apakah itu ooura_fourier_sin yang menganggap faktor ini berdosa (x) untuk melipatgandakan integrand yang diteruskan ke sana?
doetoe
Saya membuatnya bekerja. Tampaknya semua dependensinya dan hanya header, jadi saya bahkan tidak perlu menginstal atau kompilasi (kecuali untuk executable). Saya harap ini akan dimasukkan dalam boost!
doetoe
@doetoe: Ya, faktor implisit dan dimasukkan ke dalam bobot quadrature. Saya sedang mempertimbangkan untuk mempercepat dengan refactoring untuk men-cache node dan bobot; kita lihat saja nanti! sin(x)
user14717
@doetoe: Ini digabung menjadi Boost 1.71. API sedikit berbeda dari jawaban ini.
user14717