MATLAB: Simulate shock with spring ODE

odeode45

Hello everyone !
I was given the project to compute and simulate a shock between two particles, using a spring differential equation. It's a 2 dimension problem, so there are eight equations to solve:
Where and are the initial velocities of the first particle (the second one is not moving at the beginning) and are constants. But when I run my code:
SystemInit=[0,0,30,30,2,2,0,0];
time=[0,100];
[t,y]=ode45(@(t,y) odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y),time,SystemInit);
With my odefunc being:
function u=odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y)
u=zeros(8,1);
u(1)=y(1);
u(2)=y(2);
u(3)=y(3);
u(4)=y(4);
u(5)=a1*y(5)+K1x;
u(6)=a1*y(6)+K1y;
u(7)=a2*y(7)+K2x;
u(8)=a2*y(8)+K2y;
end
I get a sort of exponential law for all the components of my y vector (excepted for the two firsts). So am I doing it right ? We've just started learning how to compute differential equations in Matlab, so there might be something I didn't manage to figure out.
I hope I've been clear enough.

Best Answer

So, I managed yo get something right ! Had to rearrange the code a bit, rewrite the equations in a much simpler form (I over-complicated myself), got some help from our teacher and here it is:
m1=1;m2=1;                 %Particles' masses
radius1=0.1;radius2=0.1;    %Particles' radius
k=10000;                    %Spring's stifness parameter
L0=rayon1+rayon2;           %Spring's equilibrium's length
SystemInit=[0,0,        %Initial position of the first particle
            1,0,        %Initial position of the second particle
            1,0.1,      %Initial speed of the first particle
            0,0];       %Initial speed of the second particle
time=[0,10];
[t,y]=ode45(@(t,y) odefonc(t,y,m1,m2,k,L0),time,SystemInit);
plot(y(:,1),y(:,2),'or',y(:,3),y(:,4),'db')
axis equal
And here's the odefunction:
function u=odefonc(t,y,m1,m2,k,L0)
u=zeros(size(y));
u(1)=y(5);              %Velocity of the first particle along the X axis
u(2)=y(6);              %Velocity of the first particle along the Y axis
u(3)=y(7);              %Velocity of the second particle along the X axis
u(4)=y(8);              %Velocity of the second particle along the Y axis
d=sqrt((y(3)-y(1))^2+((y(4)-y(2))^2)); %Distance between the two particles
a1=(k/d*m1)*(L0-d);
a2=(k/d*m2)*(L0-d);
if d<L0     %i.e if there is contact 
u(5)=a1*(y(1)-y(3));    %Acceleration of the first particle along the X axis
u(6)=a1*(y(2)-y(4));    %Acceleration of the first particle along the Y axis
u(7)=a2*(y(3)-y(1));    %Acceleration of the second particle along the X axis
u(8)=a2*(y(4)-y(2));    %Acceleration of the second particle along the Y axis
end
end
What I'm plotting is the trajectory of both particles.
Anyway, thanks to you all, you've helped a lot !
Related Question