您可以更新模型以直接包含整个总和，为什么不将其参数化：

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

def model(t, x, a, b, c, As, Bs, omegas):
    rhs = np.sum([
        A * np.cos(omega * t) + B * np.sin(omega * t)
        for A, B, omega in zip(As, Bs, omegas)
    ])
    return np.array([
        x[1],
        (rhs - b * x[1] - c * x[0]) / a 
    ])

字符串
然后执行以下调用，求解IVP：

a = 2
b = -3
c = 2
As = [2, 8, 14] 
Bs = [3, 11, -2]
omegas = [20, 40, 60]
t = np.linspace(0, 100, 1000) 

solution = integrate.solve_ivp(
    model, [t.min(), t.max()], [0, 0],
    t_eval=t, args=(a, b, c, As, Bs, omegas)
)

型
由于线性，它返回与@YakovDan相同的解决方案，但只调用了solve_ivp。

你的ODE是线性的，所以你可以只对参数求和，然后求和：

import numpy as np 
import matplotlib.pyplot as plt 
import scipy as sp 
from scipy.integrate import odeint 
from scipy.integrate import solve_ivp 
#SOLVING THE GOVERNING DIFFERENTIAL EQUATION 

def get_dSdt(i):
    def dSdt(t, S): 
        y, P = S 
        dSdt = [P, (A[i]/a)*np.cos(omega[i]*t) + (B[i]/a)*np.sin(omega[i]*t) - (b/a)*P - (c/a)*y]
        return dSdt 
    return dSdt     

n = 3 # for example

t = np.linspace(0, 100, 1000) 
A = [2, 8, 14] 
B = [3, 11, -2] 
a = 2 
b = -3 
c = 2 
omega = [20, 40, 60]
y_0 = 0 
P_0 = 0 
disp_value = 0
for i in range(n):
    sol = solve_ivp(get_dSdt(i), t_span = (0, max(t)), y0 = [y_0, P_0], t_eval=t) 
    time = sol.t 
    disp_value += sol.y[0]  
plt.plot(time, disp_value)

字符串
注意，函数get_dSdt接受一个整数i并返回dSdt函数，除了为索引i定义的参数值