如何用Python扩展求解含求和项的二阶常微分方程（n取大值时）

你已经完成了n=1时的二阶常微分方程求解，要扩展到n=50这类大值场景，核心是将原方程中的单简谐激励替换为n项简谐项的求和，利用numpy向量化运算提升效率，同时保持scipy.integrate.solve_ivp的调用逻辑不变。

核心思路 #

原n=1时的激励项是(A/a)*cos(ωt) + (B/a)*sin(ωt)，扩展到n项后，需要计算所有n组(A_i/a)*cos(ω_i t) + (B_i/a)*sin(ω_i t)的和，再代入二阶ODE的右端项。使用numpy向量化运算替代Python循环，能显著提升n取大值时的计算效率。

完整实现代码 #

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

# 求和项数量
n = 50
# 时间区间与采样点
t = np.linspace(0, 100, 1000) 
# ODE参数
a = 2 
b = -3 
c = 2 

# 生成n组激励参数（可根据实际需求替换为自定义规则）
A_array = np.random.uniform(1, 3, n)  # 示例：1-3之间的随机A值
B_array = np.random.uniform(2, 4, n)  # 示例：2-4之间的随机B值
omega_array = np.random.uniform(10, 30, n)  # 示例：10-30之间的随机角频率

def dSdt(t, S): 
    y, P = S 
    # 计算n项激励的求和（向量化运算，高效快速）
    sum_term = np.sum(
        (A_array / a) * np.cos(omega_array * t) + 
        (B_array / a) * np.sin(omega_array * t)
    )
    # 二阶ODE的状态导数
    dydt = P
    dPdt = sum_term - (b/a)*P - (c/a)*y
    return [dydt, dPdt] 

# 初始条件
y_0 = 0 
P_0 = 0 

# 求解微分方程
sol = solve_ivp(
    dSdt, 
    t_span=(0, max(t)), 
    y0=[y_0, P_0], 
    t_eval=t
) 

# 可视化结果
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='位移 $y(t)$')
plt.xlabel('时间 $t$')
plt.ylabel('位移')
plt.title(f'$n={n}$ 项求和激励下的二阶ODE解')
plt.legend()
plt.grid(True)
plt.show()

关键说明 #

1. 参数自定义：如果你的求和项参数（A_i、B_i、ω_i）有特定规律（比如按等差数列、幂次规律生成），直接替换A_array、B_array、omega_array的生成代码即可，无需使用随机数。
2. 效率优化：numpy向量化运算会将求和操作转为底层C实现，比Python原生for循环快数十倍，适合n=50甚至更大的场景。
3. 状态维度不变：不管n取多大，原方程始终是二阶常微分方程，状态变量仅为[y, P]（位移和速度），因此solve_ivp的调用逻辑无需修改。

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