MATLAB: Help with Elastic Collision

circleelastic collisiongraphicshelpincorrect time displayinelasticloopmodelprojectrectanglesimulation

Hi! I am trying to write code to model elastic collisions. Right now the position vs time graphs are completely off (it doesn't even start at t=0). Second it appears to be an inelastic collision right now. Third, there seems to be in acurrate boundary conditions because the object on the right, the red, moves past x=10. I think that is because it is drawn not from a center of a circle but from the left most point of the 4-point "rectangle" with curvature (aka a circle). Any tips?

I have also tried plot (t, xx1) and (t, yy1) , etc within the loop but that also doesn't work.

Interestingly if you take out the line:

    hold on
    plot (xx1, yy1, '-r');
    plot (xx2, yy2, '-b');
    hold off

it has much trouble ploting anything for the Position vs Time graphs.

       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Modeling Elastic Collisions  %%%
%%% Final Project Matlab 225     %%%
%%% K. Gifford August 2020        %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=0;
dt=0.5;
x1=1;
y1=0;
v1x=0;
v1y=0;
m1=1;
% centers = [x1, y1];
x2=-8;
y2=0;
v2x=2;
v2y=0;
m2=1;
% centers2 = [x2, y2];
% t=0, plot position of particles
subplot (3, 1, 1)
plot ([-10, 10], [-10, 10], 'w')
% obj1 = viscircles(centers,1)
% obj2 = viscircles (centers2, 1)
r1 = rectangle ('Position', [x1,y1,1,1], 'FaceColor', 'r', 'Curvature', [1,1]);
r2 = rectangle ('Position', [x2,y2,1,1], 'FaceColor', 'b', 'Curvature', [1,1]);
axis ([-10 10 -10 10]);
subplot (3, 1, 2)
plot (t, x1)
plot (t, x2)
xlabel ('Time')
ylabel ('Position-X')
subplot (3, 1, 3)
plot  (t, y1)
plot  (t, y2)
xlabel ('Time')
ylabel ('Position-Y')
xx1 = [ ];
xx2 = [ ];
yy1 = [ ];
yy2 = [ ];
% Start moving :)
while t < 100
    if 1 < sqrt((x2-x1)^2+(y2-y1)^2)
        x1 =  x1 + (v1x * dt);
        x2 =  x2 + (v2x * dt);  
        y1 =  y1 + (v1y * dt);
        y2 = y2 + (v2y * dt);
        xx1 = [xx1, x1];
        xx2 = [xx2, x2];
        yy1 = [yy1, y1];
        yy2 = [yy2, y2];
       
    else 
        v1x = ((m1-m2)/(m1+m2))*v1x + ((2*m2)/(m1+m2)) * v2x;
        v1y = ((m1-m2)/(m1+m2))*v1y + ((2*m2)/(m1+m2)) * v2y;
        
        v2x = ((m2-m1)/(m1+m2))*v2x + ((2*m1)/(m1+m2)) * v1x;
        v2y = ((m2-m1)/(m1+m2))*v2y + ((2*m1)/(m1+m2)) * v1y;
        x1 =  x1 + (v1x * dt);
        x2 =  x2 + (v2x * dt);  
    
        y1 =  y1 + (v1y * dt);
        y2 =  y2 + (v2y * dt);
        
        
        xx1 = [xx1, x1];
        xx2 = [xx2, x2];
        yy1 = [yy1, y1];
        yy2 = [yy2, y2];
    end
    % At Boundary - Turn Direction
            if x1 ==10
                v1x= -v1x;
            end
            if x2 == 10
                v2x=-v2x;
            end
            if x1 == -10
                v1x=abs(v1x);
            end
            if x2 == -10
                v2x=abs(v2x);
            end
            if y1 == 10
                v1y=-v1y;
            end
            if y2 == 10
                v2y=-v2y;
            end
            if y1 == -10
                v1y=abs(v1y);
            end
            if y2 == -10
                v2y=abs(v2y);
            end
 
    %update the position of the rectangles 
    set(r1,'Position', [x1,y1,1,1]);
    set(r2,'Position', [x2,y2,1,1]);
    
    subplot (3, 1, 2)
    plot (t, x1)
    plot (t, x2)
    xlabel ('Time')
    ylabel ('Position-X')
    
    subplot (3, 1, 3)
    plot (t, y1)
    plot (t, y2)
    xlabel ('Time')
    ylabel ('Position-Y')
    
    hold on
    plot (xx1, yy1, '-r');
    plot (xx2, yy2, '-b');
    hold off
    drawnow
      t= t+ dt;
 end
%

Best Answer

You need to be careful about pre and post collision velocities, and the boundary checks:

t=0;
dt=0.5;
x1=1;
y1=0;
v1x=0;
v1y=0;
m1=1;
% centers = [x1, y1];
x2=-8;
y2=0;
v2x=2;
v2y=0;
m2=1;
% centers2 = [x2, y2];
% t=0, plot position of particles
subplot (3, 1, 1)
plot ([-10, 10], [-10, 10], 'w')
% obj1 = viscircles(centers,1)
% obj2 = viscircles (centers2, 1)
r1 = rectangle ('Position', [x1,y1,1,1], 'FaceColor', 'r', 'Curvature', [1,1]);
r2 = rectangle ('Position', [x2,y2,1,1], 'FaceColor', 'b', 'Curvature', [1,1]);
axis ([-10 10 -10 10]);
subplot (3, 1, 2)
plot (t, x1)
plot (t, x2)
xlabel ('Time')
ylabel ('Position-X')
subplot (3, 1, 3)
plot  (t, y1)
plot  (t, y2)
xlabel ('Time')
ylabel ('Position-Y')
xx1 = [ ];
xx2 = [ ];
yy1 = [ ];
yy2 = [ ];
% Start moving :)
while t < 100
    % Update "old" values
    v1xold = v1x; v2xold = v2x; v1yold = v1y; v2yold = v2y;
    
    if 1 < sqrt((x2-x1)^2+(y2-y1)^2)
        x1 =  x1 + (v1x * dt);
        x2 =  x2 + (v2x * dt);  
        y1 =  y1 + (v1y * dt);
        y2 = y2 + (v2y * dt);
        xx1 = [xx1, x1];
        xx2 = [xx2, x2];
        yy1 = [yy1, y1];
        yy2 = [yy2, y2];
       
    else   % Use "old" values on RHS
        v1x = ((m1-m2)/(m1+m2))*v1xold + ((2*m2)/(m1+m2)) * v2xold; 
        v1y = ((m1-m2)/(m1+m2))*v1yold + ((2*m2)/(m1+m2)) * v2yold;
        
        v2x = ((m2-m1)/(m1+m2))*v2xold + ((2*m1)/(m1+m2)) * v1xold;
        v2y = ((m2-m1)/(m1+m2))*v2yold + ((2*m1)/(m1+m2)) * v1yold;
        x1 =  x1 + (v1x * dt);
        x2 =  x2 + (v2x * dt);  
    
        y1 =  y1 + (v1y * dt);
        y2 =  y2 + (v2y * dt);
        
        
        xx1 = [xx1, x1];
        xx2 = [xx2, x2];
        yy1 = [yy1, y1];
        yy2 = [yy2, y2];
    end
    % At Boundary - Turn Direction
            if x1+dt*v1x>9
                v1x= -v1x;
            end
            if x2+dt*v2x>9
                v2x=-v2x;
            end
            if x1+dt*v1x<-10
                v1x=abs(v1x);
            end
            if x2+dt*v2x<-10
                v2x=abs(v2x);
            end
            if y1+dt*v1y>9
                v1y=-v1y;
            end
            if y2+dt*v2y>9
                v2y=-v2y;
            end
            if y1+dt*v1y<-10
                v1y=abs(v1y);
            end
            if y2+dt*v2y<-10
                v2y=abs(v2y);
            end
            
 
    %update the position of the rectangles 
    set(r1,'Position', [x1,y1,1,1]);
    set(r2,'Position', [x2,y2,1,1]);
    
    subplot (3, 1, 2)
    plot (t, x1,'ro',t,x2,'bs')
    axis([0 100 -10 10])
    xlabel ('Time')
    ylabel ('Position-X')
    
    subplot (3, 1, 3)
    plot (t, y1,'ro',t,y2,'bs')
    axis([0 100 -10 10])
    xlabel ('Time')
    ylabel ('Position-Y')
    drawnow
      t= t+ dt;
 end

Related Solutions

MATLAB: 2nd order numerical differential equation system solving

function main
y0=[0; 5; 1; 1; 1; 4; 1; 1; 2; 3; 1; 1];
t0=0;
tfinal=10;
[T Y] = ode45(@odesNew,[t0 tfinal],y0)
function dy = odesNew(t,y)
G=6.673*10^-11;
m1=1; m2=2; m3=3;
dy=zeros(12,1);
x1=y(1);
x2=y(2);
x3=y(3);
y1=y(4);
y2=y(5);
y3=y(6);
u1=y(7);
u2=y(8);
u3=y(9);
v1=y(10);
v2=y(11);
v3=y(12);
%Körper 1/Mass1
dy(1)=u1;
dy(4)=v1;
dy(7)=(((G*m1*m2)/((x2-x1)^2+(y2-y1)^2)^(3/2))*(x2-x1)+((G*m1*m3)/((x3-x1)^2+(y3-y1)^2)^(3/2))*(x3-x1))/m1;
dy(10)=(((G*m1*m2)/((x2-x1)^2+(y2-y1)^2)^(3/2))*(y2-y1)+((G*m1*m3)/((x3-x1)^2+(y3-y1)^2)^(3/2))*(y3-y1))/m1;
%Körper 2/Mass2
dy(2)=u2;
dy(5)=v2;
dy(8)=(((G*m2*m3)/((x3-x2)^2+(y3-y2)^2)^(3/2))*(x3-x2)+((G*m1*m2)/((x1-x2)^2+(y1-y2)^2)^(3/2))*(x1-x2))/m2;
dy(11)=(((G*m2*m3)/((x3-x2)^2+(y3-y2)^2)^(3/2))*(y3-y2)+((G*m1*m2)/((x1-x2)^2+(y1-y2)^2)^(3/2))*(y1-y2))/m2;
%Körper 3/Mass3
dy(3)=u3;
dy(6)=v3;
dy(9)=(((G*m3*m1)/((x1-x3)^2+(y1-y3)^2)^(3/2))*(x1-x3)+((G*m3*m2)/((x2-x3)^2+(y2-y3)^2)^(3/2))*(x2-x3))/m3;
dy(12)=(((G*m3*m1)/((x1-x3)^2+(y1-y3)^2)^(3/2))*(y1-y3)+((G*m3*m2)/((x2-x3)^2+(y2-y3)^2)^(3/2))*(y2-y3))/m3;

Best wishes

Torsten.

MATLAB: I have encountered a problem and would like to ask everyone to help solve it. The problem to be solved is how to close the following two fitted curves. Thank you very much for your enthusiastic answers.

load dataA.mat
%x1=dataA(:,1);
%y1=dataA(:,2);
f1 = fit(xAbove,yAbove,'smoothingspline','SmoothingParam',0.0001);%smoothingspline
%plot(f1,xAbove,yAbove)
% plot(f1,'-r')
xx1 = xAbove ;
yy1 = f1(xx1) ; 
hold on
load dataB.mat
%x2=dataB(:,1);
%y2=dataB(:,2);
[xx,ind] = sort(xBelow);
yy2 = smooth(xBelow,yBelow,0.5,'rloess');
yy=yy2(ind);
xAbove(1)=xx(1);
yAbove(1)=yy(1);
xx = [xx1 ; flipud(xx) ; xx1(1)] ;
yy = [yy1 ; flipud(yy); yy1(1)] ; 
plot(xx,yy,'r-')

Best Answer

Related Solutions

MATLAB: 2nd order numerical differential equation system solving

MATLAB: I have encountered a problem and would like to ask everyone to help solve it. The problem to be solved is how to close the following two fitted curves. Thank you very much for your enthusiastic answers.

Related Question