Поскольку решением сферически-симметричного уравнения Пуассона-Больцмана (ПБ) является электрический потенциал (psi), который спадает до 0 на больших расстояниях (r -> infty), я установил граничное условие так, что psi=0 на моем конечное значение расстояния дляsolve_bvp. Мое начальное условие состоит в том, что psidot распадается кулоновски при r = r0, где r0 — радиус моего центрального иона. Однако затухающая функция представляет собой гиперболический грех, который затухает до 0 при малых значениях фунтов на квадратный дюйм, но никогда до 0 на самом деле. Я считаю, что это является источником проблем с гиперболическим грехом у обоих решателей. Есть идеи, как это преодолеть?
Нажмите здесь, чтобы просмотреть уравнение ПБ
import numpy as np
from scipy.integrate import solve_bvp, odeint
import matplotlib.pyplot as plt
def model(r,x,a,b):
y = x[0]
dy = x[1]
xdot = [[],[]]
xdot[0] = dy
xdot[1] = -2/r*dy + a*sinh(b*y)
return xdot
# boundary conditions for solve_bvp:
def bc(at0, at1, a, b):
# Values at r=initial:
yat0, ydotat0 = at0
# Values at r=final:
yat1, ydotat1 = at1
# These return values are what we want to be 0:
return [ydotat0 - ydot0, yat1]
# initial condition for ODEINT
y0 = z*e/(4*pi*eps*r0)
ydot0 = -z*e/(4*pi*eps*r0**2)
r = np.linspace(r0,100*r0,201)
C1 = 0.1 #molar concentration
a = 2*e*C1*Nav*1000/(eps) #2eC/eps
b = e/(k_b*T) #e/kT
ystart = odeint(model,[y0,ydot0],r, args=(a,b,), tfirst=True)
result = solve_bvp(lambda x,r: model(x,r,a=a,b=b), lambda at0, at1: bc(at0,at1,a=a,b=b), r, ystart.T)
выход:
\anaconda3\lib\site-packages\ipykernel_launcher.py:6: RuntimeWarning: overflow encountered in sinh
\anaconda3\lib\site-packages\ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in add
\anaconda3\lib\site-packages\scipy\integrate\_bvp.py:40: RuntimeWarning: invalid value encountered in subtract
hi = y_new[i] - y[i]
Область, допустимая для вычисления с плавающей запятой экспоненты и, следовательно, гиперболического синуса, довольно мала, с верхним пределом ниже 1000. Если решатель выдает пробные решения за пределами этого диапазона, вы получаете ошибку. Замените sinh на функцию, которая для аргументов точна до 200, а затем продолжает линейно или с какой-либо другой отсечкой.