Problème N corps RK4 Fortran
Résolu
mitsuidewi
-
mys29000 Messages postés 1 Statut Membre -
mys29000 Messages postés 1 Statut Membre -
Bonjour à tous, j'ai écris un programme fortran qui me permet d'obtenir la trajectoire de N corps en interaction gravitationnelle, j'ai choisi le soleil comme centre du repère. mais en affichant les coordonnées de la terre je me suis aperçu que c'était faux, et je ne trouve pas l'erreur, quelqu'un pourrait-il m'aider?? je met mon code plus bas ainsi que les valeurs que j'ai utilisé.
[QUOTE]
PROGRAM Gravitation
IMPLICIT NONE
REAL (KIND=8), DIMENSION(:), ALLOCATABLE :: Fx, Fy, X, Y, Vx, Vy, m
REAL (KIND=8) :: t, tfin, h
REAL (KIND=8), PARAMETER :: G=6.67259d-11, A=1.49597870691d11
INTEGER (KIND=8) :: i, N
READ *, N, t, tfin, h
ALLOCATE(X(N), Y(N), Fx(N), Fy(N), Vx(N), Vy(N), m(N))
DO i=1,N
READ *, X(i), Y(i), Vx(i), Vy(i)
END DO
DO i=1,N
READ *, m(i)
END DO
DO WHILE (t<=tfin)
CALL RK4(X, Y, Vx, Vy, m, N, t, h)
t=t+h
END DO
END PROGRAM Gravitation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE Force(Fx, Fy, X, Y, m, N, t)
IMPLICIT NONE
REAL (KIND=8), DIMENSION(1:N) :: Fx, Fy, X, Y, m, Vx, Vy
REAL (KIND=8), PARAMETER :: G=6.67259d-11
INTEGER (KIND=8) :: i, j, N
REAL (KIND=8) :: t
DO i=1,N
DO j=1,N
IF (j /= i) THEN
Fx(i) = Fx(i) - (1/G)*m(j)*m(i)*(X(i)-X(j))/(((X(j)-X(i))**2+(Y(j)-Y(i))**2)**(1.5))
Fy(i) = Fy(i) - (1/G)*m(j)*m(i)*(Y(i)-Y(j))/(((X(j)-X(i))**2+(Y(j)-Y(i))**2)**(1.5))
END IF
END DO
END DO
END SUBROUTINE Force
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE RK4(X, Y, Vx, Vy, m, N, t, h)
IMPLICIT NONE
REAL (KIND=8), DIMENSION(1:N) :: Fx, Fy, X, Y, Vx, Vy, m
REAL (KIND=8), DIMENSION(1:N) :: kx1, kx2, kx3, kx4
REAL (KIND=8), DIMENSION(1:N) :: ky1, ky2, ky3, ky4
REAL (KIND=8), PARAMETER :: G=6.67259d-11
REAL (KIND=8) :: h, t
INTEGER (KIND=8) :: p, N
kx1=0; kx2=0; kx3=0; kx4=0; ky1=0; ky2=0; ky3=0; ky4=0
Fx=0; Fy=0
CALL Force(Fx, Fy, X, Y, m, N, t)
kx1=h*Fx/m
ky1=h*Fy/m
Fx=0; Fy=0
CALL Force(Fx, Fy, X+h*0.5*Vx+0.25*h*kx1, Y+h*0.5*Vy+0.25*h*ky1, m, N, t)
kx2=h*Fx/m
ky2=h*Fy/m
Fx=0; Fy=0
CALL Force(Fx, Fy, X+h*0.5*Vx+0.25*h*kx2, Y+h*0.5*Vy+0.25*h*ky2, m, N, t)
kx3=h*Fx/m
ky3=h*Fy/m
Fx=0; Fy=0
CALL Force(Fx, Fy, X+h*Vx+0.5*h*kx3, Y+h*Vy+0.5*h*ky3, m, N, t)
kx4=h*Fx/m
ky4=h*Fy/m
X=X+h*Vx+(h/6.0d0)*(kx1+kx2+kx3)
Y=Y+h*Vy+(h/6.0d0)*(ky1+ky2+ky3)
Vx=Vx+(kx1+kx2+kx2+kx3+kx3+kx4)/6.0d0
Vy=Vy+(ky1+ky2+ky2+ky3+ky3+ky4)/6.0d0
X(N)=0; Y(N)=0; Vx(N)=0; Vy(N)=0
PRINT *, t, X(2), Y(2)
END SUBROUTINE RK4
[QUOTE]
[QUOTE]
4 0 500 0.01
-1.01723985460 -0.206180767821 0.005269822149 -0.014839159853
-0.98323814725 -0.176124724912 0.002751137188 -0.016991888551
-0.98545358338 -0.176940858196 0.002966627399 -0.017595568343
0 0 0 0
1.8
3.98600440d14
4.902798d12
1.32712440018d20
[QUOTE]
[QUOTE]
PROGRAM Gravitation
IMPLICIT NONE
REAL (KIND=8), DIMENSION(:), ALLOCATABLE :: Fx, Fy, X, Y, Vx, Vy, m
REAL (KIND=8) :: t, tfin, h
REAL (KIND=8), PARAMETER :: G=6.67259d-11, A=1.49597870691d11
INTEGER (KIND=8) :: i, N
READ *, N, t, tfin, h
ALLOCATE(X(N), Y(N), Fx(N), Fy(N), Vx(N), Vy(N), m(N))
DO i=1,N
READ *, X(i), Y(i), Vx(i), Vy(i)
END DO
DO i=1,N
READ *, m(i)
END DO
DO WHILE (t<=tfin)
CALL RK4(X, Y, Vx, Vy, m, N, t, h)
t=t+h
END DO
END PROGRAM Gravitation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE Force(Fx, Fy, X, Y, m, N, t)
IMPLICIT NONE
REAL (KIND=8), DIMENSION(1:N) :: Fx, Fy, X, Y, m, Vx, Vy
REAL (KIND=8), PARAMETER :: G=6.67259d-11
INTEGER (KIND=8) :: i, j, N
REAL (KIND=8) :: t
DO i=1,N
DO j=1,N
IF (j /= i) THEN
Fx(i) = Fx(i) - (1/G)*m(j)*m(i)*(X(i)-X(j))/(((X(j)-X(i))**2+(Y(j)-Y(i))**2)**(1.5))
Fy(i) = Fy(i) - (1/G)*m(j)*m(i)*(Y(i)-Y(j))/(((X(j)-X(i))**2+(Y(j)-Y(i))**2)**(1.5))
END IF
END DO
END DO
END SUBROUTINE Force
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE RK4(X, Y, Vx, Vy, m, N, t, h)
IMPLICIT NONE
REAL (KIND=8), DIMENSION(1:N) :: Fx, Fy, X, Y, Vx, Vy, m
REAL (KIND=8), DIMENSION(1:N) :: kx1, kx2, kx3, kx4
REAL (KIND=8), DIMENSION(1:N) :: ky1, ky2, ky3, ky4
REAL (KIND=8), PARAMETER :: G=6.67259d-11
REAL (KIND=8) :: h, t
INTEGER (KIND=8) :: p, N
kx1=0; kx2=0; kx3=0; kx4=0; ky1=0; ky2=0; ky3=0; ky4=0
Fx=0; Fy=0
CALL Force(Fx, Fy, X, Y, m, N, t)
kx1=h*Fx/m
ky1=h*Fy/m
Fx=0; Fy=0
CALL Force(Fx, Fy, X+h*0.5*Vx+0.25*h*kx1, Y+h*0.5*Vy+0.25*h*ky1, m, N, t)
kx2=h*Fx/m
ky2=h*Fy/m
Fx=0; Fy=0
CALL Force(Fx, Fy, X+h*0.5*Vx+0.25*h*kx2, Y+h*0.5*Vy+0.25*h*ky2, m, N, t)
kx3=h*Fx/m
ky3=h*Fy/m
Fx=0; Fy=0
CALL Force(Fx, Fy, X+h*Vx+0.5*h*kx3, Y+h*Vy+0.5*h*ky3, m, N, t)
kx4=h*Fx/m
ky4=h*Fy/m
X=X+h*Vx+(h/6.0d0)*(kx1+kx2+kx3)
Y=Y+h*Vy+(h/6.0d0)*(ky1+ky2+ky3)
Vx=Vx+(kx1+kx2+kx2+kx3+kx3+kx4)/6.0d0
Vy=Vy+(ky1+ky2+ky2+ky3+ky3+ky4)/6.0d0
X(N)=0; Y(N)=0; Vx(N)=0; Vy(N)=0
PRINT *, t, X(2), Y(2)
END SUBROUTINE RK4
[QUOTE]
[QUOTE]
4 0 500 0.01
-1.01723985460 -0.206180767821 0.005269822149 -0.014839159853
-0.98323814725 -0.176124724912 0.002751137188 -0.016991888551
-0.98545358338 -0.176940858196 0.002966627399 -0.017595568343
0 0 0 0
1.8
3.98600440d14
4.902798d12
1.32712440018d20
[QUOTE]
6 réponses
Votre code semble bon...
Depuis le temps, avez vous trouvez la solution ?
La mécanique celeste m'intéresse pas mal, put être pourriez vous poster si jamais vous avez trouvé votre "erreur" quelques cliché photographique et votre code ?
Depuis le temps, avez vous trouvez la solution ?
La mécanique celeste m'intéresse pas mal, put être pourriez vous poster si jamais vous avez trouvé votre "erreur" quelques cliché photographique et votre code ?
Bonjour, oui mon programme fonctionne, et très bien d'ailleurs, il n'y avait pas qu'une seule erreur, cependant comme il s'agit d'un projet à rendre pour l'obtention de la licence, je ne suis pas certain de vouloir le mettre sur internet, car une simple recherche permettrait de tomber dessus.
Je pourrais éventuellement poster mes résultats sachant que je calcule la trajectoire de la Terre de la lune et de l'astéroide Apophis que nous croiserons en 2029.
Je pourrais éventuellement poster mes résultats sachant que je calcule la trajectoire de la Terre de la lune et de l'astéroide Apophis que nous croiserons en 2029.
Oui, je pense que cela peut être intéressant.
Personnellement, je suis un amateur de ce genre de travaux. Peut être pourriez vous me l'envoyer par MP ?
Parce que obtenir des résultats sans comprendre la méthode qui y amène n'est pas très intéressant.
D'ailleurs, bonne chance pour votre licence d'Astrophysique !
Personnellement, je suis un amateur de ce genre de travaux. Peut être pourriez vous me l'envoyer par MP ?
Parce que obtenir des résultats sans comprendre la méthode qui y amène n'est pas très intéressant.
D'ailleurs, bonne chance pour votre licence d'Astrophysique !
Je serai ravi de vous envoyer mes travaux, mais je ne pourrai vous les donner qu'apres mes examens, donc dans 2 semaines, etes vous d'accord?
Je precise que avec mon programme vous pouvez calculer la trajectoire de N corps en gravitation en comparant la methode RK4 et celle de Runge Kutta Fehlberg.
(le programme dans mon poste n'est qu'une partie)
Et il sagit d'une licence de physique, merci beaucou !
A bientot
Je precise que avec mon programme vous pouvez calculer la trajectoire de N corps en gravitation en comparant la methode RK4 et celle de Runge Kutta Fehlberg.
(le programme dans mon poste n'est qu'une partie)
Et il sagit d'une licence de physique, merci beaucou !
A bientot
Vous n’avez pas trouvé la réponse que vous recherchez ?
Posez votre question
Oui je veux bien. Personnellement, je fais un projet un pu similaire, dans le cadre de ma thèse en Astrophysique. C'est pour ça que cela m'intéressait au plus haut point.
Si vous dites dans deux semaines, pourquoi pas. J'aurais bien aimé jeter un petit coup d'oeil avant, histoire de pouvoir comparer avec mon programme fait sous C++ que je trouve vraiment très lourd ...
Si jamais vous avez le temps avant, cela m'arrangerait, j'ai un stage de fin d'étude à rendre pour le 3 mai.
Merci en tout cas pour vos réponses !
Si vous dites dans deux semaines, pourquoi pas. J'aurais bien aimé jeter un petit coup d'oeil avant, histoire de pouvoir comparer avec mon programme fait sous C++ que je trouve vraiment très lourd ...
Si jamais vous avez le temps avant, cela m'arrangerait, j'ai un stage de fin d'étude à rendre pour le 3 mai.
Merci en tout cas pour vos réponses !
