Saya mencoba menghitung urutan lebih tinggi (misalnya m=0
,, n=46
) momen Zernike untuk beberapa gambar. Namun, saya mengalami masalah mengenai polinomial radial (lihat wikipedia ). Ini adalah polinomial yang didefinisikan pada interval [0 1]. Lihat kode MATLAB di bawah ini
function R = radial_polynomial(m,n,RHO)
R = 0;
for k = 0:((n-m)/2)
R = R + (-1).^k.*factorial(n-k) ...
./ ( factorial(k).*factorial((n+m)./2-k) .* factorial((n-m)./2-k) ) ...
.*RHO.^(n-2.*k);
end
end
Namun, ini jelas mengalami masalah numerik dekat RHO > 0.9
.
Saya mencoba refactoring untuk polyval
berpikir mungkin ada beberapa algoritma di belakang layar yang lebih baik tapi itu tidak menyelesaikan apa pun. Mengubahnya ke perhitungan simbolis memang membuat grafik yang diinginkan tetapi sangat lambat bahkan untuk grafik sederhana seperti yang ditunjukkan.
Apakah ada cara yang stabil secara numerik untuk mengevaluasi polinomial tingkat tinggi seperti itu?
matlab
polynomials
numerical-limitations
Sanchises
sumber
sumber
Jawaban:
Dalam tulisan ini , Honarvar dan Paramesran memperoleh metode yang menarik untuk menghitung polinomial radial Zernike dengan cara rekursif yang sangat bagus. Formula rekursi secara langsung sangat mudah, tanpa pembagian atau perkalian dengan bilangan bulat besar: Saya akan merekomendasikan untuk melihat gambar 1 di Honarvar dan Paramesran kertas, yang dengan jelas menggambarkan ketergantungan antara berbagai polinomial Zernike.Rmn(ρ)=ρ(R|m−1|n−1(ρ)+Rm+1n−1(ρ))−Rmn−2(ρ)
Ini diimplementasikan dalam skrip Oktaf berikut:
Misalnya, gambar yang dihasilkan oleh kode ini menunjukkan bahwa dengan , dan , pembatalan bencana terjadi di dekat , jika polinom radial Zernike dikomputasi melalui polinomial Jacobi. Oleh karena itu, kita juga harus khawatir tentang keakuratan polinomial Zernike tingkat rendah.m=22 n=112 ρ=0.7
Metode rekursif tampaknya jauh lebih cocok untuk menghitung polinomial Zernike tingkat tinggi ini dengan cara yang stabil. Namun demikian, untuk dan , perbedaan maksimum antara Jacobi dan metode rekursif adalah (hanya?) , Yang mungkin cukup akurat untuk aplikasi Anda.m=0 n=46
1.4e-10
sumber
jacobiPD
, tidak seperti pembatalan bencana umum.JacobiPD
dari jawabannya . Ini berfungsi dengan baik untuk polinomial tingkat rendah. Misalnya, dengan , arbitrary , dan arbitrary , perbedaan antara kedua metode lebih kecil dari . Meskipun istilah individual dalam penjumlahan adalah kecil, mereka mungkin menjadi besar setelah dikalikan dengan . Selain itu mereka memiliki tanda bolak-balik, yang merupakan resep sempurna untuk pembatalan bencana.6.9e-13
JacobiPD
factorial(n+a) * factorial(n+b)
1/(factorial(s)*factorial(n+a-s)*factorial(b+s)*factorial(n-s)) * ((x-1)/2).^(n-s).*((x+1)/2).^s * factorial(n+a) * factorial(n+b)
1.4e18
-2.1
Solusi yang mungkin (disarankan oleh @ gramatester) adalah dengan menggunakan polinomial Jacobi. Ini menghindari masalah pembatalan katastropik dalam menambahkan koefisien polinom besar dengan evaluasi polinomial 'naif'.
Polinomial Zernike radial dapat diekspresikan oleh polinomial Jacobi sebagai berikut (lihat persamaan (6) )
Namun dalam MATLAB, penggunaan
jacobiP(n,a,b,x)
tidak dapat diterima lambat untuk vektor / matriks besarx=rho
. ThejacobiP
Fungsi sebenarnya adalah bagian dari Toolbox simbolik, dan evaluasi dari polinomial ditangguhkan untuk mesin simbolik, yang perdagangan kecepatan untuk presisi sewenang-wenang. Implementasi manual polinomial Jacobi dengan demikian diperlukan.Karena parameter ke fungsi Jacobi semuanya nonnegatif ( , , ), kita dapat menggunakan ekspresi berikut (lihat Wikipedia , perhatikan bahwa saya mengisi nilai untuk )α=m β=0 n∗=(n−m/2) s P(α,β)n(ρ)=(n+α)!(n+β)!⋅∑s=0n[1s!(n+α−s)!(β+s)!(n−s)!(x−12)n−s(x+12)s]
Dalam MATLAB, ini diterjemahkan menjadi (Jacobi
p olice d epartmentP olynomial, ' D ouble' implementasi)Polinomial Zernike radial yang sebenarnya dengan demikian (untuk
m=abs(m)
)Catatan: jawaban-diri ini hanya solusi praktis; jangan ragu memberi tag pada jawaban lain yang menjelaskan mengapa ini berhasil.
sumber