弹簧/阻尼器计算和绘图

Question

给定两个带有阻尼器/弹簧的系统：

第一个系统的 simulink 模型，步长为 2，最终值为 0.5：

具有相同输入的第二个系统的Simulink：

我必须找到使用 dsolve 和 ode45 的代码才能使用 Simulink 生成相同的图形。给定值是：

m1 = 500
c1 = 1200
k1 = 25000
k2 = 15000
m2 = 50

我试图找到dsolve，但无法解决。所以我不得不使用ode45，我完全迷路了。

第一个系统的微分方程：

syms x(t) y(t) 
Dy = diff(y,t); 
Dx = diff(x,t); 
D2x = diff(x,2,t); 
cond = [x(0)==0, y(0)==0, Dy(0)==0, Dx(0)==5]; 
eqn33 = D2x + (2*0.2121*0.1414*Dx) + (0.1414^2)*x==2*0.2121*0.1414*Dy+(0.1414^2)*y; 
sol33 = dsolve(eqn33,cond) 
pretty(sol33)

Answer 1

答案已更新以匹配 Simulink 模型实现

要使用ode45，您首先需要编写一个函数来计算输入向量的导数（即您的微分方程），并将该函数存储在一个单独的文件中，并以函数名作为文件名。请注意，ode 求解器只能求解一阶微分方程，因此您首先需要做一些工作将二阶微分方程转换为一阶微分方程。详情请见documentation on ode45。

根据您在 Simulink 模型中所做的，D2y 已知所有 t 的值（它是步进输入），因此我们需要将其与时间整合以获得 Dy 和 @ 987654329@。所以我们的状态向量是X = [x; Dx; y; Dy]，我们的函数看起来像（存储在diff_eqn.m）：

function dX = diff_eqn(t,X)

m1=500; 
c=1200;
k1=25000; 

dX(1) = X(2); % Dx
dX(2) = -(1/m1)*(c*(X(2)-X(4)/m1) + k1*(X(1)-X(3)/m1));; % D2x
dX(3) = X(4); % Dy

if t<2
  dX(4) = 0; % D2y
else
  dX(4) = 0.5;
end

作为dX = [Dx; D2x; Dy; D2y]。

然后，您可以在脚本或 MATLAB 命令窗口中调用 ode 求解器（根据您的 Simulink 模型，Dx、x、Dy 和 y 的初始条件均为零） ：

[t,X] = ode45(@diff_eqn,[0 20],[0; 0; 0; 0]);

调整ode 求解器选项（例如最大步长等）以获得更多数据点的结果。要获得与 Simulink 模型中相同的图，您可以处理来自 ode 求解器的结果：

D2x = diff(X(:,2))./diff(t);
D2x = [0; D2x];
D2y = zeros(size(D2x));
D2y(t>=2) = 0.5;
plot(t,[D2y 500*D2x])
grid on
xlabel('Time [s]')
legend('D2y','m1*D2x','Location','NorthEast')

它给出了以下图，与您的 Simulink 模型的结果相匹配：