Baru-baru ini saya telah membandingkan berbagai pemecah non-linear dari scipy dan sangat terkesan dengan contoh Newton-Krylov dalam Scipy Cookbook di mana mereka memecahkan persamaan persamaan diferensial orde dua dengan istilah reaksi non-linear dalam sekitar 20 baris kode.
Saya memodifikasi kode contoh untuk menyelesaikan persamaan Poisson non-linear ( juga disebut persamaan Poisson-Boltzmann , lihat halaman 17 dalam catatan ini) untuk heterostruktur semikonduktor, yang memiliki bentuk,
(Ini adalah fungsi residual yang diteruskan ke pemecah.)
Ini adalah masalah elektrostatik di mana dan p ( x , ϕ ) adalah fungsi nonlinier untuk bentuk n i ( x ) e - ( E i ( x , ϕ ) - E f ) . Detail di sini tidak penting, tetapi intinya adalah bahwa fungsi non-linear bervariasi secara eksponensial dengan ϕ sehingga fungsi residual dapat bervariasi pada rentang yang sangat besar ( 10 - 6 - 10 16 )dengan sedikit perubahan pada .
Saya secara numerik memecahkan persamaan ini dengan Scipy's Newton-Krylov, tetapi tidak akan pernah bertemu (pada kenyataannya ia akan selalu melaporkan kesalahan dengan menghitung Jacobian). Saya beralih dari pemecah Newton-Krylov ke fsolve (yang didasarkan pada MINPACK hybrd) dan berhasil pertama kali!
Apakah ada alasan umum mengapa Newton-Krylov tidak cocok untuk masalah tertentu? Apakah persamaan input perlu dikondisikan entah bagaimana?
Mungkin diperlukan lebih banyak informasi untuk berkomentar, tetapi mengapa menurut Anda fsolve bekerja dalam kasus ini?
sumber
sol = newton_krylov(func, guess, method='gmres')
) memperbaiki masalah. Tidak terlalu yakin mengapa, tetapi siapa pun dengan masalah ini mungkin mempertimbangkan untuk melakukan hal yang sama.Jawaban:
Ada dua masalah yang cenderung Anda temui.
Pengondisian buruk
Pertama, masalahnya tidak dikondisikan, tetapi jika Anda hanya memberikan residu, Newton-Krylov membuang separuh digit signifikan Anda dengan membedakan residual hingga residual untuk mendapatkan aksi Jacobian:
Perhatikan bahwa masalah yang sama berlaku untuk metode kuasi-Newton, meskipun tanpa perbedaan yang terbatas. Semua metode yang dapat diukur untuk masalah dengan operator yang tidak kompak (mis., Persamaan diferensial) harus menggunakan informasi Jacobian untuk prakondisi.
Kemungkinan besar
fsolve
tidak menggunakan informasi Jacobian atau menggunakan metode dogleg atau pergeseran untuk membuat kemajuan dengan metode "gradient descent", meskipun Jacobian dasarnya tunggal (yaitu, pembedaan terbatas akan memiliki banyak "noise" dari aritmatika presisi terbatas). Ini tidak dapat diskalakan danfsolve
kemungkinan akan semakin lambat saat Anda meningkatkan ukuran masalahnya.Globalisasi
Jika masalah linier dipecahkan secara akurat, kita dapat mengesampingkan masalah yang berkaitan dengan masalah linier (Krylov) dan fokus pada masalah yang disebabkan oleh nonlinier. Minima dan nonsmooth lokal memiliki konvergensi yang lambat atau menyebabkan stagnasi. Poisson-Boltzmann adalah model yang halus jadi jika Anda mulai dengan tebakan awal yang cukup baik, Newton akan bertemu secara kuadratik. Sebagian besar strategi globalisasi melibatkan semacam kelanjutan untuk menghasilkan tebakan awal berkualitas tinggi untuk iterasi akhir. Contohnya termasuk kelanjutan kotak (misalnya, Multigrid Penuh), kelanjutan parameter, dan kelanjutan pseudotransient. Yang terakhir ini umumnya berlaku untuk masalah steady-state, dan menawarkan beberapa teori konvergensi global, lihat Coffey, Kelley, dan Keyes (2003) . Pencarian muncul makalah ini, yang mungkin berguna bagi Anda:Shestakov, Milovich, dan Noy (2002) Solusi dari persamaan Poisson-Boltzmann nonlinier menggunakan kelanjutan pseudo-transient dan metode elemen hingga . Kelanjutan pseudotransient terkait erat dengan algoritma Levenberg-Marquardt.
Bacaan lebih lanjut
sumber