Plotting a second order system using ode45()

9 次查看（过去 30 天）
Ryan Rizzo
Ryan Rizzo 2018-3-18
评论： Star Strider 2021-12-22
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.

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2018-3-18
I believe ‘dydt(1)’ should simply be:
dydt(1) = y(2);
instead of multiplying it by ‘m’.
  5 个评论
显示 3更早的评论
Philippe Miron
Philippe Miron 2019-11-20
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');

请先登录，再进行评论。

更多回答（1 个）

enrico Rothe
enrico Rothe 2021-12-22
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
  1 个评论
Star Strider
Star Strider 2021-12-22
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by