решить_ивп от -х до + х

Question

Я пытался использовать solve_ivp или solve_bvp для решения проблемы, с которой я столкнулся, но я не добиваюсь никакого прогресса.Я думаю, что код, который я здесь, будет работать, но я не могу получить правильный диапазон.по какой-то причине я не могу понять, что диапазон всегда изменяется от 0 до x, а не от -x до x, может кто-нибудь помочь мне исправить эту часть для solve_ivp?

вот код, уменьшенный до минимума

from pylab import *
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
import numpy as np
import itertools

a=1
B=4
L= B+a
Vmax= 50
Vpot = False

N = 1000                  # number of points to take
psi = np.zeros([N,2])     # Wave function values and its derivative (psi and psi')
psi0 = array([0,1])   # Wave function initial states
Vo = 50
E = 0.0                   # global variable Energy  needed for Sch.Eq, changed in function "Wave function"
b = L                     # point outside of well where we need to check if the function diverges
x = linspace(-B-a, L, N)    # x-axis
def V(x):
    '''
    #Potential function in the finite square well.
    '''
    if -a <=x <=a:
        val = Vo
    elif x<=-a-B:
        val = Vmax
    elif x>=L:
        val = Vmax
    else:
        val = 0
    return val

def SE(z, p):
    state0 = p[1]
    state1 = 1.0*(V(z) - E)*p[0]
    return array([state0, state1])

def Wave_function(energy):
    global E
    E = energy
    #        odeint(func, y0, t)
    #     solve_ivp(fun, t_span, y0)
    psi = solve_ivp(SE, [-B-a, L], psi0, max_step = ((B+a+L)/(N)))
    return psi.y

def main():
    # main program        
    f2 = figure(2)
    plot(x, Wave_function(9.8)[0][:1000])
    grid()
    f2.show()

if __name__ == "__main__":
    main()

это то, что код дает мне в конце.с правой стороны все в порядке, но с левой стороны это не так.Я зависим от работы обеих сторон, а не для визуальных.

edit: для благотворительности это то, как должна выглядеть потенциальная функция:

окончательный график должен выглядеть аналогичнона это:

Answer 1

Ваш код в целом выглядит нормально.Однако, учитывая вашу цифру для потенциальной энергии, значение для Vo должно быть * Vo = 10 .Кроме того, в вашей функции main вы только изображаете волновую функцию как решение уравнения Шредингера.Ниже, вот что я предлагаю вам в качестве возможного решения вашей проблемы, при условии, что я правильно понял вашу проблему:

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import numpy as np

a=1
B=4
L= B+a
Vmax= 50


N = 1000                  # number of points to take
psi = np.zeros([N,2])     # Wave function values and its derivative (psi and psi')
psi0 = np.array([0,1])   # Wave function initial states
Vo = 10     # Not 50, in order to conform your figure of the potential energy
E = 0.0         # global variable Energy  needed for Sch.Eq, changed in    
                # function "Wave function"
b = L           # point outside of well where we need to check if the function diverges
x = np.linspace(-L, L, N) # linspace(-B-a, L, N)    # x-axis


def V(x):
    '''
    Potential function in the finite square well.
    '''
    if -a <=x <=a:
        val = Vo
    elif x<= -L:  # -a-B:
        val = Vmax
    elif x>=L:
        val = Vmax
    else:
        val = 0
    return val

def SE(z, p):
    state0 = p[1]
    state1 = 1.0*(V(z) - E)*p[0]
    return array([state0, state1])

def Wave_function(energy):
    global E
    E = energy
    psi = solve_ivp(SE, [-B-a, L], psi0, max_step = ((B+a+L)/(N)))
    return psi.y

def main():
    # main program       
    plt.figure()
    plt.subplot(121)
    plt.plot(x, Wave_function(9.8)[0][:1000])
    plt.grid()
    plt.title("Wave function")
    plt.xlabel(r"$ x $")
    plt.ylabel(r"$\psi(x)$")
    plt.subplot(122)

    potential = np.vectorize(V) # Make the function 'V(x)' to also work on array

    pot = potential(x) # Potential ernergy in the defined domain of 'x'
    t = [-L, -a, a, L] # the singular value of x
    y = potential(t)   # the potential energy at thos singular value of 'x'
    # But to conform your figure I'll just do y = 0 * y
    plt.plot(x, pot, t, 0*y, 'ko')
    plt.title("Potential Energy")
    plt.xlabel(r"$ x $")
    plt.ylabel(r"$V(x)$")
    plt.show()

if __name__ == "__main__":
    main()

Выходные данные следующие:

Answer 2

Не уверен, поможет ли это, но может дать подсказку.Дело не в том, что solve_ivp не работает для -x до 0, но ваша функция V может быть неправильной.Я заметил, что волна начинает появляться после того, как V уменьшается с Vmax до 0.

Этот код:

%matplotlib inline
from pylab import *
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
import numpy as np
import itertools

a=1.
B=4.
L= B+a
Vmax= 50.

N = 10000
E = 9.8

def V(x):
    if -L <= x <= -B:
        return Vmax
    else:
        return 0

def SE(z, p):
    state0 = p[1]
    state1 = (V(z) - E)*p[0]
    return array([state0, state1])

def Wave_function():
    return solve_ivp(SE, [-L, L], [0., 1.], max_step=2*L/N)

result = Wave_function()
plot(result.t, result.y[0], color='tab:blue')

дает вам "ожидаемый" результат: