Plotting a second order system using ode45()

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.

Answer 1

[Star Strider]

I believe ‘dydt(1)’ should simply be:

dydt(1) = y(2);

instead of multiplying it by ‘m’.

[Ryan Rizzo]

Many thanks for your answer.

I am not sure about this. The second order derivative has a variable 'm' tied to it which represents the mass. Should this be ignored in this case?

Either way, the displacement still settles above zero.

[Star Strider]

My pleasure.

Just to be sure, I asked the Symbolic Math Toolbox:

syms c g k m x(t) 
Eqn = m*diff(x,2) + c*diff(x) + k*x == m*g;
[DE. Subs] = odeToVectorField(Eqn, x(0) == 0)

It replied:

DE =
                       Y[2]
 -(c*Y[2] - g*m + k*Y[1])/m

with ‘Y[1]=x’, and ‘Y[2]=Dx’.

[Ryan Rizzo]

I see, so there must be another mistake form my end since I do not think the outputted graphs make sense. Thanks once again for pointing that out.

[Star Strider]

My pleasure.

[Philippe Miron]

The other error is simply the legend of the plot. y(:,1) contains the position and y(:,2) contains the velocity. So the legend arguments:

legend('Velocity', 'Displacement');

should be:

legend('Displacement', 'Velocity');

Answer 2

[enrico Rothe]

hey,

can someone explain to me why i get the error if i copy paste the code :

Not enough input arguments.

Error in springmassdamper (line 4)

dydt(1)=y(2);

i m new to matlab and appreciate any help

[Star Strider]

Somewhere on the system running that code, there is a function defined as ‘y’ with more than one argument.

Run this from a script or the Command Window —

which y -all

That should reveal the problem. Alternatively, there could be a ‘dydt’ function, however I doub thatt it would throw the same error in this instance. It could be worthwhile checking it as well. The solution is to give the function another name.