Я относительно новичок в Python и пытаюсь использовать его для решения нелинейного дифференциального уравнения второго порядка, в частности уравнения Пуассона-Больцмана в электролите.Проблемы с ODEINT в python
phi''(r) + (2/r)*phi'(r) = (k^2)*sinh(phi(r))
По существу, он описывает распад электростатического потенциала (PHI) от заряженной поверхности в электролит со скоростью распада, управляемой с помощью параметра Я к.
- фи (г) - потенциал в точке г
- DPHI (г) - производная потенциала при г
- г - расстояние от поверхности (я решения при г = 1, R = R в этот случай)
и граничные условия
- фи (1) = 5
- DPHI (R) = 0
Проблема бита кода выглядит следующим образом
from scipy.integrate import odeint
from scipy.optimize import root
from pylab import * # for plotting commands
k = 0.5
phi = 5
dphi = -10
R = 21
def deriv(z,r): # return derivatives of the array z (where z = [phi, phi'])
return np.array([
(z[1]),
((k**2)*sinh(z[0]))-((2/r)*z[1])
])
result = odeint(deriv,np.array([phi,dphi]),np.linspace(1,R,1017), full_output = 1)
Обычно при низких значениях к интеграции работает отлично и я могу использовать корень из scipy.optimize для ее решения с учетом граничных условий, однако когда я пытаюсь использовать относительно большие значения к (0.5 и выше) работает интеграция с проблемами и выводит следующее сообщение об ошибке
Excess work done on this call (perhaps wrong Dfun type).
Пробежав его full_output = 1 и взглянули на tcur, кажется, что он плавно перечитывается до определенной точки, а затем колеблется дико от очень больших до очень малых значений. В той же точке nfe число оценок функций падает до нуля. При правильной работе параметр tcur работает плавно от 1 до R. Может ли кто-нибудь просветить меня, почему это происходит? Это проблема с моей реализацией или есть неустойчивость в уравнении?
Заранее спасибо за любую помощь,
Dave