Решение двух дифференциальных уравнений второго порядка - PullRequest
0 голосов
/ 07 марта 2020

Я новичок и ПОЛЬЗОВАТЕЛЬ 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

...