% this program plots the displacement of each mass. positive displacement
% means that the objects is moved to the right, negative means to the left.
nsprings=20; % number of masses (numnber of springs is one less
% than the number of masses)
x=1:nsprings+1; % initial position of all of the masses.
vx=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
a=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
% initial velocities are all zero and so are accelerations
k=[.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1];
% spring constant of each spring
% spring "0" is the spring connecting mass(1) [the wall] and mass(2),
% and so on.
% masses are all assumed to be 1.
dt=0.01; % is this a good time step?
x(2)=2.5; % <=== here is the initial condition!
%initial condition is that one mass is stretched to a new position.
displacement=x-(1:21) % calculate displacement from initial positions
h=plot(displacement);
xlabel('mass number')
ylabel('displacement')
axis([1 21 -1 1])
set(h,'EraseMode','none','LineWidth',3)
while 1 % run the simulation until you hit cntrl-C
% first, calculate displacement of each spring on the left
% of each mass
stretch=(x(2:nsprings)-x(1:nsprings-1))-1;
% and the acceleration from Hooke's law
a(2:nsprings)=-k.*stretch; % note the matrix multiply!
% but the wall is fixed, so there can't be any acceleration at the wall.
a(nsprings)=0;
% now look at the spring to the right of each mass and add that acceleration
stretch=(x(1:nsprings-1)-x(2:nsprings))+1;
a(1:nsprings-1)=a(1:nsprings-1)-k.*stretch;
a(1)=0; % the left wall doesn't move either
oldx=x;
% calculate new positions and velocities
x=.5*a.*dt^2+vx.*dt+x; % a little simple physics here!
vx=vx+a.*dt;
displacement=x-(1:21);
set(h,'YData',displacement,'Color',[0,0,1]) % overplot a blue line
set(h,'Color',[1,0,0]) % overplot a red line
drawnow % plot it!
end