你怎么只给你的BCs喂Scipy的bvp？

Question

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

为了测试，我在一维无限势阱中观察与时间无关的薛定谔方程。这有一个简洁的解析解，通过求解DE，并插入ψ（0）=0，ψ（L）=0的边界条件，并且函数解为1，但是这个问题适用于求解我们知道的bc不是初值的任何DE。在

可以用Scipy的solve_ivp从ψ（0）=0开始进行数值求解，并用解析解适当地放置ψ'（0）。可以使用shooting方法找到合适的E值，例如上面的规范化条件。在

这是两组BCs：ψ（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)

实际上，代码看起来像这样，但我无法让它工作：

^{pr2}$

我不知道如何构建bc函数，也不知道为什么猜测是这样设置的。我不知道我怎么能猜出ψ的值，而不需要插入ψ'的猜测。（文档暗示您可以）另外，文档显示了一个示例，暗示您也可以将solve_bvp用于规范化BC，但不确定如何处理。（示例太少）

等效的工作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

有边界条件

^{pr2}$

和规范化

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()