python - 如何在solve_bvp 中设置适当的边界条件？

Question

我正在尝试解决以下边界值问题:
D*(dS^2/dz^2) - v*(dS/dz) - S/(S+k_s) = 0
其中 D、v 和 k 是常量(参见下面的代码)。边界条件为:
S(z = 0) = 23.8 和 dS/dz(z = 无穷大) = 0
我正在尝试使用 scipy.integrate.solve_bvp在 (0,100) 的 z 范围内解决问题，但我很难设置边界值。我在网上看到的所有例子都没有真正描述用于计算边界值的函数中发生了什么，作为 solve_bvp 的输入。 .我还有一个额外的复杂性，即在无穷远处有一个边界值。这是我到目前为止:

# constants
D = 200.6
v = 0.02
k = 0.20

# split second order ODE into 2 first order ODEs
def sulfate_profile(z,U):
    dUdz = [U[1], (v*U[1] + U[0]/(U[0]+k))/D]
    return dUdz

# I initiated these boundary conditions but I don't really know what they mean...
def bc(za, zb):
    return np.array([za[0], zb[1]])

z_list = np.linspace(0,100,1000)
z_a = np.zeros((2, z_list.size))

# solve SO4 concentration ODE over depth range (z_list)
sulfate_reduction = solve_bvp(sulfate_profile,bc,z_list,z_a)

这给了我一个输出，但我确定这是错误的。我真的很感激任何建议 - 提前致谢!

Answer 1

好吧，我在发布此答案之前犹豫了，因为我宁愿更彻底地检查解决方案，但正如我所见，这仍然是我发现的所有内容。请参阅下面的代码，但首先要进行一些解释:

首先，我更改了您的 ODE，使其与您问题中的 ODE 一致

我添加了一个修改后的右侧(有趣)添加参数 - 这是为了覆盖现在从转换 ODE 中改变的边界条件。我一直接近 scipy 教程，因此解释相匹配。您可能想为参数 p 尝试不同的起始值。

边界条件(bc1)相同 - 命名 bc 参数的逻辑是 Sn:n=0 边界 a，n=1 边界 b； [0] = y, [1] = y' - 所以名字覆盖了边界的边，索引为导数

因为我不知道问题的物理原理，所以我添加了一个图表来显示解决方案。这看起来似乎正确，但请注意无穷条件:现在在 b 边界(“另一侧”)处设置为 0。无穷大无法覆盖，我看到的唯一解决方案是向外扩展该限制。

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

# constants
D = 200.6
v = 0.02
k = 0.20

# split second order ODE into 2 first order ODEs
'''def sulfate_profile(U, z):
    print(U, z)
    dUdz = [U[1], (v*U[1] + U[0]/(U[0]+k))/D]
    print(dUdz)
    return dUdz'''

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html
def sulfate_profile(z, S): # compare to doc: x->z, y->S
    return np.vstack((S[1], v/D*S[1] + S[0]/(S[0]+k)/D))

def fun(z, S, p):
    k0, k1 = p[0], p[1]
    return np.vstack((S[1], k0*S[1] + k1*S[0]/(S[0]+k)))

# I initiated these boundary conditions but I don't really know what they mean...
def bc(S0, S1):
    return np.array([S0[0], S1[0]])

def bc1(S0, S1, p):
    k0 = p[0]
    return np.array([S0[0], S1[0], S0[1] - k0*23.8/(23.8+k), S1[0]]) # Sn: n=0 boundary a, n=1 boundary b; [0] = y, [1] = y',...

z_list = np.linspace(0, 100, 100)
S_a = np.zeros((2, z_list.size))
S_a[0, 0] = 23.8
S_a[0, 99] = 0 # b-boundary which is inf - expand here

# solve SO4 concentration ODE over depth range (z_list)
#sulfate_reduction = solve_bvp(sulfate_profile, bc, z_list, S_a)
sulfate_reduction = solve_bvp(fun, bc1, z_list, S_a, p=[1,1])

#print(sulfate_reduction)
x_plot = np.linspace(0, 1, 100)
y_plot_a = sulfate_reduction.sol(x_plot)[0]
plt.plot(x_plot, y_plot_a, label='sulfate_b_dpth')
plt.show()

这导致

这看起来像一个深度分布(如果是这样的话)。