RuntimeWarning: переполнение в sinh при использовании ODEINT иsolve_bvp для решения уравнения Пуассона – Больцмана.

avatar
Ilan Shumilin
9 августа 2021 в 05:57
80
0
1

Поскольку решением сферически-симметричного уравнения Пуассона-Больцмана (ПБ) является электрический потенциал (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]
Источник
Lutz Lehmann
9 августа 2021 в 07:47
1

Область, допустимая для вычисления с плавающей запятой экспоненты и, следовательно, гиперболического синуса, довольно мала, с верхним пределом ниже 1000. Если решатель выдает пробные решения за пределами этого диапазона, вы получаете ошибку. Замените sinh на функцию, которая для аргументов точна до 200, а затем продолжает линейно или с какой-либо другой отсечкой.

Ответы (0)