MatLab中的物理模拟

Question

我真的停留在这段代码上。我试图模拟一个物理问题来想象一个问题。它有2个挂钩和两个弹簧连接到x方向有50N横风的质量块。我不太清楚如何在x和y中合并限制。请任何帮助

%position of peg1 
pegX1 = 3; 
pegY1 = 10; 

%position of peg2 
pegX2 = 8; 
pegY2 = 10; 

%position of mass 
xPos = 5; 
yPos = 6; 

%velocity of mass 
xVel = 0; 
yVel = 0; 

%mass of mass 
m = 2; 

%constant of springs 
k1 = 20; 
k2 = 25; 

%equilibrium distance of spring 
eq1 = 2; 
eq2 = 1.5; 

%constant of gravity 
g = 5; 

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

%record all previous positions for plot 
positions = zeros(2,length(t)); 


for i = 1:length(t) 
    %find angle between peg1 and mass 
    theta1 = direction(xPos,yPos,pegX1,pegY1); 

    %find angle between peg2 and mass 
    theta2 = direction(xPos,yPos,pegX2,pegY2); 

    %calc X resultant force 
    FresX = -1*(eq1-dist(xPos,yPos,pegX1,pegY1,pegX2,pegY2))*k1*k2*cos(theta1+theta2); 

    %calc Y resultant force 
    FresY = -1*(eq1-dist(xPos,yPos,pegX1,pegY1,pegX2,pegY2))*k1*k2*sin(theta1+theta2) - m*g; 

    %calculate acceleration due to force 
    aX = FresX/m; 
    aY = FresY/m; 

    %change velocity based on acceleration 
    xVel = xVel + aX*dt; 
    yVel = yVel + aY*dt; 

    %change position based on velocity 
    xPos = xPos + xVel*dt; 
    yPos = yPos + yVel*dt; 

    %record position for plot 
    positions(1,i) = xPos; 
    positions(2,i) = yPos; 

    %plot position 
    p1 = plot(positions(1,1:i),positions(2,1:i),'r'); 
    axis([0 4.5 0 4.5]); 
    grid on 
    pause(0.1) 
end 

Answer 1

我不明白你的意思是什么限制。

我改变了定义合力的方式。

function m = test 
%position of peg1 
pegX1 = 3; 
pegY1 = 10; 

%position of peg2 
pegX2 = 8; 
pegY2 = 10; 

%position of mass 
xPos = 5; 
yPos = 6; 

%velocity of mass 
xVel = 0; 
yVel = 0; 

%mass of mass 
m = 2; 

%constant of springs 
k1 = 20; 
k2 = 25; 

%equilibrium distance of spring 
eq1 = 2; 
eq2 = 1.5; 

FWind = 50; 

%constant of gravity 
g = 5; 

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

%record all previous positions for plot 
positions = zeros(2,length(t)); 

figure(1) 
i = 1; 
%title([num2str(0), ' s']) 
plot(xPos, yPos, 'bo', pegX1, pegY1, 'go', pegX2, pegY2, 'ro'); 
drawnow 
d = max([distance(xPos,yPos,pegX1,pegY1), distance(xPos,yPos,pegX2,pegY2)]); 
axis([min([pegX1, pegX2])-d, max([pegX1, pegX2])+d, min([pegY1, pegY2])-d, max([pegY1, pegY2])+d]) 
ax = gca; 
ax.NextPlot = 'replaceChildren'; 

positions(1,i) = xPos; 
positions(2,i) = yPos; 

Frame(i) = getframe; 

    for i = 2:length(t) 
     %find angle between peg1 and mass 
     theta1 = direction(xPos,yPos,pegX1,pegY1); 

     %find angle between peg2 and mass 
     theta2 = direction(xPos,yPos,pegX2,pegY2); 

     %calc X resultant force 
    % FresX = -1*(eq1-dist(xPos,yPos,pegX1,pegY1,pegX2,pegY2))*k1*k2*cos(theta1+theta2); 

     FX = k1 * (distance(xPos,yPos,pegX1,pegY1) - eq1) * cos(theta1) + ... 
       k2 * (distance(xPos,yPos,pegX2,pegY2) - eq2) * cos(theta2) + ... 
       FWind; 

     %calc Y resultant force 
    % FresY = -1*(eq1-dist(xPos,yPos,pegX1,pegY1,pegX2,pegY2))*k1*k2*sin(theta1+theta2) - m*g; 

     FY = k1 * (distance(xPos,yPos,pegX1,pegY1) - eq1) * sin(theta1) + ... 
       k2 * (distance(xPos,yPos,pegX2,pegY2) - eq2) * sin(theta2) ... 
       -m * g; 

     %calculate acceleration due to force 
     aX = FX/m; 
     aY = FY/m; 

     %change velocity based on acceleration 
     xVel = xVel + aX*dt; 
     yVel = yVel + aY*dt; 

     %change position based on velocity 
     xPos = xPos + xVel*dt; 
     yPos = yPos + yVel*dt; 

     %record position for plot 
     positions(1,i) = xPos; 
     positions(2,i) = yPos; 

%  title([num2str(t(i)), ' s']) 
     plot(xPos, yPos, 'bo', pegX1, pegY1, 'go', pegX2, pegY2, 'ro'); 
     drawnow 
     Frame(i) = getframe; 
    end 

    %plot position 
    figure(2) 
    p1 = plot(positions(1,:),positions(2,:),'r'); 
    axis('tight'); 
    grid on 

    figure(1) 
    movie(Frame); 

end 

function theta = direction(x1,y1,x2,y2) 
    theta = atan2(y2 - y1, x2 - x1);  
end 

function d = distance(x1, y1, x2, y2) 
    x = x2 - x1; 
    y = y2 - y1; 

    d = sqrt(x^2 + y^2); 
end