1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab code

  1. Apr 17, 2010 #1
    View attachment 25169 1. The problem statement, all variables and given/known data

    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




    2. Relevant 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



    3. 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
    1. The problem statement, all variables and given/known data



    2. Relevant equations



    3. The attempt at a solution
     

    Attached Files:

    Last edited: Apr 17, 2010
  2. jcsd
  3. Apr 17, 2010 #2

    berkeman

    User Avatar

    Staff: Mentor

    Does your code compile and run? What results did you get?
     
  4. Apr 18, 2010 #3
    the code for start could not be compiled
     
  5. Apr 18, 2010 #4

    Mark44

    Staff: 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.
     
  6. Apr 19, 2010 #5
    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
     
  7. Apr 19, 2010 #6

    Mark44

    Staff: Mentor

    Count down to line 61. What is the code in that line?
     
  8. Apr 19, 2010 #7
    FresY1 = -1*(eq1-dist(xPos1,yPos1,pegX,pegY))*k1*sin(theta1) - m1*g;
     
  9. Apr 19, 2010 #8

    Mark44

    Staff: 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.
     
  10. Apr 20, 2010 #9
    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
     
  11. Apr 20, 2010 #10

    Mark44

    Staff: 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?
     
  12. Apr 20, 2010 #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.
     
  13. Apr 20, 2010 #12

    Mark44

    Staff: Mentor

    Last edited by a moderator: Apr 25, 2017
  14. Apr 21, 2010 #13
    Can you please check for me; whether my velocities of masses are correct .There is a sketch in the attatchement for reference
     
  15. Apr 21, 2010 #14

    Mark44

    Staff: Mentor

    As far as I'm concerned, that's really up to you. I'm just helping to get your code running.
     
  16. Apr 24, 2010 #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
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Matlab code
  1. Matlab code needed (Replies: 3)

Loading...