如何求解k=nπ对应的Sturm-Liouville方程并批量获取多阶解

你原有实现的fun和bc函数可以完全保留无需调整，只需要新增循环逻辑，为每阶特征值生成匹配的初始猜测即可，完整修改后代码如下：

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

# 原有右端项函数无需修改
def fun(x, y, p):  
    k = p[0]
    return np.vstack((y[1], -k**2 * y[0]))

# 原有边界条件函数无需修改
def bc(ya, yb, p):
    k = p[0]
    return np.array([ya[0], yb[0], ya[1] - k])

# 定义需要求解的最高阶数N
N = 10
# 存储所有求解结果
solutions = []
# 初始网格可共用，阶数较高时可适当加密网格提升收敛性
x = np.linspace(0, 1, 10)

for n in range(1, N+1):
    # 生成对应n阶解的初始y猜测，匹配解析解sin(nπx)的分布
    y = np.zeros((2, x.size))
    y[0] = np.sin(n * np.pi * x)
    # 可选：给导数赋近似初始值，加快收敛速度
    y[1] = n * np.pi * np.cos(n * np.pi * x)
    # k的初始猜测取和n成正比的数值即可，接近真实值nπ就能正常收敛
    k_guess = n * 3.0
    # 调用求解器
    sol = solve_bvp(fun, bc, x, y, p=[k_guess])
    # 验证收敛性后存储结果
    if sol.success:
        solutions.append(sol)
        print(f"n={n}，求解得到k={sol.p[0]:.4f}，预期值{n*np.pi:.4f}，误差={abs(sol.p[0]-n*np.pi):.2e}")
    else:
        print(f"n={n}求解失败：{sol.message}")

# 绘图展示所有阶数的解
x_plot = np.linspace(0, 1, 100)
plt.figure(figsize=(10,6))
for idx, sol in enumerate(solutions):
    n = idx + 1
    y_plot = sol.sol(x_plot)[0]
    plt.plot(x_plot, y_plot, label=f"n={n}, k={sol.p[0]:.3f}")
plt.xlabel("x")
plt.ylabel("y")
plt.legend(loc="upper right", bbox_to_anchor=(1.2, 1))
plt.title("Sturm-Liouville问题各阶特征解")
plt.grid()
plt.show()

关键修改说明：

• 通过循环遍历n=1到N，为每阶解单独生成初始猜测：第n阶特征解的解析形式是sin(nπx)，直接用这个表达式生成初始y值，比手动设置单点值的通用性更强，不管n多大都能保证收敛
• k的初始猜测直接取和n成正比的数值即可，和真实值nπ误差不大的情况下都能正常收敛
• 所有求解结果存储在solutions列表中，你可以后续按需取出任意阶的结果做进一步计算

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