Résolution numérique

Accueil Remonter Introduction Equation du pendule Solution de l'équation Résolution numérique Conclusion Bibliographie

3. Résolution de l'équation générale et graphes commentes de l'amplitude et de l'énergie totale.

A l'aide d'un programme en tubo-pascal édité plus loin nous calculons en fonction du temps t l'angle q ainsi que dq /dt pour obtenir des diagrammes de phase. Avec toutes ces données, nous estimons l'énergie totale du système. Pour tracer les courbes nous utiliserons gnuplot.

L'équation introduite dans le programme est valable pour un nombre indéterminé de période et est:

elle est sans dimension afin que nous puissions comparer les différents coefficients.

L’équation de l’énergie est donnée par :

et s’écrit, en remplaçant t par t’ :

Etot=1/2 mlg(dq /dt’)²+mgl(1-cosq )

Nous poserons par besoin de simplification :

m= 1 kg

l= 1 m .

A l’intérieur du programme, l’équation est résolue de façon numérique à l’aide de la méthode de Runge-Kutta, on résoud simultanément deux équations du premier ordre :

dq /dt = p = f1(q ,pt)

q /dt² =dp/dt = f2(q ,pt)

A partir de la valeur initiale q 0, on calcul la solution des équations point par point.

Ici, on utilise la méthode d'ordre 4, l’erreur d’approximation est donc localement en h5, h étant le pas d’intégration qui sera d’autant plus faible que l’intervalle de temps du calcul sera faible.

Program pendule;

uses crt;

var z,p,x,t,k1,k2,k3,k4 :real;

l1,l2,l3,l4 :real;

a,b,h,fa,fb,c,e,m,r,amort,tau:real;

i:integer;

tab:array[0..500,1..3] of real;

(* cette procédure initialise le tableau à enregistrer dans le fichier'*)

procedure init;

begin

clrscr;

for i:= 0 to 500 do

begin

tab[i,1]:=0;

tab[i,2]:=0;

end;

end;

(*définition de la fonction annexe p = x' *)

function f1(t,x,p:real):real;

begin

f1:= p;

end;

(* expression de x" *)

function f2(t,x,p:real):real;

begin

(* tau:=sqrt(r/9.8); *)

(* fa:=amort*tau/(m*r); *)

if p>0

then

f2:= (-fb*p*p-amort*p-sin(x)-c)

else

f2:= (fb*p*p-amort*p-sin(x)+c);

end;

(* entrée des données *)

procedure entree;

begin

writeln('entrez la masse du pendule');

readln(m);

writeln('entrez la longueur du pendule');

readln(r);

writeln('entrez le coefficient quadratique ');

readln(fb);

writeln('entrez le coefficient d''amortissement');

readln(amort);

writeln('entrez les bornes de l''intervalle ');

readln(a,b);

writeln('entrez le coefficient de frottement solide ');

readln(c);

writeln(' le nb de points (limité à 500)');

readln(i);

h:=(b-a)/i;

writeln(' et les conditions initiales x(t0) et p(t0)');

readln(x,p);

end;

(* résolution du sytème couplé des 2 équations différentielles du 1er ordre *)

(* algorithme de Runge Kutta d'ordre 4 *)

procedure calcul;

begin

writeln(' t x p');

t:=a; i:=1;

while t<b-h do

begin

k1:= h*f1(t,x,p);

l1:= h*f2(t,x,p);

k2:=h*f1(t+(h/2),x+(k1/2),p+l1/2);

l2:=h*f2(t+(h/2),x+(k1/2),p+l1/2);

k3:= h*f1(t+(h/2),x+(k2/2),p+l2/2);

l3:= h*f2(t+(h/2),x+(k2/2),p+l2/2);

t:=t+h;

k4:= h*f1(t,x+k3,p+l3);

l4:= h*f2(t,x+k3,p+l3);

x:= x + (k1+2*k2+2*k3+k4)/6;

p:= p + (l1+2*l2+2*l3+l4)/6;

e:=(0.5*m*r*p*p*9.81)+(9.81*m*r*(1-cos(x)));

tab[i,1]:=t;

tab[i,2]:=x;

tab[i,3]:=e;

writeln (t:8:4,x:16:8,p:16:8,e:16:8);

i:=i+1;

end;

end;

procedure crefich1;

(* ce programme crée un fichier de nombres réels générés par un tableau*)

(* les nombres sont stockés en mémoire centrale dans le tableau tab *)

(* et sur le disque dur avec un nom à préciser *)

var

f1 : string[12];

f:text;

i,j :integer;

begin

writeln ('entrez le nom du fichier à créer');

readln (f1);

(* ouverture de la zone tampon sous le nom donné par l'utilisateur *)

assign(f,f1);

(* reinitialisation du fichier et positionnement du pointeur sur l'élément 0 *)

rewrite (f);

(* écriture des éléments dans le fichier*)

for i:= 0 to 500 do

begin

write (f,tab[i,1],' ',tab[i,2]);

writeln(f);

end;

close (f);

(* fermeture du fichier *)

end;

procedure energie;

(*ce programme crée un fichier pour les valeurs de l'énergie*)

var

f3: string[12];

g:text;

i,j :integer;

begin

writeln('entrez le nom du fichier énergie');

readln (f3);

assign(g,f3);

rewrite (g);

for i:=0 to 500 do

begin

write (g,tab[i,1],' ',tab[i,3]);

writeln(g);

end;

close (g);

end;

 

 

begin

init;

entree;

calcul;

crefich1;

energie;

end.

Un programme plus simple permet d'obtenir le diagramme de phase.

Sauf mention spéciale le pendule théorique utilise pour tracer les courbes a une masse de 1 Kg et une longueur de 1 m.

Nom

a

b

c

Duré (s)

angle initial(rad)

amplia/energa

0,1

0

0

50

2,007

amplib/energb

0

0,1

0

50

2,007

amplic/energc

0

0

0,1

50

2,007

amc0/enerc0

0,1

0,1

0

50

2,007

amb0/enerb0

0,1

0

0,1

50

2,007

ama0/enera0

0

0,1

0,1

50

2,007

amtot/enertot

0,1

0,1

0,1

50

2,007

peta/petea

0,1

0

0

160

0,5

petb/peteb

0

0,1

0

160

0,5

petc/petec

0

0

0,1

160

0,5

phas'a

0,1

0

0

50

2,007

phas'c

0

0

0,1

50

2,007

hary1/haro1

0,1

0

0,1

30

2,007

hary2/haro2

0,1

0

0,01

30

2,007

hary3/haro3

0,01

0

0,1

30

2,007

a1/ae1/pha1

0

0

0,014

50

0,075

a2/ae2/pha2

0,011

0

0,014

50

0,075

a3/ae3/pha3

0,38

0

0,014

50

2,009

a4/ae4/pha4

0,34

0

0,013

50

2,009

a5/ae5/pha5

0,34

0

0,013

50

0,94

En observant les courbes représentant l’énergie totale on a l’impression que la pente de l’énergie s’annule pour certaines valeurs de ,nous allons donc tenter de le vérifier.

On a :

Em=1/2mv²+mgl(1-cos)

=1/2ml²(d/dt)²+mgl(1-cos)

 

On cherche donc à résoudre :

dEm/dt=ml²(/dt².d/dt)+mgl(d/dt )sin =0

On obtient deux solutions :

-soit d/dt0 ,on a d²/dt²=-(g/l)sin,

On retrouve l’équation de l’oscillateur non amorti.

-soit d/dt=0 ,

Pour le cas qui nous intéresse, nous ne retiendrons que la deuxième solution qui nous indique que la variation d’énergie mécanique est nulle lorsque la vitesse du pendule s’annule.

Le graphe ‘ampli’, ’energ’ illustre bien ce qui vient d’être dit, en effet on constate que lorsque la vitesse s’annule c’est à dire lorsque l’amplitude est maximum, la pente de l’énergie est nulle.

 

Commentaires sur le graphe amplia4, amplib4, amplic4.

Ce graphe permet de bien visualiser la forme de la décroissance que fait subir chacun des coefficients de frottements à l’amplitude. Cette décroissance est linéaire pour le coefficient solide alors qu’elle est exponentielle pour les coefficients quadratique et linéaire.

Commentaires sur les graphes 'amplia,b,c' et 'energa,b,c' et 'phas'a,c'.

Ces schémas ont été tracés pour une seule valeur non nulle des coefficients A, B, C (0.1) de l'équation (I,3). Dans ce cas nous lâchons le pendule d'un angle de 2 rad.

Nous constatons que pour les toutes premières oscillations l'amortissement est plus important pour A¹ 0 puis B et enfin C qui parait le moins efficace. Mais passé la troisième demi-période le phénomène s’inverse. Ce paradoxe s'explique bien par le fait que la force exercée par coefficient C est indépendant de la vitesse et va donc agir de manière linéaire. En revanche les forces exercées par A et B sont fonction respectivement de la vitesse au carre et de la vitesse c'est pourquoi elles agissent d'autant plus fort que la célérité est grande. Lorsqu'elle s'affaiblit (à partir de 15s) la force issue de C écrase de manière constante l'amplitude et la conduit bien plus vite que les autres vers 0 (à t=25s). La force issue de A décrois plus vite que toute autre puisqu'elle est quadratique. Le

diagramme de phase montre que tous les systèmes tendent vers un état stable mais qu'au bout de quinze secondes seul celui avec C non nul l'atteint.

Les calculs ayant été fait avec une équation dédimensionnée et des valeurs identiques pour les coefficients non nuls on peut donc conclure que la force de frottement solide devient la plus forte à mesure que la vitesse décrois.

 

On constate aussi que la dérivée de l'énergie est toujours négative ou nulle ce qui confirme que les frottements prennent de l'énergie au système. La pente horizontale correspond donc à une vitesse nulle ce qui revient à dire que lorsque le pendule est immobile il n’interagit pas avec le fluide.

Par l'observation du graphe de l'énergie on constate que la partie "horizontale" des courbes est de largeur différente pour chacun des coefficients de frottements. On peut donc penser que le temps passe par le pendule a très faible vitesse va dépendre des différents coefficients de frottement. A priori ce temps est plus faible pour le coefficient solide, et plus élevé pour le coefficient quadratique.

 

Commentaires sur les graphes 'amc0,a0,b0,tot'.

On a pris pour chaques graphe un coefficient nul et les autres égaux à 0.1(sauf pour les graphes amtot et enertot où tout les coefficients étaient égaux à 0.1).L'énergie tout comme l'amplitude diminue plus rapidement lorsque tous les coefficients de frottements sont présents.

En comparent 'amplic' (a=0,b=0) et 'amc0' (c=0) on constate que l'amplitude arrive à 0 plus rapidement avec le coefficient de frottement solide (c=0.1) qu'avec deux coefficients de frottement visqueux réunis (0.1 chacun). Ceci est sans doutes dû au fait que notre coefficient solide est démesurément grand.

Graphes 'peta,b,c'.

Ces graphes représentent les oscillations du pendule aux petites amplitudes ( on lâche le pendule a 0.5 rad ) en changeant a chaque fois la nature du frottement de coefficient 0.1 .

Les graphes 'hary1,2,3' et 'haro1,2,3' .

Le graphe 1 tend vite vers 0, le graphe 2 va en revanche osciller plus longtemps, c'est donc le fait que le coefficient c est dix fois moins fort que dans la première courbe et dix fois moins que a. Le fait que l'on ait dédimensionné l'équation mous permet de dire qu'un facteur dix entre deux coefficients va fortement influencer le pendule. Ainsi si nous avions un rapport de cent entre a et c il serrait possible de négliger c.

Le graphe a1 est un bon exemple de la décroissance linéaire due au seul coefficient c.

Le graphe pha3 est particulièrement intéressant car il montre que la vitesse maximum est atteinte avant q =0. Ce paradoxe apparent s'explique par le fait que le coefficient a=0.38 est relativement élevé par rapport aux autres. Comme il agit avec le carre de la vitesse plus celle ci augmente plus la force fluide freine. Dans le vide la vitesse maximum serrait atteinte à q =0 mais dans un fluide a cet endroit les forces de dissipations ont déjà commencé à ralentir le pendule ce qui explique cette fausse bizarrerie.

 Précédente ] Accueil ] Remonter ] Suivante ]

Page mise à jour le 25/11/06