Как вы кормите Свипи бвп только у вас есть БК?

Question

Единственный пример / документы, которые я могу найти, находятся на странице Scipy docs .

Чтобы проверить, я смотрю на не зависящий от времени уравнение Шреда в 1d бесконечной потенциальной яме. У этого есть точное аналитическое решение, найденное, решая DE и вставляя граничные условия 0 (0) = 0, ψ (L) = 0, и что функция soln к 1, но этот вопрос относится к решению любого DE, где BCs мы знаем, не для начального значения.

Вы можете решить эту проблему численно с помощью метода execute_ivp Сципи, начав с ψ (0) = 0, и выбрав соответствующий код для размещения ψ '(0), используя аналитическое решение. Можно использовать метод съемки, чтобы найти подходящее значение E, например, приведенное выше условие нормализации.

Это два набора BC: ψ (0) = 0 для обоих, нормализация для обоих и второе значение ψ для аналитического подхода и начальное значение ψ для подхода ivp. Кажется, решение Scipy solve_bvp предлагает решение, использующее первый набор чисел BC (поскольку мы обманываем, вставляя ψ '), но я не могу заставить его работать. Этот псевдокод описывает проблему и как я ожидаю, что API будет вести себя:

bcs = {0: (0, None), L: (0, None)} # Two BCs on ψ; no BCs on derivative
x_span = (0, L)

sol = solve_bvp(rhs, bcs, x_span)

На самом деле код выглядит примерно так, и я не могу заставить его работать:

def bc(ψ_a, ψ_b):
    return np.array([ψ_a[0], ψ_b[0]])

x_span = (0, L)
x_eval = np.linspace(x_span[0], x_span[1], int(1e5))

x_guess = np.array([0, L])
ψ_guess = np.array([[0, 1], [0, -1]])

res = solve_bvp(rhs_1d, bc, x_guess, ψ_guess)

Я не знаю, как построить функцию bc, и не знаю, почему догадки настроены так, как они есть. И не знаете, как я могу угадать значение ψ, не вставляя предположение для guess '. (Документы подразумевают, что вы можете) Также обратите внимание, что в документах показан пример, подразумевающий, что вы можете использовать solve_bvp также для нормализации BC, но не знаете, как подойти. (Пример слишком скудный)

Эквивалентный и рабочий код ivp, для ref: (Сравните с моим псевдокодом solve_bvp)

Код Python:

ψ_0 = (0, sqrt(2/L) * n*π/L)
x_span = (0, L)

sol = solve_ivp(rhs_1d, x_span, ψ_0)

Answer 1

Для задачи на собственные значения

-u''+V(x)u = c*u

с граничными условиями

u(0)=0=u(L)

и нормализация

int(u(x)^2, x=0 to L)=1 

установить интеграл в качестве третьего компонента. С собственным значением в качестве параметра это 4 измерения, учитывающие 4 граничных условия, дополнительные 2 состоят в том, что интеграл в 0 равен нулю и что интеграл в L имеет значение 1.

# some length
L = 10;

# some potential function
def V(x): return 1+(2*x-L)**2;

# the ODE function
def odesys(x,y,p):
    u,v,S = y; c=p[0]
    return [v, (V(x)-c)*u , u**2 ]

# the boundary conditions
def boundary(y0, yL, c):
    return [ y0[0], yL[0], y0[2], yL[2]-1 ]

Исходным предположением вы приблизительно выбираете, какую собственную функцию / собственное значение вы получите, более или менее.

n=11;
w = (np.pi*n)/L
x_init = np.linspace(0,L,4*n+1);
u_init = np.sin(w*x_init);
v_init = np.cos(w*x_init)*w;
y_init = [ u_init, v_init, x_init/L ]

Нет необходимости вкладывать слишком много точек в предположение, достаточно того, чтобы структура первого компонента была точно представлена.

Затем вызовите решатель с подготовленными данными, обратите внимание, что допуск по умолчанию составляет 1e-3, если вы хотите лучше, вы должны разрешить более точное подразделение. Если все работает нормально, нарисуйте решение.

res = solve_bvp(odesys, boundary, x_init, y_init, p=[w**2], max_nodes=10000, tol=1e-6)
print res.message

if res.success:
    x_disp = np.linspace(0,L,3001)
    y_disp = res.sol(x_disp)
    plt.plot(x_disp, y_disp[0])
    plt.title("eigenfunction to eigenvalue $\lambda=%.6f$"%res.p[0]);
    plt.grid(); plt.show()