Discussion :

[.Ophiuchus.]

Bonsoir tout le monde,
Etudiants en deuxieme annee de physique nous devons ecrire un programme sur Octave permettant de calculer la période d'un oscillateur voilà une ébauche de programme mais il ne fonctionne pas cause de nombreuses erreurs sauf que je ne vois pas ou elles sont votre aide est la bienvenue, merci d'avance!
"

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
clear all;
clf;
% Initialisation of alpha and beta
alpha =1.0;
beta=0.1;
% Input the initial position
x1=input('Enter the initial value of the displacement in meters:');
% set up the initial velocity, mass and time step
Nstep = 601
v1=0; % Initial velocity in meters/s.
mass = 1.0 % Mass of particle
step = 0.02 %time step in seconds
 
for istep = 1: Nstep % make the number of steps necessary
t1 (istep) = (istep -1) * step; % time
Force1 = alpha * x1 - 4 * beta * x1^3; % compute the force on the particle
accel1 = Force1 / mass; % compute the accel. on the particle
% compute the velocity and position using Euler algo.
x1 = position1(istep) = x1 + v1 * step; 
v1 = vitesse1(istep) = v1 + accel1 * step; 
 
if (vitesse1(istep) = 0)
x=position1(istep);
disp(x)
t1 (istep) = t;
disp(t)
endif
endfor
 
Jstep=istep+0.02;
for jstep = Jstep:Nstep
t2 (jstep) = (jstep -1) * step; % time
Force2 = alpha * x2 - 4 * beta * x2^3; % compute the force on the particle
accel2 = Force2 / mass; % compute the accel. on the particle
x2 = position2(jstep) = x2 + v2 * step; 
v2 = vitesse2(jstep) = v2 + accel2 * step;
 
 
if (position2(jstep)==x)
Periode=t2 (jstep) - t;
disp (Periode)
endif
endfor

"

[tubaas]

salut
- en Matlab il n'y a pas de endif, endfor, simplement end
- le signe = (en programmation en général) n'a rien à voir avec le signe mathématique mais correspond à l'assignation de variable : tu ne peux pas écrire a=b=c
- dans les tests conditionnels il faut utiliser == pour l'égalité
- si tu utilises clf alors qu'aucune fenêtre n'est ouverte ça en ouvre une ce qui est pas top
Sachant que tu n'as pas l'air de te servir des variables position et vitesse voici ce que je te propose (pas d'erreur mais je ne pense pas que les résultats soient bons, il y a des lignes vraiment bizarres en orange)

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
clear all;
clf;
% Initialisation of alpha and beta
alpha =1.0;
beta=0.1;
% Input the initial position
x1=input('Enter the initial value of the displacement in meters:');
% set up the initial velocity, mass and time step
Nstep = 601
v1=0; % Initial velocity in meters/s.
mass = 1.0 % Mass of particle
step = 0.02 %time step in seconds

for istep = 1: Nstep % make the number of steps necessary
t1 (istep) = (istep -1) * step; % time
Force1 = alpha * x1 - 4 * beta * x1^3; % compute the force on the particle
accel1 = Force1 / mass; % compute the accel. on the particle
% compute the velocity and position using Euler algo.
x1 = x1 + v1 * step; 
v1 = v1 + accel1 * step; 

if v1 == 0
x=x1;
disp(x)
t1 (istep) = t; % t n'est pas déclarée précédemment
disp(t)
end
end

Jstep=istep+0.02;
for jstep = Jstep:Nstep % Jstep=1.02 à la première itération : ne peut être un indice de boucle
t2 (jstep) = (jstep -1) * step; % time
Force2 = alpha * x2 - 4 * beta * x2^3; % compute the force on the particle
accel2 = Force2 / mass; % compute the accel. on the particle
x2 = x2 + v2 * step; 
v2 = v2 + accel2 * step;


if (x2==x)
Periode=t2 (jstep) - t;
disp (Periode)
end
end

[.Ophiuchus.]

Bonsoir, merci beaucoup pour votre aide Tubaas! J'ai pris compte de vos remarques et de vos conseils et j'ai modifié mon programme en conséquence voyez un peu. Le problème est qu'il ne fonctionne toujours pas... En fait j'arrive à évaluer un semblant de période mais comme la vitesse est comprise entre 0 et 1 et que je me base sur sa partie entière (toujours égale à 0), les résultats sont aléatoires...
J'ai bien penser à une autre méthode consistant à calculer les valeurs de ("istep") pour lesquelles la vitesse ("vitesse (istep)") change de signe mais pour cela il faudrait que je mette au point une gestion dynamique de mes variables vitesse pour pouvoir les comparer et je n'y parviens pas. De plus mon graphique est faux je ne sais plus trop quoi faire pour le coup :/

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
clear all;
 
% Initialisation of alpha and beta
alpha =-1.0;
beta = 0.0;
 
% set up the initial velocity, mass, position and time step
Nstep = 601;
 
% Initialisationde la variable x à 1 pour caller l'échelle des x avec les valeurs de la question 1 du TP
x = 1;
v = 0; % Initial velocity in meters/s.
mass = 1.0; % Mass of particle
step = 0.02; %time step in seconds
 
for istep = 1 : Nstep % make the number of steps necessary
 
	t = temps (istep) = (istep -1) * step; % time
	Force = alpha * x - 4 * beta * x^3; % compute the force on the particle
	accel = Force / mass; % compute the accel. on the particle
 
	% compute the velocity and position using Euler algo.
	x = position (istep) = x + v * step; 
	v = vitesse (istep) = v + accel * step; 
 
		% La vitesse étant une fonction dérivée de la position, alors si celle-ci s'annule nous obtenons la valeur maximum ou minimum de la position
		% en recherchant les deux extréma les plus proches, on obtient par la différence de temps correspondante la période d'oscillation du système.
		% Ici on prend la partie entière puisque la fonction que nous utilisons pour tracer la courbe représentative de la vitesse est une
		% approximation donc on n'a peu de chance de tomber sur une valeur de v égale à 0
		if (floor (vitesse (istep)) == 0) 
 
			I1 = istep; % On récupère la valeur de istep pour laquelle la vitesse s'annule
 
			% On l'intègre dans une autre variable qui va nous servir à initialiser le pas de recherche de l'autre valeur de X; la multiplication par deux 				
			% est nécessaire sinon les valeurs pour lesquelles v=0 sont trop proches les unes des autres dans le temps pour que les valeurs de X 				
			% correspondantes soient distantes d'une période
			global S3 = 2 * I1; 
 
			global X1 = position (I1); % On recherche la position correspondante à I1
			T1 = temps (I1); % Et également le temps correspondant
 
			% Affichage des éléments de réponse
			str = [" A la position X1 = ", num2str(X1), " la vitesse s'annule ce qui permet d'affirmer que la trajectoire du point admet un maximum ou un minimu à T1 = ", 				num2str(T1)] 
			break % On arrète les calculs
		endif
endfor	
 
	% On remplace 1 par S3 afin de ne pas refaire le même calcul que précédemment	
	for jstep = S3 : Nstep % make the number of steps necessary
 
		temps (jstep) = (jstep -1) * step; % time
                Force = alpha * x - 4 * beta * x^3; % compute the force on the particle
                accel = Force / mass; % compute the accel. on the particle
 
                % compute the velocity and position using Euler algo.
                position (jstep) = x = x + v * step; 
                vitesse (jstep) = v = v + accel * step; 			
 
	% Même condition que dans le 1er cas sauf que l'on ajoute également que la nouvelle valeur de x doit avoir le même signe que X1 dans le cas contraire on aurait 		% une demi-période seulement
	if (floor (vitesse (jstep)) == 0)
		I2 = jstep;
		X2 = position (I2);
		elseif (sign (X2) == sign (X1)) 
			T2 = temps (I2);
			str = [" A la position X2 = ", num2str(X2), " la vitesse s'annule ce qui permet d'affirmer que la trajectoire du point admet un maximum ou un minimum à T2 = ", 				num2str(T2)]
 
		% Calcul de la période totale
               Periode = T2 - T1;
               str = [" La période totale est de ", num2str(Periode), " secondes"]	
               break
	endif
endfor
 
warning off;
clf;
% Tracé des graphes en échelle semi-logarithmique
Ostep = floor (Periode) + 2;
for kstep = 1 : Ostep % make the number of steps necessary
	t = temps (kstep) = (kstep -1) * step; % time
        Force = alpha * x - 4 * beta * x^3; % compute the force on the particle
        accel = Force / mass; % compute the accel. on the particle
 
        % compute the velocity and position using Euler algo.
        position (kstep) = x = x + v * step; 
        vitesse (kstep) = v = v + accel * step;
 
	% Tracé de la position
	XlimitM =max(position) + 0.1;	
	Xlimitm =min(position) - 0.1;
 
	axis([0, Nstep * step, Xlimitm, XlimitM]); % Set axis limits
	semilogy (temps, position, 'b+', temps, vitesse, 'r-'); % plot the position and the velocity	
	xlabel('Time (s)'); % x-axis label
	ylabel('Position (m), Velocity (m.s¯¹)'); % y-axis label
	title('X=f(t) sur une période')
	legend( 'Position (m)' , 'Velocity (m.s¯¹)')
endfor