|
Menu |
APPWiki /
H08Cours01ProjetSessionProjet 1 : Résolution numérique d'un système complexeIndex
IntroductionCe projet d'envergure vise à favoriser l'intégration de plusieurs concepts de physiques tant sur le plan théorique, numérique qu'expérimental. Il se déroulera durant les 11 premières semaines de la session. Le thème du projet sera choisi parmi la liste ci-dessous. Chaque étudiant devra répartir sur ce formulaire?, 100 unités qu'il attribuera au(x) projet(s) de son choix pour établir son intérêt face à un thème ou l'autre. Le professeur s'assurera ensuite de former les équipes en veillant à avoir des groupes optimisant l'intérêt face au thème. Cette démarche déterminera les équipes pour tous les autres projets du cours. Les équipes ne devront pas comporter plus de 4 membres. Liste des thèmes et hyperliens pertinents
Les systèmes complexes sont souvent impossibles à résoudre analytiquement faute d'outils mathématiques suffisamment avancés. Pour contourner ce problème, un certain nombre de méthodes de résolution numérique ont été élaborées. Ces méthodes ont démontré tout leur intérêt avec l'avènement d'ordinateurs puissants accessible à tous. Méthodes numériquesLa méthode numérique se résume à effectuer l'intégrale temporelle d'une fonction. Dans le cas qui nous intéresse la fonction à intégrer sera l'accélération. La première étape de votre projet est donc d'élaborer une fonction donnant l'accélération par rapport aux variables comme la vitesse (p.ex dans le cas de frottement fluide), la position (p.ex. dans le cas de force de gravitation), ou tout autre paramètre physique pertinent. Il existe plusieurs façons de réaliser cette intégration numérique. La plus simple est la méthode de Euler qui a de désavantage de générer une erreur plus grande. La méthode Runge-Kutta est par contre très précise mais a le désavantage non négligeable d'être assez complexe à mettre en place. Aussi dans le cadre de votre projet je vous suggère fortement d'opter pour une méthode de précision intermédiaire nommée méthode de Verlet. Si vous désirez obtenir de plus amples informations sur les méthodes numériques je vous suggère de lire ce document synthèse que j'ai écris document. Méthode de VerletDans le texte qui suit, un variable suivie de l'apostrophe indique la dérivée par rapport au temps de cette variable. S'il y a deux apostrophes alors il s'agit de la dérivée seconde par rapport au temps. Par exemple si la variable est x, x' represente la vitesse selon x et x'' l'accélération selon x. L'indice sous la variable indique quant à lui le numéro du pas de calcul ( {# \Delta t #} ). Par exemple x0 indique la position à t=0 et xi, la position à {# t = i \Delta t #}. La méthode de Verlet s'exprime comme suit: {## x_{i+1} = x_{i} + x'_{i} \Delta t + \frac{1}2 x''_{i} {\Delta t}^2 ##} {## x'_{i+1} = x'_{i} + \frac{1}2 (x_{i+1} + x_{i}) {\Delta t} ##} Pour être implémentée cette méthode requiert que l'accélération à t=0 et {# t = \Delta t #}. Ces valeurs devront être estimées par vous au meilleur de vos connaissances. Il y aurait aussi moyen d'utiliser la méthode de Euler pour trouver l'accélération à {# t = \Delta t #}. La détermination de l'accélération à t=0 ne pose pas de problème puisque vous devrez poser les valeurs initiales de position et vitesse. Votre premier défi consiste à trouver une expression de votre accélération en x et en y. Pour y arriver vous exprimez la force sur la première masse m1 de la façon suivante: {## a_{1x} = \frac{\sum F_{1x}}{m_1} ##} {## a_{1y} = \frac{\sum{F_{1y}}}{m_1} ##} Vous avez donc à déterminer l'équation donnant la somme de force exercées sur la masse 1. Procédez de la même façon pour chaque masse impliquée dans votre problème. Application: Orbite d'un satelliteAppliquons cette méthode à un problème simple qui consiste à résoudre le mouvement d'un satellite autour d'un très grande masse telle que le Soleil. Pour simplifier nous prendrons pour acquis que le mouvement du Soleil ne sera pas perturbé par le satellite (dans vos projets vous ne pourrez pas faire cette supposition). Notre exemple se réduit donc à trouver l'accélération en x et y d'une seule masse soit cette du satellite. L'unique force agissant sur le satellite est la force de gravité provoquée par la présence du Soleil. Cette force est donnée par l'équation suivante: {## F_x = \frac{G {m_s} m (x-x_s)}{r^3} ##} Ici x est la position du satellite en x, xs est la position du Soleil en x et r est la distance séparant le satellite et le Soleil. Une équation analogue peut être écrite en y. r peut être calculé avec le théorème de Pytagore: {## r = \sqrt{(x-x_s)^2+(y-y_s)^2} ##} Nous connaissons donc maintenant l'accélération en x et en y car elles sont données par la force divisée par la masse du satellite m. L'étape suivante consiste à insérer toutes ces équations dans un tableau excel. L'orbite du satellite dans excelLe tableau ci-dessous représente ce que vous devrez entrer dans excel pour résoudre ce problème.
Si vous entrez ces données dans excel ou openoffice, vous devriez avoir un résultats semblable à ceci:  Sur cette figure, les cases jaunes sont des données fixes (des constantes de la nature ou des conditions initiales) alors que les cases orangées sont de simple copies (Ctrl-c Ctrl-v) de la case qui précède verticalement. Notez que le graphique de x en fonction de y (c.-à-d. la forme de la trajectoire) donne un cercle. En fait dans ce problème j'ai simulé un objet de 1 kg en orbite autour du Soleil à la distance Terre-Soleil qui aurait une période orbitale de 365,25 jours. Dans cet exemple nous avons utilisé un pas de calcul de 100000 sec (à peu pres 1 jour) mais qu'est-ce qui se passe si nous définissons un pas de calcul beaucoup plus long. Pour répondre à cette question, j'ai fixé le pas de calcul à 1500000 sec, soit près de deux semaines. Voici le résultat:  Vous remarquez clairement que la trajectoire n'est plus du tout circulaire ce qui démontre qu'un choix de pas de calcul trop grand génère des erreurs numériques graves. Il faudra donc porter une attention particulière à ce choix pour votre projet. Dans notre exemple, j'ai calculé la vitesse initiale de la masse de sorte qu'elle puisse faire un tour en un an sur une orbite circulaire. Voyons ce qui se passe si nous augmentons de 20% la vitesse initiale (nous revenons aussi au pas de calcul de 100000sec).  Notez que les limites des axes a été augmentée afin que pouvoir toute la trajectoire. Il apparait que si la vitesse n'es pas bien choisie, l'orbite n'est plus circulaire mais plutot elliptique tel que la loi de Kepler l'énonçait. Modèle de page excel pour le problème à 3 corps soumis à la gravitationProgrammes fortran pour le mouvement des astresPour démontrer certains effets plus subtils, il faut absolument être en mesure de simuler à haute résolution temporelle un très grand nombre d'orbites. Excel ne permet pas de réaliser cette expérience, c'est pourquoi j'ai écris ce petit programme fortran. Ce programme a été conçu pour expliquer la formation de la division de Cassini dans les anneaux de Saturne mais il pourrait être modifié pour tout autre problème. Compilateur fortran sous windows
Programme fortran pour le problème à 3 corpsCe code produit des fichiers de position en fonction du temps pour chacun des corps. c c programme pour simuler la division de cassini c programme par Martin Aube a partir des equations developpees par c l equipe de Jonathan Bechette c utilisation de la methode de verlet c Ce programme peut servir à simuler tout autre problème à 3 corps c dont l orbite de chacun est centre sur le centre de masse du systeme c Le programme pourrait etre transforme pour simuler un systeme multiple c sous reserve de modifier la section de determination des vitesses c orbitales initiales c c ===== c Declaration des variables c double precision xoa,yoa,xon,yon,xsa,ysa,xsn double precision ysn,xma,yma,xmn,ymn,rmi double precision rsi,roi,rso,rmo,rms,dt,dti,G double precision ms,mm,mo,Tm,To,Ts,noo,nom,pi double precision vxoa,vyoa,vxon,vyon,vxsa,vysa double precision vxsn,vysn,vxma,vyma double precision vxmn,vymn,axoa,ayoa,axon,ayon double precision axsa,aysa,axsn,aysn,axma double precision ayma,axmn,aymn integer step,nloop,i,mout,n,ne c c ===== c Identification des variables c c ms=masse de Saturne c mo=masse d une particule de l anneau c mm=masse de Mimas c rmi=distance de Mimas par rapport au centre de masse (cm) c roi=distances de la particule pr rapport au cm c rsi=distances de Saturne pr rapport au cm c rms=distance entre la particule et Saturne c rmo=distance entre Mimas et la particule c rso=distance entre Saturne et la particule c xoa et yoa sont les position de la particile a t_i c xon et yon sont les position de la particile a t_n c xma et yma sont les position de Mimas a t_i c xmn et ymn sont les position de Mimas a t_n c xsa et ysa sont les position de Saturne a t_i c xsn et ysn sont les position de Saturne a t_n c vxoa et vyoa sont les vitesses de la particile a t_i c vxon et vyon sont les vitesses de la particile a t_n c vxma et vyma sont les vitesses de Mimas a t_i c vxmn et vymn sont les vitesses de Mimas a t_n c vxsa et vysa sont les vitesses de Saturne a t_i c vxsn et vysn sont les vitesses de Saturne a t_n c axoa et ayoa sont les accelerations de la particile a t_i c axon et ayon sont les accelerations de la particile a t_n c axma et ayma sont les accelerations de Mimas a t_i c axmn et aymn sont les accelerations de Mimas a t_n c axsa et aysa sont les accelerations de Saturne a t_i c axsn et aysn sont les accelerations de Saturne a t_n c dt=pas de calcul c dti=pas de calcul initial c Tm,To,Ts periodes orbitales c nloop=nombre de boucles de calcul c mout=saut pour l impression de sorties c c ===== c Sorties du modele c c mimas.dat -> position x et y de Mimas c rmimas.dat -> rayon orbital de Mimas en fonction du temps c objet.dat -> position x et y de la particule c robjet.dat -> rayon orbital de la particules en fonction du temps c saturne.dat -> position x et y de Saturne c rsaturne.dat -> rayon orbital de Saturne en fonction du temps c c ===== c Valeurs initiales c c roi=1.16857E+8 c dti=0.01 dt=0.01 n=0 ne=0 pi=3.14159 G=6.6725985E-11 ms=5.6851E+026 mo=1. mm=3.7493E+022 rmi=1.855E+8 roi=1.16857E+8 rsi=(rmi*mm+roi*mo)/ms rms=rsi+rmi rmo=dabs(rmi-roi) rso=rsi+roi xoa=roi yoa=0. xsa=-rsi ysa=0. xma=rmi yma=0. c c ===== c Calcul des vitesses orbitales initiales c print*,'Calcul des vitesses initiales...' vxoa=0. vyoa=(dabs(roi*(G*ms*(xsa-xoa)/rso**3.+G*mm*(xma-xoa)/rmo**3.))) +**0.5 vxma=0. vyma=(dabs(rmi*(G*ms*(xsa-xma)/rms**3.+G*mo*(xoa-xma)/rmo**3.) +))**0.5 vxsa=0. vysa=-(dabs(rsi*(G*mm*(xma-xsa)/rms**3.+G*mo*(xoa-xsa)/rso** +3.)))**0.5 c c ===== c Estimation des périodes orbitales c print*,'Calcul des périodes orbitales...' Tm=2.*pi*rmi/vyma To=2.*pi*roi/vyoa Ts=2.*pi*rsi/vysa c c ===== c c nloop=2000000000 mout=int(1000./dt) if (mout.gt.100000) mout=100000 print*,"Periode de l objet=",To/3600./24.,"Periode de Mimas=", +Tm/3600./24.," jours" noo=dble(nloop)*dt/To nom=dble(nloop)*dt/Tm print*,"Orbites de l objet=",noo,"Orbites de Mimas=",nom open(unit=1,file="objet.dat",status="unknown") open(unit=2,file="mimas.dat",status="unknown") open(unit=3,file="saturne.dat",status="unknown") open(unit=4,file="robjet.dat",status="unknown") open(unit=5,file="rmimas.dat",status="unknown") open(unit=10,file="rsaturne.dat",status="unknown") print*,'a' do i=1,nloop rso=((xoa-xsa)**2.+(yoa-ysa)**2.)**0.5 rmo=((xma-xoa)**2.+(yma-yoa)**2.)**0.5 rms=((xma-xsa)**2.+(yma-ysa)**2.)**0.5 axon=G*(ms*(xsa-xoa)/rso**3.+mm*(xma-xoa)/rmo**3.) ayon=G*(ms*(ysa-yoa)/rso**3.+mm*(yma-yoa)/rmo**3.) axmn=G*(ms*(xsa-xma)/rms**3.)+G*(mo*(xoa-xma)/rmo**3.) aymn=G*(ms*(ysa-yma)/rms**3.)+G*(mo*(yoa-yma)/rmo**3.) axsn=G*(mm*(xma-xsa)/rms**3.)+G*(mo*(xoa-xsa)/rso**3.) aysn=G*(mm*(yma-ysa)/rms**3.)+G*(mo*(yoa-ysa)/rso**3.) if (i.eq.1) then c c ===== c Ajuste la premiere valeur prededente a la valeur actuelle c dt=dti/100. axoa=axon ayoa=ayon axsa=axsn aysa=aysn axma=axmn ayma=aymn else dt=dti endif vxon=vxoa+0.5*(axoa+axon)*dt vyon=vyoa+0.5*(ayoa+ayon)*dt vxsn=vxsa+0.5*(axsa+axsn)*dt vysn=vysa+0.5*(aysa+aysn)*dt vxmn=vxma+0.5*(axma+axmn)*dt vymn=vyma+0.5*(ayma+aymn)*dt xon=xoa+vxoa*dt+0.5*axoa*dt**2. yon=yoa+vyoa*dt+0.5*ayoa*dt**2. xsn=xsa+vxsa*dt+0.5*axsa*dt**2. ysn=ysa+vysa*dt+0.5*aysa*dt**2. xmn=xma+vxma*dt+0.5*axma*dt**2. ymn=yma+vyma*dt+0.5*ayma*dt**2. roi=(xon**2.+yon**2.)**0.5 rmi=(xmn**2.+ymn**2.)**0.5 rsi=(xsn**2.+ysn**2.)**0.5 n=n+1 c c ===== c Impression des resultats c if (n.eq.mout) then n=0 ne=ne+1 print*,ne,"/",nloop/mout, + " temps(h)=",(dble(i-1)*dti+dti/100.)/3600. write(4,*),dble(i-1)*dti+dti/100.,roi write(5,*),dble(i)*dti+dti/100.,rmi write(7,*),dble(i)*dti+dti/100.,rsi write(1,*) xon,yon write(2,*) xmn,ymn write(3,*) xsn,ysn endif c c ===== c store la valeur precedentes c axoa=axon ayoa=ayon axsa=axsn aysa=aysn axma=axmn ayma=aymn vxoa=vxon vyoa=vyon vxsa=vxsn vysa=vysn vxma=vxmn vyma=vymn xoa=xon yoa=yon xsa=xsn ysa=ysn xma=xmn yma=ymn enddo close(unit=1) close(unit=2) close(unit=3) close(unit=4) close(unit=5) close(unit=10) stop end Programme multi objetCe code produit la position et vitesse finale du système (fichier cf.txt). Tous les calculs sont effectués dans le plan xy. c c programme pour simuler le probleme du mouvement à N corps sous c l action de la force de gravite c programme par Martin Aube c utilisation de la methode de verlet c c ===== c Declaration des variables c double precision xa(1001),xn(1001),ya(1001),yn(1001) double precision dt,dti,G double precision m(1001),pi,ray double precision vxa(1001),vya(1001),vxn(1001),vyn(1001) double precision axa(1001),aya(1001),axn(1001),ayn(1001) integer step,nloop,i,j,mout,n,ne,nm,nn,ii character*30 nomvect c c ===== c Identification des variables c c m(n)=masse de la masse n c nm=nombre de masses c xa(n) ya(n) les position de la particile n a t_i c xn(n) et yn(n) sont les position de la particile n a t_n c vxa(n) et vya(n) sont les vitesses de la particile n a t_i c vxn(n) et vyn(n) sont les vitesses de la particile n a t_n c axa(n) et aya(n) sont les accelerations de la particile n a t_i c axn(n) et ayn(n) sont les accelerations de la particile n a t_n c ray=rayon entre deux masses c dt=pas de calcul c dti=pas de calcul initial c Tm,To,Ts periodes orbitales c nloop=nombre de boucles de calcul c mout=saut pour l impression de sorties c c c Commentaires: c Le programme attend en entree un fichier nomme ci.txt (conditions initiales) c et retourne les conditions finales cf.txt dans le meme format c le format est le nombre de masses sur la premiere ligne suivie sur les lignes c suivantes (une ligne par masse) des x y vx vy m c c Le programme fait 2E9 boucles de calculs temporels. La variable dti (a ajuster c avant de compiler) est le pas de calcul. Pour connaitre le temps apres les c 2E9 boucles il faut faire tf = 2E9 x dti - 0.99dti c Dans sa version originale dti=1jour ce qui implique un tf>5E6 annees c c ===== c Valeurs initiales c c dti=86400*7 => une semaine c dti=86400.*0.5 nloop=100000 mout=10000 n=0 pi=3.14159 G=6.6725985E-11 dt=dti c c ===== c Lecture des positions et vitesses initiales c open(unit=1,file='ci.txt',status='unknown') read(1,*) nm do i=1,nm read(1,*) xa(i),ya(i),vxa(i),vya(i),m(i) enddo close(unit=1) c c ===== c c c ouverture du fichier de sortie open(unit=2,file="cf.txt",status="unknown") do n=1,nloop do i=1,nm do j=1,nm if (i.ne.j) then ray=((xa(j)-xa(i))**2.+(ya(j)-ya(i))**2.)**0.5 axn(i)=axn(i)+G*m(j)*(xa(j)-xa(i))/ray**3. ayn(i)=ayn(i)+G*m(j)*(ya(j)-ya(i))/ray**3. endif enddo if (n.eq.1) then c c ===== c Ajuste la premiere valeur prededente a la valeur actuelle c dt=dti/100. axa(i)=axn(i) aya(i)=ayn(i) else dt=dti endif vxn(i)=vxa(i)+0.5*(axa(i)+axn(i))*dt vyn(i)=vya(i)+0.5*(aya(i)+ayn(i))*dt xn(i)=xa(i)+vxa(i)*dt+0.5*axa(i)*dt**2. yn(i)=ya(i)+vya(i)*dt+0.5*aya(i)*dt**2. nn=nn+1 c c ===== c Impression des resultats c if (nn.eq.mout) then nn=0 print*," temps(jours)=",(dble(n-1)*dti+dti/100.)/3600./24. + ,'/',dble(nloop)*dti/3600./24., + nint(100.*(dble(n-1)+1/100.)/dble(nloop)),'%' endif nomvect='toto' if (nloop/10.eq.n) nomvect='vecteur01.txt' if (2*nloop/10.eq.n) nomvect='vecteur02.txt' if (3*nloop/10.eq.n) nomvect='vecteur03.txt' if (4*nloop/10.eq.n) nomvect='vecteur04.txt' if (5*nloop/10.eq.n) nomvect='vecteur05.txt' if (6*nloop/10.eq.n) nomvect='vecteur06.txt' if (7*nloop/10.eq.n) nomvect='vecteur07.txt' if (8*nloop/10.eq.n) nomvect='vecteur08.txt' if (9*nloop/10.eq.n) nomvect='vecteur09.txt' if (nloop.eq.n) nomvect='vecteur10.txt' c c ===== c Ecrire les donnees en format vector pour gnuplot c x,y,deltax,deltay if (nomvect.ne.'toto') then open(unit=3,file=nomvect,status="unknown") write(3,*) nm do ii=1,nm write(3,*) xa(ii),ya(ii),vxa(ii)*1.e9,vya(ii)*1.e9 enddo close(unit=3) endif c c ===== c store la valeur precedentes c axa(i)=axn(i) aya(i)=ayn(i) vxa(i)=vxn(i) vya(i)=vyn(i) xa(i)=xn(i) ya(i)=yn(i) enddo enddo c c ===== c Ecrire les donnees de sortie c write(2,*) nm do i=1,nm write(2,*) xa(i),ya(i),vxa(i),vya(i),m(i) enddo close(unit=2) stop end Pistes pour la rédaction du rapportUn rapport complet devra être remis. Il devra comprendre:
Notes importantes:
|