Python使用数组参数调用bvp求解器求解四阶边值问题

核心思路 #

你预定义的一维t数组中每个元素对应独立的边值问题，直接遍历t数组对每个参数单独调用solve_bvp求解即可，建议用solve_bvp自带的fun_args参数传参替代全局变量，代码更易维护。

可运行示例代码 #

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

# 预定义的t数组，可根据你的需求修改
t_arr = np.linspace(1, 5, 10)

# 微分方程定义，增加t作为入参
def fun(x, y, t):
    return np.vstack([y[1], y[2], y[3], t/(t-y[0])])

# 边界条件不变
def bc(ya, yb) :
    return np.array([ya[0], ya[1], yb[0], yb[1]])

xmesh = np.linspace(0, 10, 100)
# 存储所有t对应的求解结果
sol_results = []

for t in t_arr:
    # 初始化猜测值
    y_init = np.zeros((4, len(xmesh)))
    # 传入当前t参数求解
    sol = solve_bvp(fun, bc, xmesh, y_init, fun_args=(t,))
    # 校验求解状态
    if sol.success:
        sol_results.append(sol)
    else:
        print(f"t={t:.2f}求解失败：{sol.message}")

# 批量绘制所有结果
plt.figure(figsize=(8, 5))
for idx, sol in enumerate(sol_results):
    plt.plot(xmesh, sol.sol(xmesh)[0], label=f't={t_arr[idx]:.2f}')
plt.xlabel('x')
plt.ylabel('w(x)')
plt.legend()
plt.show()

注意事项 #

• 如果你的t是随x变化的数组（即t是x的函数，不是全局常量），只需要在fun内部根据当前x值索引获取对应的t值即可，比如你预先将t和x的对应关系存为t_x = np.array([...])，则在fun中可以用current_t = t_x[np.argmin(np.abs(xmesh - x))]获取对应位置的t值。
• 注意t的取值范围不要过小，避免出现t - y[0]趋近于0的奇点，导致求解报错。

内容的提问来源于stack exchange，提问作者Anweshan