Я новичок и ПОЛЬЗОВАТЕЛЬ 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')
Пожалуйста, кто-нибудь может мне помочь и сказать, что это неправильно?
![https://i.stack.imgur.com/zXnHv.png](https://i.stack.imgur.com/zXnHv.png)