У меня сложная задача БВП (почти уравнение Шредингера).Мой код работает, но решение полностью неверно или равно нулю.Что я делаю неправильно?Я также получил правильное решение от Maple, проблем не было.Также я не понимаю, почему нет проблемы с графиком, когда это должна быть комплексная функция.
import numpy as np
from scipy.integrate import odeint
from scipy.integrate import solve_bvp as bvp
import matplotlib.pyplot as plt
mp = 938.2720813
mn = 939.5654133
mu = (mn + mp)/4
h2m = hbar**2/(2*mu)
V0 = 20
Rv = 1.5
Q0 = 1.5
Rq = 4.5
EIm = 0.3
ERe = 1
V = lambda x : -V0*np.exp(-x/Rv)
Q = lambda x : -Q0*np.exp(-x/Rq)
def fun(x, y):
return np.vstack((y[1], -( Q(x)/ h2m ) - ((ERe + 1j * EIm) *y[0]/ h2m ) + V(x)*y[0]/h2m - (2/y[0])* y[1]))
def bc(ya, yb):
return np.array([ya[0], yb[0]])
x = np.linspace(0, 1000, 10000)
y_a = np.zeros((2, x.size), dtype=np.complex)
# print(x.size)
i = 0
while i < x.size - 1:
i = i + 1
y_a[0, i] = 1000* 1j
y_a[1, i] = 1j
from scipy.integrate import solve_bvp
res_a = solve_bvp(fun, bc, x, y_a)
x_plot = np.linspace(0, 1000, 10000)
y_plot_a = res_a.sol(x_plot)[0]
import matplotlib.pyplot as plt
plt.plot(x_plot, y_plot_a, label='y_a')
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.show()
Upd: Я исправил ошибку в уравнении.Результат все еще не верен.Но есть и другая ошибка - деление на ноль.Как этого избежать?Например, если я выберу x = np.linspace (0.1, 1000, 10000), это не поможет.