Spring/Damper Calculation & Plotting

Question

Given two systems with damper/spring:

First system's simulink model with step time 2, final value 0.5:

Simulink of the second system with same input:

I have to find the code using dsolve and ode45 to generate the same graph with Simulink. Given values are:

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

I tried to find dsolve but it couldn't solve it. So I got to use ode45, and I am totally lost.

Differential equation of the first system:

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

Answer updated to match Simulink model implementation

To use ode45, you first need to write a function that computes the derivative of you input vector (i.e. your differential equation), and store that function in a separate file with the function name as the filename. Please note that the ode solvers can only solve first-order differential equations, so you first need to do a bit of work to convert your second-order differential equation to a first-order one. For more details, see the documentation on ode45.

Based on what you have done in your Simulink model, D2y is known for all values of t (it's the step input), so we need to integrate it with respect to time to get Dy and y. So our state vector is X = [x; Dx; y; Dy] and our function looks like (stored in 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

as dX = [Dx; D2x; Dy; D2y].

In your script or your MATLAB command window, you can then call the ode solver (initial conditions all being equal to zero for Dx, x, Dy and y, as per your Simulink model):

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

Adjust the ode solver options (e.g. max step size, etc...) to get results with more data points. To get the same plot as in your Simulink model, you can then process the results from the ode solver:

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')

Which gives the following plot, matching the results from your Simulink model: