求解两个二阶微分方程

Question

我是初学者和MATLAB用户。我试图用Python用odeint解两个微分二阶方程组，但我的图与预期的不同。我认为我的函数中有些东西定义得很糟糕。 这是我的代码：

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

m1=100
m2=4.15
k1=120
k2=5
A1=1
A2=0
wn1=math.sqrt(k1/m1)
wn2=math.sqrt(k2/m2)
wf1=wn1
wf2=0
chi1=0.05
chi2=0.001
c1=2*chi1*wn1*m1
c2=2*chi2*wn2*m2

winicio=0.05
wpaso=0.01
wfin=2.5

T=2*math.pi/winicio
paso=0.05*T
tinicio=0
tfinal=T

yopos1=0
yovel1=0
yopos2=0
yovel2=0

def fy(y, t):
     return [y[2],
             y[3], 
             (-(c1+c2)*y[2] + c2*y[3] - (k1+k2)*y[0] + k2*y[1] + A1*math.sin(winicio*t))/m1,
             (-c2*(y[3]-y[2]) - k2*(y[1]-y[0]) + A2*math.sin(wf2*t))/m2]

tspan=np.linspace(tinicio, tfinal, 1000)
y0 = [yopos1, yopos2, yovel1, yovel2]
Ty = odeint(fy, y0, tspan)
ys1 = Ty[:,0]
ys2 = Ty[:,1]
vs1 = Ty[:,2]
vs2 = Ty[:,3]

plt.figure(1)          
plt.subplot(2,1,1)
plt.xlabel("time")
plt.ylabel("position")
plt.title("GDL 1:Position-Time")
plt.plot(tspan,ys1,'r');
plt.legend()
plt.subplot(2,1,2)
plt.xlabel("time")
plt.ylabel("position")
plt.title("GDL 1:Position-Time")
plt.plot(tspan,ys2,'r');

plt.figure(2)          
plt.subplot(2,1,1)
plt.xlabel("Time")
plt.ylabel("Velocity")
plt.title("GDL 1:Velocity-Time")
plt.plot(tspan,vs1,'r',);
plt.legend()
plt.subplot(2,1,2)
plt.xlabel("Time")
plt.ylabel("Velocity")
plt.title("GDL 2:Velocity-Time")
plt.plot(tspan,vs2,'r');

但是我在Matlab中使用了相同的方法和相同的函数，我的图形工作得非常好：

fy=@(t,y) [y(3);
          y(4);
          (-(c1+c2)*y(3)+c2*y(4)-(k1+k2)*y(1)+k2*y(2)+A1*sin(winicio*t))/m1;
          (-c2*(y(4)-y(3))-k2*(y(2)-y(1))+A2*sin(wf2*t))/m2];

tspan=tinicio:paso:tfin; % Intervalo de tiempo
yopos=[yopos1;yopos2;yovel1;yovel2]; % Condiciones iniciales
[t,y]=ode45(fy,tspan,yopos); % Resolución ec. dif
n=length(t);

figure(1) % Gráficas 1.Posición-t y Velocidad-t

subplot(2,1,1)
plot(t,y(:,1),'b'); % GDL 1: posición-tiempo
xlim([0 t(n)]);
hold on; grid on; grid minor
xlabel('Tiempo')
ylabel('Posición')
title('Gdl 1: Posición-Tiempo')

subplot(2,1,2)
plot(t,y(:,2),'r'); % GDL 2: posición-tiempo
xlim([0 t(n)]);
hold on; grid on,grid minor
xlabel('Tiempo')
ylabel('Posición')
title('Gdl 2: Posición-Tiempo')

有人能帮我说一下是怎么回事吗