Simulate shock with spring ODE
显示 更早的评论
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:
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.
13 个评论
darova
2020-3-7
Everything looks correct. One thing:
% your case: u = y(2), du = y(1)
u(1)=y(1);
u(5)=a1*y(5)+K1x;
% usually: u = y(1), du = y(2)
u(1)=y(5);
u(5)=a1*y(1)+K1x;
You know what i mean?
Ameer Hamza
2020-3-7
Actually, it is necessary to use the order as suggested by @darova, otherwise, you will get the wrong answer (as you mentioned exponential curve in your case).
Also, is the 2nd last equation correct? It again mentions 
.
.
darova
2020-3-7
Eugène PLANTEUR
2020-3-8
Ameer Hamza
2020-3-8
If you can provide the values of a1,a2,K1x,K1y,K2x,K2y, then we can try to run the code and find the issue.
darova
2020-3-8
Try to reduce time
time=[0,10];
Eugène PLANTEUR
2020-3-8
Well, one of the final goals of the project is to study those constants (that depend on other constants) would change the phenomenon, and I've been told to first manage to solve the equations with matlab (so using any value for the constants) and then worry about them. But there one major problem, that I thought I'd solve after managing to set up my main code:
There are basically 3 steps in the problem. At the beginning, the first particle moves towards (or not) the second partcile. If, at any given time, the distance between them is lesser or equal to the sum of the radius of the two particles (i.e. they collide), then we have to solve the differential equations at that point in time. But, each K constant depends on the value of the position. For example:
Where k is the stifness constant of the spring, 
the mass of the first particle, 
the equilibrium length of the spring and d the distance bewtween the two particles. So I don't either know how I could compute that :/
the mass of the first particle, the equilibrium length of the spring and d the distance bewtween the two particles. So I don't either know how I could compute that :/
Ameer Hamza
2020-3-8
It appear from your equation tha 
is a function of 
, so you cannot ignore that.
is a function of , so you cannot ignore that.
Ameer Hamza
2020-3-8
The values of parameters are also important. For example, if I run, (negative values of 
and 
).
and ).a1=-1;
a2=-1;
K1x=1;
K1y=1;
K2x=1;
K2y=1;
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);
I don't get exponential curve.
The reason will become clear if you study analytical solution of 2nd order ODEs.
Remember:
function u=odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y)
u=zeros(8,1);
u(1)=y(5);
u(2)=y(6);
u(3)=y(7);
u(4)=y(8);
u(5)=a1*y(1)+K1x;
u(6)=a1*y(2)+K1y;
u(7)=a2*y(3)+K2x;
u(8)=a2*y(4)+K2y;
end
Eugène PLANTEUR
2020-3-8
Ameer Hamza
2020-3-8
Even in that case, you need to include the effect of 
, for example I changed the odefunc as follow
, for example I changed the odefunc as follow function u=odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y)
u=zeros(8,1);
u(1)=y(5);
u(2)=y(6);
u(3)=y(7);
u(4)=y(8);
u(5)=a1*y(1)-K1x*y(3);
u(6)=a1*y(2)-K1y*y(4);
u(7)=a2*y(3)-K2x*y(1);
u(8)=a2*y(4)-K2y*y(2);
end
I lumped the effect of 
into K1x. I made the guess about the values of 
,
and 
. In that case, the system can converge, even with the positive values of a1 and a2
into K1x. I made the guess about the values of , and . In that case, the system can converge, even with the positive values of a1 and a2a1=0.5;
a2=0.5;
K1x=1;
K1y=1;
K2x=1;
K2y=1;
SystemInit=[10,0,10,0,1,2,1,2];
time=[0,100];
[t,y]=ode45(@(t,y) odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y),time,SystemInit);
plot(t,y)
Although it still diverges for some value of initial conditions, i guess that is actually the property of this system of differential equations.
Eugène PLANTEUR
2020-3-8
Ameer Hamza
2020-3-8
Glad to help.
采纳的回答
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
