Discussion :

[fospdjk]

Bonjour à tous!
Je suis nouveau sur le forum et j'ai lu un peu comment poster des articles, mais veuillez me pardonner si ce n'est pas exactement comme il fallait!
J'utilise depuis peu Matlab à l'université pour modéliser diverses situations simples, mais je suis un peu bloqué sur un problème en ce moment alors je viens demander un petit coup de main!
Voici ce que j'ai écris jusqu'ici:
(j'ai mi en "code" tout le script que j'ai écris, donnée comprise!

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
%Donnée (exerice 12)
%A Particle is traveling with a speed of 1m/sec. Its starting point is at X=0. 
%At X=20 and X=40, there are two other particles that initially do not move. The first particle has
%a radius of 2m and the other ones have a radius of 3m. The last particle is fixed (will never 
%move). The elastic constant of the collision is 100[N]. The mass of every particle is 1kg.
%Take dt = 1e-1. What is the position of the particles after 1 minute? Plot the position of the particles (xaxis) versus time (y-axis).
 
clear,clc,clf
%Physics
X       = 0:20:40;          % initial positions (x) of the particles
R       = [2 3 3];          % radius of the three particles, [m]
m       = [1 1 1];          % mass of the particles, [kg]
V       = [1 0 0];          % velocity of the three particles, [m/s]
F       = [0 0 0];          % Initial forces applied
k       = 100;              % elastic contant, [N]
dist    = R(1)+R(2);        % distance between particle1 and particle2
dist2   = R(2)+R(3);        % distance between particle 2 and particle 3
%Numerics
t       = 0;
dt      = 1e-2;
nt      = 10000;
%Action
for it=1:nt
    dX = abs(diff(X));
    if (dX(dX < dist) < dist)
        E = (dist-dX)./dist;
        F = -k*E;
        a = [-F(1)/m(1) F(2)/m(2) 0];
 
    else
        a = [0 0 0];
    end
 
 
 
    V   = V + a*dt;
    X   = X + V.*dt;
    t = t+dt;
plot(X,t,'o')
xlabel('position')
ylabel('time [s]')
title('Elastic collision')
axis([-5 120 -5 60])
drawnow
end

(J'ai volontairement changé dt et nt afin de voir ce qu'il se passe au ralenti !)
Je vois bien ce qui ne va pas, mais je n'arrive pas à résoudre le problème; la collision entre les 2 premières particules à l'air correcte, seulement la première particule est à nouveau affectée par la force -F(1)/m(1) lorsque la deuxième particule arrive à la 3ème...
J'ai essayé d'écrire ceci :
dX = abs(diff(X));
if (dX(1)(dX(1) < dist) < dist)
E = (dist-dX(1))./dist;
F = -k*E;
a = [-F(1)/m(1) F(2)/m(2) 0];
afin de préciser quelle particule il fallait considérer, puis j'ai fait un else if pour if(dX(2).... afin qu'après le premier choc, seules les particules 2 et 3 s'influencent. Mais il ne se passe plus rien sur mon graphe...

Aussi, pensez-vous que la méthode pour représenter le temps qui passe soit bien ? J'hésite à laisser les particules immobiles sur l'axe Y et représenter une particule sur l'axe Y justement qui monterait verticalement en fonction du temps qui passe...
Lorsque je tape "E" dans la console, il n'y a qu'une seule valeur retournée alors qu'il devrait y en avoir deux... Dois-je créer un tableau pour les mettre ?
Navré si les questions sont un peu élémentaires ou si ça n'a pas vraiment de sens, mais j'ai pas mal de peine à bien tout comprendre pour le moment...
En vous remerciant de votre aide, je vous souhaite une belle journée !

Fospdjk

[le fab]

Salut

déjà
tu aura peut être remarqué que dX est un vecteur à deux éléments
il ne contiendra les distances qu'entre les particules 1 et 2 d'une part et 2 et 3 d'autre part
1 et 3 ne peuvent interagir entre elles
pour qu'elle le puisse tu pourrait faire un dX = abs(diff([X X(1)])) peut être

ensuite
ton problème vient de a = [-F(1)/m(1) F(2)/m(2) 0]; : la particule 3 ne sera jamais impacté par les autres

enfin tu devrais :
- créer ton graph avant la boucle avec la fonction figure
- créer les titres et étiquettes des axes hors de la boucle (avant)
- récupérer le handle de ton plot afin de modifier les données de celui ci plutôt que de retracer, comme ici
cela rendrait ton animation plus fluide

Fabien