Решение сложной проблемы BVP с помощью scipy solve_bvp - PullRequest
0 голосов
/ 20 сентября 2019

У меня сложная задача БВП (почти уравнение Шредингера).Мой код работает, но решение полностью неверно или равно нулю.Что я делаю неправильно?Я также получил правильное решение от 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), это не поможет.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...