你如何只给Scipy的bvp提供你拥有的BCs？

Question

我能找到的唯一示例/文档都在Scipy docs page上。

为了测试，我看的是一维无限势阱中与时间无关的Schrod方程。这有一个简洁的解析解，通过求解DE，并插入边界条件ψ(0) = 0，ψ(L) = 0，函数soln为1，但这个问题适用于求解任何我们知道的BCs不是初值的DE。

您可以使用Scipy的solve_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函数，也不知道为什么猜测是这样设置的。并且不确定如何在不插入对ψ‘的猜测的情况下猜测ψ的值。(文档暗示您可以)同样值得注意的是，文档显示了一个示例，暗示您也可以将solve_bvp用于规范化BC，但不确定如何实现。(示例过于稀疏)

ref:的等效有效的ivp代码(与我的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处的积分为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()

​

​