# Matlab code

View attachment 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

## The Attempt at a Solution

#### Attachments

• 300.6 KB Views: 199
Last edited:

Related Engineering and Comp Sci Homework Help News on Phys.org
berkeman
Mentor
Does your code compile and run? What results did you get?

Does your code compile and run? What results did you get?
the code for start could not be compiled

Mark44
Mentor
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 dont know if there are more

Mark44
Mentor
Count down to line 61. What is the code in that line?

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;

Mark44
Mentor
That's not all that's on that line.
FresY1 = -1*(eq1-dist(xPos1,yPos1,pegX,pegY))*k1*sin(theta1) - m1*g; clear all;
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

Mark44
Mentor
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?

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.

Mark44
Mentor
Last edited by a moderator:
Can you please check for me; whether my velocities of masses are correct .There is a sketch in the attatchement for reference

Mark44
Mentor
As far as I'm concerned, that's really up to you. I'm just helping to get your code running.

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