Discussion :

[jphysic]

Bonjour, j'ai un petit problème dans la programmation en fortran90. j'ai fait un code qui calcule les trajectoires d'un système de vortex qui est décrit par le système d'équations différentielles suivant:

dx(i)/dt=Sum  {[omega(i)* (y(j)-y(i))] / [((x(i)-x(j))**2+(y(i)-y(j))**2)]}

dy(i)/dt=Sum  {[omega(i)* (x(i)-x(j))] / [((x(i)-x(j))**2+(y(i)-y(j))**2)]}

la somme se fait sur les i=1,N, sauf les cas ou (i différent de j). et le N étant un naturel différent de zéro.

lorsque je compile (ctrl+F8) et je build (shift+F8) , il n y a pas d'erreurs. Mais lorsque j'exécute une fenêtre d'erreurs apparait, et voici le message:

run-time error M6101: Math
-floating-point error: invalid
Press any key to continue

je comprend pas, alors je demande de l'aide de l'un des spécialistes du fortran car moi , je suis qu'un débutant.

voici le code que j'ai réalisé:

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
PROGRAM VortexEuler
Implicit None 
INTEGER :: N,i,j,k,iter
REAL, Dimension(10,10):: x,y,som1, som2
!DOUBLE PRECISION :: somme1,somme2
DOUBLE PRECISION :: dt,w
!LOGICAL :: fini
! Initialisation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
N=3
iter=10
dt= 1.0e-2
som1(1,1)=0.d0
som2(1,1)=0.d0
w=1.0e-2
x(1,1)=1.d0
y(1,1)= 2.d0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
open(1,File='Vortex euler.txt')
 
do k=1,iter
 do j=1,N
  do i=1,N
    if (i.ne.j) then
 
    som1(i,k)=som1(i,k)+  dt*((w*(y(i,k)-y(j,k)))/(((x(i,k)-x(j,k))**2)+ ((y(i,k)-y(j,k))**2)))
    !
    som1(i,k)=som1(i,k)
    !                                                          
    som2(i,k)=som2(i,k)+ dt*((w*(x(j,k)-x(i,k)))/(((x(i,k)-x(j,k))**2)+ ((y(i,k)-y(j,k))**2)))
        !
    som2(i,k)=som2(i,k)
 
    x(j,k+1)=    x(j,k) + dt* som1(i,k)
    y(j,k+1)= y(j,k) + dt* som2(i,k)
    !
       write(1,*) 'x(j,k+1)=',x(j,k+1),'y(j,k+1)=',y(j,k+1),'iter=',iter
    !
       end if
   enddo
 
  enddo
enddo
 
 
 
close(1)
 
 
END PROGRAM VortexEuler

Merci d'avance.

[LittleOwl]

Bonjour,

Je pense que vous obtenez une erreur d'éxécution arithmétique car votre calcul fait intervenir des valeurs qui n'ont pas été initialisées.
Le premier calcul effectué correspond au cas (k = 1, j = 1, i = 2), le cas (k =1, j = 1, i = 1) n'étant pas calculé puisque i = j. Or les quantités som1(2,1), som2(2,1), y(2,1), x(2,1) n'ont pas été initialisées. Les résultats produits sont alétoires, ce qui est probablement à l'origine de l'erreur d'éxécution au bout d'un certains nombre d'itérations. De mon point de vue, le problème vient de l'algorithme (l'ordre d'imbrication des boucles est-il le bon?)

cordialement.

[François L.]

En effet, il y a le problème de l'initialisation, mais pas que: la boucle sur k amène k à la valeur iter, soit 10. Or, il y a un

Code :

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

x(j,k+1)= x(j,k) + dt* som1(i,k)

dans le code, qui fait que l'on essaiera d'attaquer la 11ème valeur du second indice de x... Donc un futur problème de dépassement de tableau.

Par ailleurs, il y a effectivement un problème d'imbrication des boucles, le calcul de x(j,k+1) ne doit pas être fait dans la boucle sur i, dont le rôle n'est que de calculer la somme. Dans le même ordre d'idées, les variables som ne devraient pas être indicées par i, si elles représentent la somme pour i allant de 1 à N... (l'indice i est muet...). Ou alors je me trompe, i n'est pas l'indice de sommation (mais dans ce cas c'est j, et la remarque s'applique aussi). Vous ne calculez donc visiblement pas ce que vous souhaiteriez calculer...

Ajoutons que numéroter les fichiers à 1 n'est pas une bonne pratique; commencez à partir de 10 (donc:

Code :

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

open(10,File='Vortex euler.txt')

Ca n'est pas le problème actuellement, mais bon...

Enfin, d'un strict point de vue numérique, je n'aime pas la méthode d'Euler explicite, qui a souvent tendance à diverger, sauf à mettre des pas d'intégration très faibles... Et ça, ça peut très bien faire partie du problème actuel...

Je vous conseille de tout reprendre à zéro, en écrivant le système pour un temps donné tout d'abord, puis en l'insérant dans l'intégrateur. Pourquoi ne pas faire une subroutine qui calculerait les dérivées? Ca permettrait d'utiliser n'importe quel intégrateur...

Bon courage.