MATLAB: Plotting a second order system using ode45()

differential equationsfunctionMATLABode45plotplotting

The system I am modelling is a spring-mass-damper where m is mass, k is spring stiffness and c is the damping coefficient.
First of all I have a basic function splitting the second order differential equation:
function dydt=springmassdamper(t,y,m,c,k)
dydt=zeros(2,1);
dydt(1)=m*y(2);
dydt(2)=-c*y(2)-k*y(1) + m*9.81;
Next I am allowing user input and plotting using ode45() as follows:
clear all;
close all;
clc;
disp('Spring-Mass-Damper Calculator');
try
    mass = input('Enter mass m (kg) between 0.5kg and 20kg ');
    if (mass >=0.5 && mass <= 20)
        m = mass;
    else
        mass = input('Enter mass m (kg) between 0.5kg and 20kg');
        m = mass;
    end
      spring_damping = input('Enter spring damping c (Ns/m)');
      if (spring_damping > 0)
          c = spring_damping;
      else
          disp("Enter a value larger than 0");
      end
      spring_stiffness = input('Enter spring stifness k (N/m)');
      if (spring_stiffness > 0)
          k = spring_stiffness;
      else
          disp("Enter a value larger than 0");
      end
%initial conditions [x(0) and x*(0)] [displacement velocity]


    [t,y] = ode45(@(t,y) springmassdamper(t,y,m,c,k), [0 1], [0 0]);
      plot(t,y(:,1),'r','LineWidth',1);
      hold on;
      plot(t,y(:,2),'k','LineWidth',1);
      legend('Velocity', 'Displacement');
      xlabel('Time (s)');
      ylabel('Displacement (m)');
      title('Spring-Mass-Damper');
catch ME
warning('Can not use characters, please enter a number instead');  
end
Thus, with 0 initial conditions for both velocity and displacement I get the following graph which does not make sense to me:
The issue I am encountering is that, when I have any type of input, one of the graphs (displacement) does not settle at zero. I do not understand how this is possible. Is there an issue with the parameters I am plotting or is there something wrong in one of the functions?
Example User input:
Enter mass m (kg) between 0.5kg and 20kg: 10, Enter spring damping c (Ns/m): 50, Enter spring stifness k (N/m): 1000
The red graph does not seem correct to me. It looks like I am exceeding Hooke's Law and deforming the spring.
I also notice that if I increase the velocity initial conditions enough such that
    %initial conditions [x(0) and x*(0)] [displacement velocity]
    [t,y] = ode45(@(t,y) springmassdamper(t,y,m,c,k), [0 1], [0 15]);
the red graph (displacement) approaches zero. With the exact same user Input I get
When I input initial conditions for displacement to be 10 and velocity at 0, I get the mass settling at 0 displacement which is correct.
      %initial conditions [x(0) and x*(0)] [displacement velocity]
      [t,y] = ode45(@(t,y) springmassdamper(t,y,m,c,k), [0 1], [10 0]);
Any suggestions on what I am doing wrong and/or flaws in my reasoning would be appreciated. I must be missing something crucial here.

Best Answer

I believe ‘dydt(1)’ should simply be:
dydt(1) = y(2);
instead of multiplying it by ‘m’.
Related Question