求解具有求和项的二阶 ODE

Question

我有一个微分方程，我正在尝试用数值方法求解。

我正在尝试寻找当 n 很大（比如 50）时使用 scipy 中的循环或其他方法在 python 中解决这个问题的方法。

当n=1时我能够解决它。

这是我用来执行此操作的代码。

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 
t = np.linspace(0, 100, 1000) 
A = 2 
B = 3 
a = 2 
b = -3 
c = 2 
omega = 20 
def dSdt(t, S): 
    y, P = S 
    dSdt = [P, (A/a)*np.cos(omega*t) + (B/a)*np.sin(omega*t) - (b/a)*P - (c/a)*y] 
    return dSdt 
y_0 = 0 
P_0 = 0 
sol = solve_ivp(dSdt, 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)

但是我如何将其扩展到更一般的情况？

Answer 1

您的 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 dSdt(t, S): 
    y, P = S 
    dSdt = [P, (A/a)*np.cos(omega*t) + (B/a)*np.sin(omega*t) - (b/a)*P - (c/a)*y] 
    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(dSdt, 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)