Simulating System of Masses Connected by Springs

  • Thread starter Thread starter 200917717
  • Start date Start date
  • Tags Tags
    Springs System
AI Thread Summary
The discussion focuses on simulating a system of masses connected by springs, detailing the parameters for each mass and spring, including their positions, masses, spring constants, and equilibrium distances. The simulation aims to track the movement of the masses over 100 seconds in 0.1-second intervals, requiring the calculation of forces, accelerations, and velocities based on spring dynamics without considering friction or gravity. Participants encountered issues with code compilation and plotting, particularly in ensuring that all masses' paths were displayed correctly in different colors on a single plot. Suggestions were made to modify the plotting commands and ensure proper force calculations for accurate simulation results. The conversation emphasizes the importance of debugging and refining the code to achieve the desired output.
200917717
Messages
8
Reaction score
0
https://www.physicsforums.com/attachments/251691. Homework Statement

assigned system of masses connected by springs with the parameters specified
Mass M0 details:
Stationary mass, does not move
Starting Positions in meters:
X: 5.8
Y: 4.2
-----------------------------------------
Mass M1 details:
Mass 14 kg
Starting Positions in meters:
X: 5.8
Y: 3
-----------------------------------------
Mass M2 details:
Mass 16 kg
Starting Positions in meters:
X: 5.4
Y: 5.6
-----------------------------------------
Mass M3 details:
Mass 14 kg
Starting Positions in meters:
X: 4.2
Y: 4.8
-----------------------------------------
#####################################################################
Spring S0 details:
K: 0.78
Equilibrium Distance (meters):
Eq: 0.9
-----------------------------------------
Spring S1 details:
K: 0.7
Equilibrium Distance (meters):
Eq: 1.6
-----------------------------------------
Spring S2 details:
K: 1.02
Equilibrium Distance (meters):
Eq: 1.3
-----------------------------------------
Spring S3 details:
K: 0.78
Equilibrium Distance (meters):
Eq: 1
-----------------------------------------
Spring S4 details:
K: 0.66
Equilibrium Distance (meters):
Eq: 0.8
Using model it must simulate the movement of all the masses from 0 to 100 seconds (including 100) with a time step of 0.1 seconds. Record the position of each mass (X and Y position) at every time step in an array. You will then plot the path of each mass in a different colour on a single plot with the X Axis to run from 0 to 10 meters and the Y axis from 0 to 7.5 meters.No friction or gravity




Homework Equations



F=-kx where k is the spring constant and x is the displacement from the springs equilibrium point.
Fx=Fcos(θ) where Theta is the angle at which the spring is pulling.
Fy=Fsin(θ)
You add up all the X component forces and Y component forces for a specific mass to find the result X and Y component force acting on that mass. Once you have the resultant force acting on a specific mass you can find its Acceleration from the following:
Fres = ma thus a = m/Fres
a=∆v/∆t thus ∆V=a∆t
∆V= ∆s/∆t thus ∆s= ∆v∆t



The Attempt at a Solution


code for start
clear all;
close all;
clc;

%position of peg
pegX = 5.8;
pegY = 4,2;

%position of mass
xPos1 = 5.8;
yPos1 = 3;
xPos2 = 5.4;
yPos2 = 5.6;
xPos3 =4.2 ;
yPos3 = 4.8;

%velocity of mass
xVel1 = 0;
yVel1 = 2;
xVel2 = 1.67;
yVel2 = 2.79;
xVel3 = -0.64;
yVel3 = -3.19;
xVel4 = -6;
yVel4 = 5;
%mass of mass
m1 = 14;
m2 = 16;
m3 = 14;
%constant of spring
kpeg = 0.78;
k1 = 0.7;
k2 = 1.02;
k3 = 0.78;
%equilibrium distance of spring
egpeg = 0.9;
eq1 = 1.6;
eq2 = 1.3
eq3 = 1;
%constant of gravity
g = 0.5;

%define time domain
dt = 0.1;
t = 0:dt:100;

%record all previos positions for plot
positions = zeros(4,length(t));
for i = 1:length(t)
%find angle between peg and mass
theta1 = direction(xPos1,yPos1,pegX,pegY);
theta2 = direction(xPos3,yPos3,xPos1,yPos1);
theta3 = direction(xPos2,yPos2,xPos3,yPos3);
theta4 = direction(xPos1,yPos1,xPos2,yPos2);
theta5 = direction(pegX,pegY,xPos2,yPos2);
%calc X resultant force
FresX1 = -1*(eq1-dist(xPos1,yPos1,pegX,pegY))*k1*cos(theta1);
FresX2 = -1*(eq2-dist(xPos2,yPos2,xPos1,yPos1))*k2*cos(theta2);
FresX3 = -1*(eq3-dist(xPos3,yPos3,xPos2,yPos2))*k2*cos(theta3);
%calc Y resultant force
FresY1 = -1*(eq1-dist(xPos1,yPos1,pegX,pegY))*k1*sin(theta1) - m1*g;`clear all;
FresY3 = -1*(eq3-dist(xPos3,yPos3,xPos1,yPos1))*k3*sin(theta2) - m3*g;`clear all;
FresY2 = -1*(eq2-dist(xPos2,yPos,xPos3,yPos3))*k2*sin(theta3) - m2*g;`clear all;
FresY4 = -1*(eq1-dist(xPos1,yPos1,xPos2,yPos2))*k1*sin(theta4) - m1*g;`clear all;
FresY5 = -1*(eqpeg-dist(xPospegx,yPospegy,xPos2,yPos2))*kpeg*sin(theta5) - mpeg*g;`clear all;
%calculate acceleration due to force
aX1 = FresX1/m1;
aY1 = FresY1/m1;
aX2 = FresX2/m2;
aY2 = FresY2/m2;
aX3 = FresX3/m3;
aY3 = FresY3/m2;
%change velocity based on acceleration
xVel1 = xVel1 + aX1*dt;
yVel1 = yVel1 + aY1*dt;
xVel2 = xVel2 + aX2*dt;
yVel2 = yVel2 + aY2*dt;
xVel3 = xVel3 + aX3*dt;
yVel3 = yVel3 + aY3*dt
%change position based on velocity
xPos1 = xPos1 + xVel1*dt;
yPos1 = yPos1 + yVel1*dt;
xPos2 = xPos2 + xVel2*dt;
yPos2 = yPos2 + yVel2*dt;
xPos3 = xPos3 + xVel3*dt;
yPos2 = yPos3 + yVel3*dt;
%record position for plot
positions(1,i) = xPos1;
positions(2,i) = yPos1;
positions(3,i) = xPos2;
positions(4,i) = yPos2;
positions(5,i) = xPos3;
positions(6,i) = yPos3;
end
%plot position
p1 = plot(positions(1,:),positions(2,:),'r');
p2 = plot(positions(3,:),positions(4,:),'b');
p2 = plot(positions(5,:),positions(6,:),'g');
grid on;
Code for direction:

function [ a ] = direction(x1,y1,x2,y2 )
%DIRECTION find direction of line between the 2 points
a = atan2(y2-y1,x2-x1);
end
code for distance:
function [ d] = dist( x1,y1,x2,y2)
%DIST calculate euclidian distance between 2 points
d = sqrt((x1-x2)^2+(y1-y2)^2);
end
 

Attachments

Last edited:
Physics news on Phys.org
Does your code compile and run? What results did you get?
 
berkeman said:
Does your code compile and run? What results did you get?
the code for start could not be compiled
 
Did you get any error messages? Help us out here. Don't make us pick through 50 or so lines of code to see what the problem might be.
 
yes,These was the first errow"? Error: File: E:\prac.m Line: 61 Column: 73
Missing variable or function."I don't know if there are more
 
Count down to line 61. What is the code in that line?
 
Mark44 said:
Count down to line 61. What is the code in that line?

FresY1 = -1*(eq1-dist(xPos1,yPos1,pegX,pegY))*k1*sin(theta1) - m1*g;
 
That's not all that's on that line.
FresY1 = -1*(eq1-dist(xPos1,yPos1,pegX,pegY))*k1*sin(theta1) - m1*g; `clear all;[/color]
I believe that part in red is causing your problem. Use % for comments. There's another form of comment that uses ' followed by an ellipsis (...), but you don't have the three dots.
 
I tried the following code,it compiled but does not display the accurate plot:
clear all;
close all;
clc;

%position of peg
pegX = 5.8;
pegY = 4,2;

%position of mass
xPos1 = 5.8;
yPos1 = 3;
xPos2 = 5.4;
yPos2 = 5.6;
xPos3 =4.2 ;
yPos3 = 4.8;

%velocity of mass
xVel1 = 0;
yVel1 = 2;
xVel2 = 1.67;
yVel2 = 2.79;
xVel3 = -0.64;
yVel3 = -3.19;
xVel4 = -6;
yVel4 = 5;
%mass of mass
m1 = 14;
m2 = 16;
m3 = 14;
%constant of spring
kpeg = 0.78;
k1 = 0.7;
k2 = 1.02;
k3 = 0.78;
%equilibrium distance of spring
egpeg = 0.9;
eq1 = 1.6;
eq2 = 1.3
eq3 = 1;

%define time domain
dt = 0.1;
t = 0:dt:100;

%record all previos positions for plot
positions = zeros(4,length(t));
for i = 1:length(t)
%find angle between peg and mass
theta1 = direction(xPos1,yPos1,pegX,pegY)+direction(xPos1,yPos1,peg3,peg3)+direction(xPos1,yPos1,peg2,peg2);
theta2 = direction(xPos3,yPos3,xPos1,yPos1)+direction(xPos3,yPos3,xPos2,yPos2);
theta3 = direction(xPos2,yPos2,xPos3,yPos3)+ direction(xPos2,yPos2,xPos1,yPos1)+direction(xPos2,yPos2,pegX,pegy);
%calc X resultant force
FresX1 = -1*(eq1-dist(xPos1,yPos1,pegX,pegY))*k1*cos(theta1);
FresX2 = -1*(eq2-dist(xPos2,yPos2,xPos1,yPos1))*k2*cos(theta2);
FresX3 = -1*(eq3-dist(xPos3,yPos3,xPos2,yPos2))*k2*cos(theta3);
%calc Y resultant force
FresY1 = -1*(eq1-dist(xPos1,yPos1,pegX,pegY))*k1*sin(theta1);
FresY3 = -1*(eq3-dist(xPos3,yPos3,xPos1,yPos1))*k3*sin(theta2);
FresY2 = -1*(eq2-dist(xPos2,yPos2,xPos3,yPos3))*k2*sin(theta3);
FresY4 = -1*(eq1-dist(xPos1,yPos1,xPos2,yPos2))*k1*sin(theta3);

%calculate acceleration due to force
aX1 = FresX1/m1;
aY1 = FresY1/m1;
aX2 = FresX2/m2;
aY2 = FresY2/m2;
aX3 = FresX3/m3;
aY3 = FresY3/m3;
%change velocity based on acceleration
xVel1 = xVel1 + aX1*dt;
yVel1 = yVel1 + aY1*dt;
xVel2 = xVel2 + aX2*dt;
yVel2 = yVel2 + aY2*dt;
xVel3 = xVel3 + aX3*dt;
yVel3 = yVel3 + aY3*dt
%change position based on velocity
xPos1 = xPos1 + xVel1*dt;
yPos1 = yPos1 + yVel1*dt;
xPos2 = xPos2 + xVel2*dt;
yPos2 = yPos2 + yVel2*dt;
xPos3 = xPos3 + xVel3*dt;
yPos3 = yPos3 + yVel3*dt;
%record position for plot
positions(1,i) = xPos1;
positions(2,i) = yPos1;
positions(3,i) = xPos2;
positions(4,i) = yPos2;
positions(5,i) = xPos3;
positions(6,i) = yPos3;
end
%plot position
p1 = plot(positions(1,:),positions(2,:),'r');
p2 = plot(positions(3,:),positions(4,:),'g');
p3 = plot(positions(5,:),positions(6,:),'r');
grid on
 
  • #10
So at least now your code compiles. Can you be more specific about what you mean when you say it doesn't display the accurate plot? You're calling the plot function three times. Are all three not displaying correctly? What does each one display?
 
  • #11
programme does not plot the paths of each mass on a single plot in different colour,it just plot a sigle line for the last plot 3 only.
 
  • #12
Last edited by a moderator:
  • #13
Can you please check for me; whether my velocities of masses are correct .There is a sketch in the attatchement for reference
 
  • #14
As far as I'm concerned, that's really up to you. I'm just helping to get your code running.
 
  • #15
I tried the these code,but it plots for mass3 only
function [c]=tall
clear all;
close all;
clc;

%position of peg
pegX = 5.8;
pegY = 4.2;

%position of mass
xPosm1 = 5.8;
yPosm1 = 3;
xPosm2 = 5.4;
yPosm2 = 5.6;
xPosm3 =4.2 ;
yPosm3 = 4.8;

%velocity of mass
xVelm1 = 0;
yVelm1 = 0;
xVelm2 = 0;
yVelm2 = 0;
xVelm3 = 0;
yVelm3 = 0;
%mass of mass
m1 = 14;
m2 = 16;
m3 = 14;
%constant of spring
kpeg = 0.78;
k1 = 0.7;
k2 = 1.02;
k3 = 0.78;
k4 = 0.66;
%equilibrium distance of spring
eqpeg = 0.9;
eq1 = 1.6;
eq2 = 1.3;
eq3 = 1;
eq4 = 0.8;

%define time domain
dt = 0.1;
t = 0:dt:100;

%record all previos positions for plot
positionsm1= zeros(2,length(t));
positionsm2= zeros(2,length(t));
positionsm3= zeros(2,length(t));

for i = 1:length(t)
%find angle between peg, mass1, mass2 and mass3
thetam1peg = direction(xPosm1,yPosm1,pegX,pegY);
thetam13 = direction(xPosm1,yPosm1,xPosm3,yPosm3);
thetam12 = direction(xPosm1,yPosm1,xPosm2,yPosm2);
%find angle between peg, mas1, mass2 and mass3
thetam2peg = direction(xPosm2,yPosm2,pegX,pegY);
thetam21 =direction(xPosm2,yPosm2,xPosm1,yPosm1);
thetam23 = direction(xPosm2,yPosm2,xPosm3,yPosm3);
%find angle between mas1,mass2 and mass3
thetam31 = direction(xPosm3,yPosm2,xPosm1,yPosm1);
thetam32 =direction(xPosm3,yPosm2,xPosm2,yPosm2);
%calc X resultant force of m1
FresXm1 = -1*(eqpeg-dist(xPosm1,yPosm1,pegX,pegY)*kpeg*cos(thetam1peg))-1*(eq4-dist(xPosm1,yPosm1,xPosm3,yPosm3)*k4*cos(thetam13)) -1*(eq2-dist(xPosm1,yPosm1,xPosm2,yPosm2)*k2*cos(thetam12));
%calc X resultant force of m2
FresXm2 = -1*(eq1-dist(xPosm2,yPosm2,pegX,pegY)*k1*cos(thetam2peg))-1*(eq2-dist(xPosm2,yPosm2,xPosm1,yPosm1)*k2*cos(thetam21)) -1*(eq3-dist(xPosm2,yPosm2,xPosm3,yPosm3)*k3*cos(thetam23));
%calc X resultant force of m3
FresXm3 = -1*(eq4-dist(xPosm3,yPosm3,xPosm1,yPosm1)*k4*cos(thetam31))-1*(eq3-dist(xPosm3,yPosm3,xPosm2,yPosm2)*k3*cos(thetam32));
%calc Y resultant force of m1
FresYm1 = -1*(eqpeg-dist(xPosm1,yPosm1,pegX,pegY)*kpeg*sin(thetam1peg))-1*(eq2-dist(xPosm1,yPosm1,xPosm2,yPosm2)*k2*sin(thetam12)) -1*(eq4-dist(xPosm1,yPosm1,xPosm3,yPosm3)*k4*sin(thetam13));
%calc Y resultant force of m2
FresYm2 = -1*(eq1-dist(xPosm2,yPosm2,pegX,pegY)*k1*sin(thetam2peg))-1*(eq2-dist(xPosm2,yPosm2,xPosm1,yPosm1)*k2*sin(thetam21)) -1*(eq3-dist(xPosm2,yPosm2,xPosm3,yPosm3)*k3*sin(thetam23));
%calc Y component force of m3
FresYm3 = -1*(eq4-dist(xPosm3,yPosm3,xPosm1,yPosm1)*k4*sin(thetam31))-1*(eq3-dist(xPosm3,yPosm3,xPosm2,yPosm2)*k3*sin(thetam32));
%calculate acceleration due to force
aXm1 = FresXm1/m1;
aYm1 = FresYm1/m1;
aXm2 = FresXm2/m2;
aYm2 = FresYm2/m2;
aXm3 = FresXm3/m3;
aYm3 = FresYm3/m3;
%change velocity based on acceleration
xVelm1 = xVelm1 + aXm1*dt;
yVelm1 = yVelm1 + aYm1*dt;
xVelm2 = xVelm2 + aXm2*dt;
yVelm2 = yVelm2 + aYm2*dt;
xVelm3 = xVelm3 + aXm3*dt;
yVelm3 = yVelm3 + aYm3*dt
%change position based on velocity
xPosm1 = xPosm1 + xVelm1*dt;
yPosm1 = yPosm1 + yVelm1*dt;
xPosm2 = xPosm2 + xVelm2*dt;
yPosm2 = yPosm2 + yVelm2*dt;
xPosm3 = xPosm3 + xVelm3*dt;
yPosm3 = yPosm3 + yVelm3*dt;
%record position for plot
positionsm1(1,i) = xPosm1;
positionsm1(2,i) = yPosm1;
positionsm2(3,i) = xPosm2;
positionsm2(4,i) = yPosm2;
positionsm3(5,i) = xPosm3;
positionsm3(6,i) = yPosm3;
end
%plot position
p1 = plot(positionsm1(1,:),positionsm1(2,:),'r');
hold all
p2 = plot(positionsm2(3,:),positionsm2(4,:),'g');
hold all
p3 = plot(positionsm3(5,:),positionsm3(6,:),'b');
grid on
 

Similar threads

Distribution of Energy when work is done on a system of 2 masses connected by a spring
Replies
17
Views
1K
Two masses connected by spring, find period of oscillation
Replies
3
Views
5K

Hot Threads

Recent Insights

Back
Top