具有第二个成员列表的非均匀二阶微分方程

Question

你好。我正在从事学校项目，而我仍然是Python的初学者。我需要解决这个微分方程，但我不知道如何，因为第二个成员是列表。 这是我的代码

import numpy as np
import scipy.integrate as integr
t=np.linspace(0,799,96)
Q=np.array([60.50,113.93,199.59,292.04,376.85,459.97,526.13,570.23,604.16,623.66,626.21,627.06,628.75,622.81,600.76,585.50,582.10,567.69,536.30,501.53,489.66,480.33,451.49,411.63,372.61,332.75,304.76,276.77,242.85,211.47,178.39,107.99,24.03,-36.18,-13.28,19.79,39.3,51.17,41.84,28.27,18.94,19.79,27.42,38.45,40.99,31.66,21.48,12.16,4.52,-2.25,-10.73,-17.52,-20.06,-17.52,-14.98,-14.98,-16.67,-25.15,-34.48,-36.18,-31.94,-27.02,-24.30,-20.06,-16.67,-14.13,-15.82,-20.06,-26.85,-33.63,-37.03,-35.33,-31.09,-32.79,-37.03,-39.57,-40.42,-37.03,-31.09,-25.15,-26.85,-32.79,-40.42,-44.66,-46.36,-47.20,-48.05,-51.45,-53.99,-51.45,-44.66,-36.18,-29.39,-24.30,-20.06,-17.52]
)
Pmes=np.array([70.10,77.17,85.26,87.62,90.33,92.02,93.41,94.89,96.07,96.79,97.56,98.63,99.45,99.96,99.71,99.45,99.5,99.09,98.58,97.91,97.25,96.58,96.02,95.56,95.15,94.43,94.07,94.33,94.69,94.53,93.61,91.97,90.49,85.67,81.88,83.11,89.87,89.77,89.31,89.05,89.1,89.31,89.62,89.31,88.9,88.64,88.33,87.82,87.41,87.26,87.21,86.85,86.23,85.82,85.41,84.95,84.13,83.47,83.21,83.11,82.85,82.44,81.88,81.47,81.01,80.6,80.34,80.09,79.93,79.68,79.37,78.91,78.45,78.14,77.78,77.73,77.52,77.32,77.12,76.76,76.24,75.83,75.63,75.37,74.91,74.14,73.48,73.02,72.71,72.25,71.74,71.38,71.12,70.86,70.61,70.81])
#Conditions initiales
P0=70
dP0=70
#valeur de la résistance périphérique
Rp=0.63
#valeur de Rp4,Rc4,C4 et L pour le modèle de Windkessel à 4 éléments
Rc4=0.045
C4=2.53
L=0.0054
n=96
c=800/n
h=10**(-2)
A=np.array([[0,1],[-Rc4/(C4*L*Rp),-Rc4/L -1/(C4*Rp)]])
a=np.matrix(A)
#dérivée de Q
def D(F,i):
    if i!=0 and i!=n-1:
        return (F[i+1]-F[i-1])/2*h
    elif i==0:
        return (F[i+1]-F[i])/h
    else:
        return (F[i]-F[i-1])/h
#dérivé seconde de Q        
def DD(F,i):
    if i!=0 and i!=n-1:
        return (F[i+1]+F[i-1]-2*F[i])/(2*(h*h))
    elif i==0:
        return (F[i+2]-2*F[i+1]+F[i])/(h*h)
    else:
        return (F[i]-2*F[i-1]+F[i-2])/(h*h)
#Second membre de l'équation différentielle        
QQ=[DD(Q,i)*Rc4+D(Q,i)*(1/C4 +Rc4/(C4*Rp))+Q[i]*Rc4/(C4*Rp) for i in range(n)]
Pp=np.empty(n)
Pp[0]=[P0,dP0]
    for i in range(n-1):
        Pp[i+1]=c*a*np.matrix(Pp[i])+QQ[i]
P4WK=[Pp[i][0] for i in range(n)]
plt.plot(t,P4WK,'b',label='4WK')
plt.plot(t,Pmes,'g',label='Mesurée')
plt.title('La pression de l\'aorte ascedante pour le modèles de Windkessel à trois éléments')
plt.xlabel('temps(ms)')
plt.ylabel('Pression aortique(mmHg)')
plt.legend()
plt.grid()
plt.show()```