Возникли проблемы с ODEINT в python

Я относительно новичок в Python и пытаюсь использовать его для решения нелинейного дифференциального уравнения второго порядка, в частности уравнения Пуассона-Больцмана в электролите.

phi''(r) + (2/r)*phi'(r) = (k^2)*sinh(phi(r))

По сути, он описывает спад электростатического потенциала (phi) вдали от заряженной поверхности в электролите со скоростью спада, определяемой параметром k.

  • phi(r) - потенциал на r
  • dphi(r) - производная потенциала по r
  • r - расстояние от поверхности (в данном случае я решаю от r = 1 до r = R)

и граничные условия

  • фи (1) = 5
  • дфи (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)

Как правило, для низких значений k интеграция работает нормально, и я могу использовать root из scipy.optimize для ее решения с учетом граничных условий, однако, когда я пытаюсь использовать относительно большие значения k (0,5 и выше), интеграция сталкивается с проблемами и выдает следующую ошибку

Excess work done on this call (perhaps wrong Dfun type).

Запустив его с full_output = 1 и взглянув на параметр tcur, кажется, что он плавно подсчитывает до определенного момента, а затем резко колеблется от очень больших до очень маленьких значений. В той же точке nfe количество вычислений функции падает до нуля. При правильной работе параметр tcur плавно изменяется от 1 до R. Может ли кто-нибудь просветить меня, почему это происходит? Это проблема с моей реализацией или в уравнении есть нестабильность?

Заранее благодарю за любую помощь,

Дэйв


scipy poisson differential-equations odeint nonlinear-functions
person David Gillespie    schedule 12.12.2014    source источник

Ответы (1)


arrow_upward
0
arrow_downward

ODE, вероятно, нестабилен. Связанное уравнение

phi''(r) = (k^2)*sinh(phi(r))

имеет решение для k=1 (для соответствующих начальных условий при r=1):

phi(r) = 2 arctanh(sin(r))

Решение имеет особенность в точке r=pi/2. Численный решатель ОДУ не сможет интегрировать уравнение за пределами этой точки. Вполне вероятно, что аналогичное уравнение с членом первой производной (которым в любом случае должно быть пренебрежимо мало вблизи сингулярностей) ведет себя аналогично.

Фактическая проблема, которая у вас есть, заключается в том, что метод стрельбы с использованием решателя ОДУ не является хорошим способом решения проблем с граничными значениями - вы должны использовать методы словосочетания, которые достаточно стабильны. См., например. http://www.scholarpedia.org/article/Boundary_value_problem и ссылки в нем.

Для программного обеспечения Python см. https://pypi.python.org/pypi?%3Aaction=search&term=boundary+value+problem&submit=search

Обычно также очень легко написать такие решатели самостоятельно, так как единственный необходимый шаг — это дискретизация задачи до набора (многих) уравнений, после чего root может их решить.

person pv.    schedule 14.12.2014