如何在Python中求解中间点固定值的ODE系统？

你的问题本质是边值问题（BVP），而solve_ivp是初值问题（IVP）求解器，仅支持设置区间一端的初始值。在Python中可以通过两种常见方式实现中间点固定的ODE求解：

方法一：拆分区间+反向积分 #

把原区间[-1, 1]拆分为[-1, 0]和[0, 1]两段：

• 对于[0, 1]：直接以t=0处的(E₀, phi₀)为初始值，用solve_ivp正向积分到t=1。
• 对于[-1, 0]：将时间变量反向（令τ = -t），把区间转化为[0, 1]，同样以t=0处的(E₀, phi₀)为初始值，积分后再把时间变量转换回原t轴。

具体代码实现：

import numpy as np
from scipy.integrate import solve_ivp

def ode(t, Ephi, g, l):
    return [
        2 * (Ephi[0]**2) * g * np.sin(2 * Ephi[1]),
        Ephi[0] * (l + 2 * g * np.cos(2 * Ephi[1]))
    ]

# 反向积分用的ODE（τ = -t，dτ/dt = -1，所以dy/dτ = -dy/dt）
def ode_rev(τ, Ephi, g, l):
    dydt = ode(-τ, Ephi, g, l)
    return [-dydt[0], -dydt[1]]

l_value = 1
g_value = 1
t_min = -1
t_max = 1
num_points = 1000

E_0 = 1
phi_0 = 0

# 求解右半区间 [0, 1]
sol_right = solve_ivp(
    ode,
    (0, t_max),
    (E_0, phi_0),
    t_eval=np.linspace(0, t_max, num_points//2),
    args=(g_value, l_value)
)

# 求解左半区间 [-1, 0]，通过反向积分τ∈[0,1]对应t∈[-1,0]
sol_left = solve_ivp(
    ode_rev,
    (0, -t_min),  # τ从0到1，对应t从0到-1
    (E_0, phi_0),
    t_eval=np.linspace(0, -t_min, num_points//2),
    args=(g_value, l_value)
)
# 转换回原t轴
sol_left.t = -sol_left.t
# 反转数组顺序，让t从-1到0递增
sol_left.t = sol_left.t[::-1]
sol_left.y = sol_left.y[:, ::-1]

# 合并两段解
t_full = np.concatenate([sol_left.t[:-1], sol_right.t])
E_full = np.concatenate([sol_left.y[0, :-1], sol_right.y[0]])
phi_full = np.concatenate([sol_left.y[1, :-1], sol_right.y[1]])

# 验证t=0处的值
print(f"t=0处E的值: {E_full[np.argmin(np.abs(t_full))]}")
print(f"t=0处phi的值: {phi_full[np.argmin(np.abs(t_full))]}")

方法二：使用边值问题求解器solve_bvp #

scipy.integrate.solve_bvp专门用于求解边值问题，支持自定义边界条件。我们可以将问题转化为：在区间[-1, 1]上，约束t=0处的E和phi等于给定值。

具体代码实现：

import numpy as np
from scipy.integrate import solve_bvp

def ode(t, y, g, l):
    E, phi = y
    dE_dt = 2 * E**2 * g * np.sin(2 * phi)
    dphi_dt = E * (l + 2 * g * np.cos(2 * phi))
    return np.vstack([dE_dt, dphi_dt])

# 边界条件函数：残差为0时满足约束
def bc(ya, yb, y_mid, E0, phi0):
    # ya是t=-1处的状态，yb是t=1处的状态，y_mid是t=0处的状态
    return np.array([
        y_mid[0] - E0,  # t=0处E=E0
        y_mid[1] - phi0   # t=0处phi=phi0
    ])

l_value = 1
g_value = 1
t_min = -1
t_max = 1
num_points = 1000

E_0 = 1
phi_0 = 0

# 构造包含t=0的求解网格
t_grid = np.linspace(t_min, t_max, num_points)
# 找到t=0对应的索引
mid_idx = np.argmin(np.abs(t_grid))

# 初始猜测解（可以用常数猜测）
y_guess = np.zeros((2, num_points))
y_guess[0, :] = E_0
y_guess[1, :] = phi_0

# 调用solve_bvp，传递中间点的状态到边界条件函数
sol = solve_bvp(
    lambda t, y: ode(t, y, l_value, g_value),
    lambda ya, yb: bc(ya, yb, y_guess[:, mid_idx], E_0, phi_0),
    t_grid,
    y_guess
)

# 检查求解是否成功
if sol.success:
    print("求解成功")
    # 验证t=0处的值
    t0_idx = np.argmin(np.abs(sol.x))
    print(f"t=0处E的值: {sol.y[0, t0_idx]}")
    print(f"t=0处phi的值: {sol.y[1, t0_idx]}")
else:
    print(f"求解失败: {sol.message}")

注意：solve_bvp对初始猜测的质量比较敏感，如果初始猜测不合适可能导致求解失败，可以根据问题的物理意义调整初始猜测。

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