你如何只喂你拥有的 BCs 的 Scipy 的 bvp？

Posted 2023-02-18

【中文标题】你如何只喂你拥有的 BCs 的 Scipy 的 bvp？

【英文标题】：How do you feed Scipy's bvp only the BCs you have?

【发布时间】：2019-04-02 21:05:24

【问题描述】：

我能找到的唯一示例/文档位于 Scipy docs page。

为了测试，我正在查看一维无限势阱中与时间无关的 Schrod eq。这通过求解 DE 找到了一个简洁的解析解，并插入 ψ(0) = 0，ψ(L) = 0 的边界条件，并且 func soln 为 1，但这个问题适用于求解任何 BC 的 DE我们知道不是为了初始值。

您可以使用 Scipy 的solve_ivp 以数值方式求解它，方法是从ψ(0) = 0 开始，并使用解析soln 作弊以适当地放置ψ'(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，但不确定如何处理。 （例子太稀疏了）

等效且有效的 ivp 代码，用于参考：（与我的 solve_bvp 伪代码比较）

Python 代码：

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

sol = solve_ivp(rhs_1d, x_span, ψ_0)

【问题讨论】：

"...我无法让它工作"。出了什么问题？你有错误吗？如果是这样，显示 complete 错误消息。或者你是否得到了你知道不正确的输出？

我怀疑你的 x_guess 应该有更多的点——对于某个整数 N，类似于 x_guess = np.linspace(0, L, N)。然后 ψ_guess 必须是一个形状为 (2, N) 的数组，其中包含一个粗略的对边值问题的解的猜测。查看solve_bvp 文档字符串中的“Bratu”示例；那里的猜测非常粗略！

我没有收到错误，但得到我知道是错误的输出。这可能是因为我不确定如何将 BC（例如，我在两点求解的 fn 的已知值，在我的第一个代码 sn-p 中描述）扩展为 solve_bvp api 使用的复杂参数函数。

【参考方案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()