Cette page est destinée aux débutants en programmation et en calcul numérique: étudiants en sciences, élèves de classes prépa. et tous ceux qui veulent s'initier aux joies de la physique numérique. Elle est orientée FORTRAN, parce que c'est encore un des langages de programmation enseigné en fac de sciences. C'est aussi et surtout un des principaux langages utilisés dans les labo. Pour ceux qui préfèrent le C (et ses extensions C++), je les invite dans la page Initiation au C.
Je ne suis pas certain que les étudiants arrivant en taupe ou en fac de sciences aient entendu parler de FORTRAN avant d'assister à leur premier cours d'informatique. Ce n'est pas étonnant. Aujourd'hui, des langages comme C, C++ et même Java sont bien plus "visibles". FORTRAN est pourtant le langage le plus utilisé dans la communauté scientifique (en 2006 !), bien que C/C++ se rapproche. Les raisons de cette longévité sont nombreuses:
Si vous êtes ici, c'est que vous avez entendu parler de FORTRAN ou même êtes vous entrain de "sécher" sur votre premier programme! Je vous propose donc d'examiner les principales caractéristiques de ce langage et d'en étudier quelques applications simples au calcul scientifique.
Cette page n'a pas la prétention de faire du lecteur un programmeur confirmé mais de donner quelques indications pour aider les débutants en FORTRAN. Nous n'aborderons que les notions essentielles pour écrire des programmes de physique numérique. Elles seront introduites à mesure du besoin.
Le premier compilateur FORTRAN fut conçu et implémenté par John Backus (IBM Lab, San José California) en 1954. Il s'agissait du FORTRAN I. L'objectif de Backus était de concevoir un langage de programmation qui permettrait de réaliser des programmes en utilisant une syntaxe la plus proche possible de la syntaxe mathématicienne. Le terme FORTRAN traduit cette ambition, puisqu'il signifie FORmula TRANslation.
FORTRAN I était un peu primaire, mais quand même très supérieur aux assembleurs et langages machines utilisés au début des années 50. Il évolua rapidement en FORTRAN II et III.
FORTRAN prit son essor dans les années 60, surtout à partir de 1966, date à laquelle l'ANSI (organisme de standardisation américain) publia la norme ANSI66 qui normalisa le FORTRAN IV. A partir de cette date, tous les compilateurs utilisaient la même syntaxe. Un chercheur pouvait donc utiliser un code source écrit en FORTRAN IV sur n'importe quelle machine (bien qu'à l'époque, en calcul scientifique, en dehors des machines IBM, il n'y avait pas vraiment de salut...)
Par la suite, FORTRAN subit plusieurs évolutions, marquées par la publication par l'ANSI de normes:
Le FORTRAN est donc un langage qui a un demi-siècle d'existence, mais qui reste un des plus utilisés en calcul scientifique (malheureusement, dirons certains!).>
Comme tous les programmes, un programme FORTRAN est un ensemble organisé d'instructions données par le programmeur à la machine. Il n'est pas dans mon intention de faire ici un cours d'informatique (il en existe de très bon sur le Web), je supposerai donc que le lecteur sait faire la différence entre un code source (écrit avec un éditeur, dans un langage de programmation, ici le FORTRAN) d'un code exécutable (construit par le compilateur et le linker, en langage machine). Nous parlerons ici essentiellement de code source (ou programme source, on dit aussi "source").
Un programme FORTRAN, du moins son code source, le plus simple est organisé comme suit:
PROGRAM ProgName
<liste d'instructions>
END
Les mots PROGRAM et END sont les mots clé FORTRAN qui délimitent le programme principal. Le nom de ce programme sera ProgName.
La liste d'instructions bornée par PROGRAM et END forme un bloc fonctionnel. Tel quel (en enlevant <liste d'instructions>, tout de même!), ce bloc fonctionnel forme un programme FORTRAN tout à fait correct qui pourrait être compilé et exécuté. Evidement, il ne ferait pas grand chose..
On peut adjoindre à ce programme principal deux autres types de blocs fonctionnels:
SUBROUTINE SubName (Arguments)
<liste d'instructions>
END
et
type FUNCTION FunctionName (Arguments)
<liste d'instructions>
END
Ces deux blocs sont utilisés pour alléger le programme principal de calculs ou de traitements répétitifs. Nous verrons plus loin comment les utiliser.
Un programme FORTRAN classique ressemble donc à quelque chose comme ça:
PROGRAM ProgName
< liste d'instructions >
END
SUBROUTINE SubName (Arguments)
<liste d'instructions>
END
type FUNCTION FunctionName (Arguments)
<liste d'instructions>
END
Note: il est possible de disposer les blocs SUBROUTINE et FUNCTION avant le bloc PROGRAM. C'est une affaire d'habitude et de convenance personnelle. Une seule chose compte: la lisibilité de votre code (on y reviendra dans les exemples plus loin). Par exemple, évitez d'écrire un bloc PROGRAM de 5 000 lignes!
Il reste maintenant le plus gros du travail: identifier et comprendre ce que l'on entend par <liste d'instructions>
Comme tous les langages informatiques, FORTRAN utilise des données et des opérateurs. Nous commencerons par examiner les données, leur type et leur forme.
Les principaux types de donnée simples (nous verrons les tableaux plus loin) que nous utiliserons sont:
Il en existe d'autres (CHARACTER,LOGICAL) que nous verrons plus loin.
Selon le type de la constante son écriture sera différente (attention, c'est un piège classique):
On peut (c'est même une bonne habitude) nommer ses constantes en utilisant le mot clé PARAMETER comme dans :
REAL pi
PARAMETER (pi = 3.14)
Notez que j'ai d'abord précisé que pi est un REAL
Pour construire une variable dans un programme, on énonce tout simplement le type de la variable et son nom. Par exemple:
INTEGER toto
COMPLEX z
REAL x,y,z
On définit les variables d'un programme, d'une SUBROUTINE ou d'une fonction entre le mot clé de début de bloc (PROGRAM,par exemple) et la première instruction exécutable.
Il y a plusieurs règles à connaître à propos des variables, dont certaines sont spécifiques à FORTRAN:
FORTRAN supporte les opérateurs arithmétiques suivants:
Attention à la conversion de type dans une opération! Par exemple, si vous écrivez a**2, le résultat sera du type de a, car 2 est une constante entière. Si vous écrivez b**2 - 4*a*c et que a et c sont REAL alors que b est DOUBLE PRECISION, le résultat sera DOUBLE PRECISION. Mais le calcul 4*a*c sera fait en REAL. Ce genre de piège, qui se traduit par une perte de précision, est assez courant. Attention donc à la cohérence de vos typages de variables.
FORTRAN supporte un grand nombre de fonctions mathématiques: SIN, COS ( en REAL, DOUBLE PRECISION et même COMPLEX), TAN, ATAN, ASIN, ACOS, LOG10, LOG, EXP, SQRT, etc.. On les retrouvera dans nos programmes.
Signalons aussi les fonctions MIN et MAX, la fonction ABS, les fonctions INT (partie entière) et MOD (reste de la division entière), dont les noms varient selon le type des variables.
Voyons les instructions qui nous permettent d'afficher un résultat sur l'écran ou de saisir un paramètre depuis notre clavier. A ce stade nous en resterons à ce genre d'entrée/sortie. Nous parlerons des fichiers plus loin.
Pour lire les informations saisies depuis un clavier, il existe deux méthodes:
- la simple (et toujours applicable sauf aux chaînes de caractères), qui est la lecture non formatée. La donnée saisie prend le type de la variable dans laquelle on la stocke. Cela s'écrit:
REAL x
READ(*,*) x
- un peu plus compliquée, la lecture formatée. On indique au programme un certaine nombre d'informations en particulier le nombre de chiffres avant et après la virgule et bien d'autres choses. cela s'écrit, par exemple:
REAL x
READ(*,'(e12.5)') x
ce qui indique qu'on lit un réel en précisant que l'on veut 12 caractères, dont 5 chiffres pour la mantisse. Nous ne nous servirons pas de cette possibilité pour l'instant.
Le premier * dans le READ indique le canal de lecture. Si on n'en utilise qu'un seul (le clavier, par exemple), il est inutile de le préciser. Sinon, on indique un nombre entier compris entre 10 et 99, qui identifie le canal pour chaque opération d'entrée/sortie (voir plus loin la cas d'utilisation des fichiers).
Les principes sont les mêmes. On utilisera l'écriture non formatée sur l'écran. Pour ce faire, la bonne instruction est:
WRITE(*,*) x
Deux instructions à connaître pour saisir une information depuis le clavier de manière sympathique:
WRITE(*,'(a,$)') 'Rayon de la sphère : '
READ(*,*) rayon
Notez la syntaxe du WRITE. Comme je l'ai dit plus haut, la lecture et l'écriture non formattée d'une chaîne de caractères est interdite. Vous utiliserez donc le format indiqué.
Que faut-il savoir de plus pour commencer à écrire un petit programme FORTRAN? Quelques spécificités propres à FORTRAN....
Historiquement, le code FORTRAN était saisi sur des cartes perforées. Ce qui induisait certaines contraintes qui ont survécues jusqu'en FORTRAN 77. Signalons les principales:
J'ai dit tout à l'heure que FORTRAN n'était pas sensible à la casse (minuscule ou majuscule). Mais vous avez sans doute remarqué aussi que j'écrivais tous les mots clé en majuscule. Ce n'est pas obligatoire, mais c'est une coutume en FORTRAN (elle provient des premiers standards de FORTRAN)! Et puis, cela nous distingue de ces vulgaires programmeurs de C ou java....
Ces règles ont évolué avec les dernières moutures de FORTRAN (90 et 95). Dans ces versions, elles sont moins contraignantes. Je les cite parce que j'utilise ici un compilateur 77, qui est encore très largement utilisé, mais aussi parce qu'elles font partie du folklore...
Tout ça à l'air un peu compliqué et ringard, mais on s'y fait très vite. D'autant que les éditeurs de textes, pour la plupart, connaissent et gèrent les règles de FORTRAN.
Pour notre premier programme, nous calculerons le volume d'une sphère en fonction de son rayon que nous indiquerons au programme. J'ai préféré un programme de calcul, même trivial, plutôt que le classique "Hello the word" de C ou Java. Nous sommes des scientifiques que diable! Et puis, ce genre de programme nous permet de mettre en oeuvre les notions introduites ci-dessus.
Le code source de ce programme:
C Programme de calcul du volume d'une sphere
PROGRAM VolumeSphere
IMPLICIT NONE
C Definition des parametres
REAL PI
PARAMETER (PI=3.14159265)
C Definition des variables
REAL rayon, volume
C Saisie du rayon de la spher
WRITE(*,'(a,$)')'Rayon de la sphere : '
READ(*,*) rayon
C Calcul du volume
volume = (4*PI/3)*rayon**3
C Affichage
WRITE(*,*) 'Volume de la sphere = ', volume
END
Il commence par une ligne de commentaire indiquant l'objet du programme. Bonne habitude à prendre que de commenter son code.
Puis le mot clé PROGRAM suivi de son nom (Attention, là aussi seuls les 6 premiers caractères sont significatifs).
Notez ensuite l'instruction IMPLICIT NONE, qui m'oblige à déclarer toutes mes variables. Ainsi, je n'en n'oublierai pas.
Viennent les déclarations des constantes et autres paramètres: déclarer ici tous vos paramètres symboliques (PI, constante de Plank, densité de l'eau, distance terre-lune, etc...). Attention à la syntaxe et à déclarer le type des constantes.
Déclaration des variables: simple et de bon goût.
Et enfin le corps du programme:
qui se termine par END, comme il se doit.
Voilà un tout petit programme, mais qui contient tous les germes d'un grand code. Voyons comment le saisir, le compiler et l'exécuter.
Pour créer et exécuter un programme FORTRAN il faut un minimum d'outils. A commencer par un ordinateur.
Là, survient mon premier dilemme: dois-je m'orienter vers le monde Windows ou Linux? Personnellement je préfère et je travaille sous Linux (Mandriva). Mais je suppose que la plupart d'entre vous sont plutôt sous Windows. Il faut savoir que FORTRAN se moque du système d'exploitation (OS): la syntaxe est la même sous Linux ou Windows. Seuls les outils différent, et encore! Les environnements de développement FORTRAN existent souvent pour les deux OS. Quant au compilateur, j'utilise G77, un compilateur disponible à la fois sous Windows et sous Linux.
Bref, j'ai opté, pour les exemples de cette page, de travailler sous Windows. Les linuxiens me pardonneront (je suis des leurs!) et se consoleront en se disant que tout ce que j'écris ici est directement transposable (sauf les commandes système, évidement... mais les linuxiens sont rarement des informaticiens débutants!).
Donc, il faut rassembler quelques programmes, que nous allons trouver gratuitement sur le net:
Puis je vous laisse installer ces deux programmes.
Il serait pratique de créer un répertoire pour nos petits travaux de programmation. Le mien s'appelle numlab\labphysique.
Lancer le programme VFort. La vue principale de l'IDE s'affiche.
Pour saisir notre programme, il faut:
Et maintenant on peut compiler et linker le fichier source, en faisant Project/Build. Si vous n'avez fait aucune erreur de syntaxe, vous devez voir apparaître les messages suivants (ou approchant..) dans la fenêtre inférieure:
Compiling (command: g77 -c) ...
E:\NumLab\LabPhysique\VolumeSphere\VolumeSphere.for
Linking (command: g77) ...
Code generation: VolumeSphere.exe
Tout c'est bien passé. Sinon, corrigez vos erreurs indiquées par le compilateur et recommencez à l'étape Project/Build.
Vous noterez que le programme executable porte le nom donné dans l'instruction PROGRAM. Il se trouve dans le sous-répertoire Output de votre répertoire de travail (ici NumLab\LabPhysique\VolumeSphere\).
Pour le lancer, sous Windows, il faut ouvrir une fenêtre de commande (cmd.exe) et déclarer le répertoire output comme répertoire courant. Puis tapez VolumeSphere et Return. Le programme s'exécute.
Un conseil, gardez la fenêtre de commande ouverte, cela vous évitera des opérations fastidieuses.
Et voilà, votre premier programme fonctionne! Passons à des choses un peu plus compliquées...
Dans cette partie un peu longue, on va examiner:
Pour aborder les conditions, il faut d'abord parler d'expression logique.
Il existe en FORTRAN un type de variable LOGICAL. Une variable de type LOGICAL peut prend deux valeurs: .TRUE. ou .FALSE. Il faut savoir que .TRUE. vaut 1 et .FALSE. vaut 0.
FORTRAN définit les expressions logiques suivantes, dans lesquelles x et y sont deux variables :
Il définit aussi les opérateurs logiques unaires:
Munis de ces expressions logiques, on peut définir les instructions conditionnelles suivantes:
IF (condition) THEN
<liste d'instructions>
ENDIF
La signification de ce jeu d'instruction est simple: si l'expression condition est vraie, par exemple x .LT. 0, alors la liste d'instructions contenue entre le IF et le ENDIF est exécutée. Sinon, le programme continu et exécute l'instruction immédiatement suivante le ENDIF.
Il existe plusieurs variantes d'instructions conditionnelles. Par exemple, la condition alternative:
IF (condition) THEN
<liste d'instructions 1>
ELSE
<liste d'instructions 2>
ENDIF
Si l'expression condition est vraie, alors la liste d'instructions 1 est exécutée (et pas la liste d'instructions 2). Si l'expression condition est fausse, c'est la liste d'instruction 2 qui est exécutée. Il s'agit en fait de la forme complète de l'instruction conditionnelle IF.
Le langage C possède une instruction très sympathique, le 'switch', qui permet de faire des traitements différents en fonction de la valeur d'une variable. Cette instruction n'existe pas en FORTRAN 77, mais on peut la simuler en utilisant le jeu d'instructions suivant:
IF (condition1) THEN
<liste d'instructions 1>
ELSE IF (condition2) THEN
<liste d'instructions 2>
ELSE IF (condition3) THEN
<liste d'instructions 3>
ELSE
<liste d'instructions n>
ENDIF
Vous en devinez l'utilisation! Supposons que je doive exécuter un traitement différent selon qu'une variable x est égale à 1, 2, 3 ou différente de 1, 2, et 3. Je poserai:
IF (x .EQ. 1) THEN
<liste d'instructions 1>
ELSE IF (x .EQ. 2) THEN
<liste d'instructions 2>
ELSE IF
(x .EQ. 3) THEN
<liste d'instructions 3>
ELSE
C cas de x différent de 1,2 et 3
<liste d'instructions 4>
ENDIF
Si x = 1, seule la liste d'instructions 1 est exécutée. Si x=2, seule la liste d'instructions 2 est exécutée, etc...
C'est un jeu d'instruction très employé en traitement des données, mais qui peut générer des programmes un peu fouillis! Les commentaires sont indispensables.
Sans doute, ce sont les instructions les plus utilisées en calcul. FORTRAN 77 comporte deux types de boucles.
Sa syntaxe
DO var= debut, fin, pas
<liste d'instructions>
ENDDO
où var est la variable de boucle, qui évolue de debut à fin selon un incrément de pas. Notons que pas peut être omis, dans ce cas, il vaudra 1.
Son mode de fonctionnement
La liste d'instructions est exécutée un certain nombre de fois, déterminé par la valeur de debut et de fin. En principe, on ne doit pas sortir de la boucle autrement que lorsque var > fin, qui est la sortie naturelle par le ENDDO. C'est techniquement possible mais à proscrire rigoureusement.
var, debut, fin et pas sont des variables de type INTEGER.
pas peut être négatif. Dans ce cas debut doit être plus grand que fin
Je vous laisse découvrir ce qui se passe si fin > debut, selon la valeur de pas (>0 ou <0)
Exemple
somme = 0
DO i = 1, 100
somme = somme + i
ENDDO
Cette boucle calcule la somme de 100 premiers entiers.
Curiosité: on peut écrire une boucle infinie en écrivant:
DO
<liste d'instructions>
ENDDO
Dans ce cas, on sort de la boucle en utilisant l'instruction exit en un lieu approprié de la liste d'instructions. Ce genre de manip est à éviter dans un premier temps...
Sa syntaxe
DO WHILE(condition)
<liste d'instructions>
ENDDO
où condition est une expression logique.
Son mode de fonctionnement
Le programme entre dans la boucle si condition = .TRUE. Il reste dans la boucle tant que (i.e. la liste d'instruction est exécutée autant de fois que) condition vaut .TRUE.
Lorsque condition vaut .FALSE. l'exécution se poursuit à la première instruction qui suit le ENDDO. Cela implique qu'il existe dans la liste d'instructions une instruction susceptible de faire passer condition à .FALSE. Sinon le programme bouclera indéfiniment, ce qui rarement le but recherché (mais cela arrive parfois dans les traitements d'acquisition).
Exemple
C déclaration des variables
INTEGER n
DOUBLE PRECISION somme, epsilon
LOGICAL fin
C Initialisation des variables
epsilon = 1.d-10
somme = 0d0
fin = .FALSE.
C Boucle de calcul. La condition est .NOT. fin de telle sorte que la condition soit vraie à l'entrée
DO WHILE (.NOT. fin)
somme = somme + 1d0/n
n = n +1
C On sort de la boucle si fin devient TRUE
fin = (1d0/n .LT. epsilon*somme)
ENDDO
Cette boucle calcule la suite des inverses jusqu'à ce que le terme 1/n soit inférieur à epsilon fois la somme partielle.
Attention, dans cet exemple à la gestion de la condition. C'est un cas très répandu mais un peu tordu. On pose fin = .FALSE. et la condition de test de la boucle est .NOT. fin. Cela signifie qu'on entre dans la boucle, car au premier tour (.NOT. fin) = .TRUE
Lorsque la condition mathématique de convergence est remplie fin devient .TRUE. Le test de boucle nous donne donc .NOT. fin = .FALSE. et l'on sort de la boucle. Tout ça pour que le test est un sens mathématique usuel...
Le dernier élément indispensable pour faire du calcul numérique! Je ne sais pas ce que l'on ferait sans tableau... J'ai l'habitude, comme beaucoup de mes collègues de les appeler des 'vecteurs' ou des 'matrices'. Vous allez vite comprendre!
On déclare un tableau en FORTRAN de la façon suivante:
type Nom(dim1, dim2, dim3, ...)
type désigne le type des variables stockées dans le tableau nommé Nom. Le nombre de constantes dim indique la dimension du tableau. Un vecteur est donc un tableau de 1 dimension...
dim1, dim2, ... déterminent la taille du tableau. Ce sont des constantes entières. FORTRAN 77 ne gère pas les tableaux dynamiques, ni même le redimensionnement dynamique des tableaux. Il est de bon ton de fixer les valeurs dim1, dim2, ... sous forme de paramètres.
Par exemple:
INTEGER n
PARAMETER (n = 1000)
REAL toto(n
On vient de construire un vecteur de 1000 réels. Notons qu'un tel tableau occupe 4000 octets. Où encore une matrice carrée de 1000*1000 DOUBLE PRECISION:
DOUBLE PRECISION matcarre(n,n)
qui occupe elle 1000*1000*8 octets.....
Voilà pour la déclaration des tableaux. Voyons maintenant l'utilisation de ceux-ci:
Pour stocker une variable dans un tableau, il faut préciser la position du tableau dans laquelle on veut écrire. Par exemple :
i = 2
toto(i) = cos(x)
ou
i = 1
j = 1
matcarre(i,j) = toto(i)
Cette dernière écriture signifie: je prends la valeur contenue dans la case i du tableau toto et je l'assigne à la case (i,j) de matcarre.
Voici quelques exemples superclassiques d'utilisation de tableaux.
La somme de deux matrices carrées:
INTEGER n
PARAMETER (n = 1000)
DOUBLE PRECISION m1(n,n), m2(n,n), m3(n,n)
DO i=1, n
DO j=1, n
m3(i,j) = m1(i,j) + m2(i,j)
ENDDO
ENDDO
Le produit scalaire de deux vecteurs:
INTEGER n
PARAMETER (n = 1000)
DOUBLE PRECISION u(n), v(n)
DOUBLE PRECISION Scalar
Scalar = 0d0
DO i=1, n
Scalar = Scalar + u(i)*v(i)
ENDDO
Les fonctions et les subroutines sont deux notions importantes de programmation. On emploie une fonction ou une SUBROUTINE lorsqu'on veut coder un traitement qui est souvent répété dans un programme. Il peut s'agir:
Les subroutines servent aussi fréquemment à matérialiser la décomposition fonctionnelle d'un problème. On découpe un problème compliqué en plusieurs blocs fonctionnels différents. Cela permet de mieux comprendre le code, mais aussi de pouvoir le faire réaliser par plusieurs programmeurs. Dans ce cas, on est vite confronté au besoin de partager des données entre subroutines.
ATTENTION: une fonction (resp. une SUBROUTINE) est un bloc fonctionnel différent du programme principal. Chaque bloc fonctionnel possède ses propres variables locales qui ne sont pas connues des autres. Ainsi, une variable peut être nommée toto dans le programme principal et une autre variable porter le même nom dans une fonction (resp. une SUBROUTINE) sans pour autant qu'il s'agisse de la même variable.
De même, le nom des paramètres dans l'en-tête d'une fonction (resp. d'une SUBROUTINE) n'a absolument pas de signification en dehors de cette fonction (resp.SUBROUTINE).
La structure générale d'une fonction est la suivante:
type FUNCTION NomFonction(par1, par2, ...)
type par1
type par2
C déclaration des variables locales
C Traitement effectué par la fonction
<liste d'instructions dont une instruction sera de la forme NomFonction = ... >
RETURN
END
Depuis un programme, on appelle une fonction comme n'importe quelle fonction prédéfinie de FORTRAN. Par exemple:
PROGRAM toto
REAL R, V
R = 1.
C Notez que la variable dans toto s'appelle R, alors que la paramètre formel de VolumSphere s'appelle rayon
V = VolumeSphere(R)
END
C définition de la fonction VolumeSphere
REAL FUNCTION VolumeSphere(rayon)
REAL rayon
C Definition des parametres - Il s'agit de variables locales, PI est inconnu du programme toto
REAL PI
PARAMETER (PI=3.14159265)
VolumeSphere = (4*PI/3)*rayon**3
RETURN
END
Attention, lorsque vous appelez une fonction depuis votre programme (on va voir comment), vous devez vous assurer que les paramètres utilisés dans l'appel sont bien strictement du même type que les paramètres déclarés dans l'en-tête de définition de la fonction. Sinon, bug assuré!
La structure générale d'une subroutine est la suivante:
SUBROUTINE NomSub(par1, par2, ...)
type par1
type par2
....
C déclaration des variables locales
C Traitement effectué par la fonction
<liste d'instructions>
RETURN
END
Depuis un bloc fonctionnel (programme, autre SUBROUTINE, fonction), on appelle une SUBROUTINE par un CALL, comme dans l'exemple ci-dessous. ATTENTION, FORTRAN 77, contrairement à FORTRAN 90, ne supporte pas les appels récursifs (une SUBROUTINE ne peut pas s'appeler elle-même).
PROGRAM toto
REAL R,V
R = 1.
CALL VolumeSphere(R,V)
END
C définition de la SUBROUTINE VolumeSphere
SUBROUTINE VolumeSphere(rayon, volume)
REAL rayon, volume
C Definition des parametres - Il s'agit de variables locales, PI est inconnu du programme toto
REAL PI
PARAMETER (PI=3.14159265)
volume = (4*PI/3)*rayon**3
RETURN
END
Attention, comme pour les fonctions, lorsque vous appelez une SUBROUTINE depuis votre programme, vous devez vous assurer que les paramètres utilisés dans l'appel sont bien strictement du même type que les paramètres déclarés dans l'en-tête de définition de la SUBROUTINE. Sinon, bug assuré!
Nous verrons dans l'exemple du programme Euler l'utilisation de ces différents blocs fonctionnels.
C'est un vaste sujet! Nous nous limiterons sur cette page à examiner les instructions nécessaires à créer, lire et écrire les fichiers qui conviennent pour stocker nos résultats de calcul ou nos données d'entrée.
On ne considérera donc que les seuls fichiers séquentiels formattés (autrement appelés fichier texte). Ces fichiers ont la propriété intéressante de pouvoir être lus et modifiés avec n'importe quel éditeur de texte.
La syntaxe de l'instruction pour ouvrir un fichier est:
OPEN(unité,
FILE= chaine,
FORM = chaine,
STATUS = chaine
ERR = numlabel)
où
Exemple : OPEN(10, FILE = 'd:\NumLab\LabPhysique\data.txt') ouvre le fichier data.txt, qu'il existe ou non (par défaut). En cas d'erreur, le programme se plante...
La syntaxe de l'instruction de lecture est:
READ(unité,
format,
ERR = numlabel,
END = numlabel) liste de données
où
Exemple : READ (10,*,ERR=90, END =100) toto lit une donnée dans le fichier data.txt et la stocke dans la variable toto.
La syntaxe de l'instruction de lecture est:
WRITE(unité,
format,
END = numlabel) liste de données
où
Exemple : WRITE (10,*) toto écrit la variable toto dans le fichier data.txt.
Facile... On utilise l'instruction CLOSE(unité).
En principe, tous les fichiers sont fermés automatiquement à la fin d'un programme. Mais c'est une bonne habitude à prendre que de les fermer par un CLOSE lorsqu'on a fini de s'en servir.
Supposons que nous ayons stocké les résultats d'un calcul dans deux vecteurs x et y de dimension n. Nous voulons sauvegarder ces valeurs dans un fichier pour les plotter avec un logiciel comme Matlab, Scilab ou GnuPlot.
Voilà le bout de code correspondant. On supposera que le fichier data.txt n'existe pas.
INTEGER n, i
PARAMETER (n = 100)
DOUBLE PRECISION x(n), y(n)
CHARACTER*80 nomfic
nomfic = 'data.txt'
OPEN(10, FILE=nomfic)
DO i=1,n
WRITE(10,*) REAL(x(i)), REAL(y(i))
ENDDO
CLOSE(10)
END
Il y a plusieurs nouveautés dans ce code:
Nous voilà maintenant en possession de tous les éléments pour réaliser un programme de calcul en FORTRAN.
Pour illustrer cette présentation, j'ai choisi la méthode d'Euler. C'est la méthode la plus simple, à défaut d'être la plus efficace, de résolution d'une EDO (équation différentielle ordinaire ou ODE comme disent les anglophones). Cette méthode est présentée dans tous les cours d'analyse numérique. Rappelons là brièvement.
Soit à résoudre une EDO avec une condition initiale de type:
y' = f(x, y)
y(x0) = y0
La méthode d'Euler consiste à estimer l'accroissement de y entre les abscisses xi et xi+1 par un développement de Taylor d'ordre 1 (on dit que c'est une méthode à un pas):
yi+1 = yi + h*f(xi, yi) où h = xi+1 - xi
Cette méthode n'est pas très précise. De plus, elle génère et amplifie les erreurs. Mais dans les cas simples, si f n'est pas trop variable et h très petit, elle sera suffisante en première approximation.
Pour traiter des EDO d'ordre supérieur, on procéde à un changement de variables qui aboutit à un système d'équations différentielles. Par exemple, si l'on considère l'EDO:
y'' = f(x,y,y')
avec y(x0) = y0 et y'(x0) = y'0
On peut la décomposer en un système de 2 équations du premier ordre en posant y' = u(x), d'où:
y' = u(x)
u'(x) = f(x,y,u)
et résoudre ce système par la méthode d'Euler.
Nous allons illustrer l'utilisation de la méthode d'Euler par l'étude de la dynamique du pendule simple. Comme vous le savez sans doute, l'équation différentielle de son mouvement est :
d2(theta)/dt2 = -(g/l)sin(theta)
où l = longueur du fil de suspension. Je ne présente pas g...
Pour appliquer la méthode d'Euler, je vais décomposer cette EDO en un système à 2 équations en posant:
d(theta)/dt = omega la vitesse angulaire
a(theta) = d(omega)/dt l'accelération angulaire, qui vaut -(g/l)*sin(theta)
L'application du schéma d'Euler me donne les équations:
thetai+1 = thetai + h*omegai (1)
omegai+1 = omegai + h*ai = omegai + h*(-(g/l)sin(thetai)) (2)
Nous retrouverons tout ce petit monde dans le programme ci-dessus:
C Programme de demonstration de l'usage de la methode d'EULER
C Calcul du mouvement d'un pendule simple
C Dominique Lefebvre mars 2006
PROGRAM PenduleEuler
IMPLICIT NONE
C Declaration des constantes
REAL g, PI, n, P
PARAMETER (g = 9.81, PI = 3.141592654)
PARAMETER (n = 2) ! Nb de cycles de calcul
PARAMETER (P = 1000) ! Nb de pas de calcul par cycle
C Declaration des variables
INTEGER i
REAL l, Theta, h
REAL X1(n*P), X2(n*P), t(n*P) ! resp. l'angle, la vitesse angulaire, le temps
C Initialisation des variables
l = g/(4*PI*PI) ! je choisis une longueur arbitraire qui m'arrange
Theta = 2*PI*SQRT(l/g) ! la periode vaut 1 pour l'approximation des petits angles
h = Theta/P ! pas temporel pour le calcul
C Declaration des conditions initiales d'integration
X1(1) = 60*PI/180 ! angle initial en radian (par exemple...)
X2(1) = 0 ! vitesse angulaire initiale
t(1) = 0.0 ! Origine des temps
C Calcul du mouvement selon l'algorithme d'Euler
DO i=1, n*P-1
t(i+1) = t(i)+h ! incrementation du temps
C Le systeme differentiel descriptif du pendule
f1 = X2(i)
f2 = -(g/l)*SIN(X1(i))
C L'algorithme d'EULER
X1(i+1) = X1(i)+f1*h ! Equation (1) du systeme
X2(i+1) = X2(i)+f2*h ! Equation (2) du systeme
ENDDO
C Sauvegarde des donnees dans le fichier Euler.dat
C On sauvegarde les données pour tracer la courbe theta = f(t) mais on pourrait
C tracer d'autre courbes, en particulier omega = f(theta). Il suffit de bien
C choisir les données à stocker
OPEN(10, FILE='Euler.dat')
DO i=1, n*P
WRITE(10,*) REAL(t(i)),
REAL(X1(i)*180./PI) ! notez la conversion en degres<
ENDDO
CLOSE(10)
END
ou encore sous une forme plus structurée et plus facilement modifiable:
C Programme de demonstration de l'usage de la methode d'EULER
C Calcul du mouvement d'un pendule simple
C Dominique Lefebvre mars 2006
PROGRAM PenduleEuler
IMPLICIT NONE
C Declaration des constantes
REAL g, PI, n, P
PARAMETER (g = 9.81, PI = 3.141592654)
PARAMETER (n = 2) ! Nb de cycles de calcul
PARAMETER (P = 1000) ! Nb de pas de calcul par cycle
C Declaration des variables
INTEGER i
REAL l, Theta, h
REAL X1(n*P), X2(n*P), t(n*P) ! resp. l'angle, la vitesse angulaire, le temps
C declaration des zone commune
COMMON /Systeme/l
C Initialisation des variables
l = g/(4*PI*PI) ! je choisis une longueur arbitraire qui m'arrange
Theta = 2*PI*SQRT(l/g) ! la periode vaut 1 pour l'approximation des petits angles
h = Theta/P ! ! pas temporel pour le calcul
C Declaration des conditions initiales d'integration
X1(1) = 60*PI/180 ! angle initial en radian
X2(1) = 0 ! vitesse angulaire initiale
t(1) = 0.0 ! Origine des temps
C Calcul du mouvement selon l'algorithme d'Euler
DO i=1, n*P-1
t(i+1) = t(i)+h ! incrementation du temps
C Le systeme differentiel descriptif du pendule
CALL Derivee(X1(i), X2(i), f1, f2)
C L'algorithme d'EULER
CALL Euler(X1,X2,f1,f2,i,h)
ENDDO
C Sauvegarde des donnees dans le fichier Euler.dat
C On sauvegarde les données pour tracer la courbe theta = f(t) mais on pourrait
C tracer d'autre courbes, en particulier omega = f(theta). Il suffit de bien
C choisir les données à stocker
OPEN(10, FILE='Euler.dat')
DO i=1, n*P
WRITE(10,*) REAL(t(i)), REAL(X1(i)*180./PI)
ENDDO
CLOSE(10)
END
C Subroutine de calcul du systeme differentiel
C En cas d'utilisation d'un autre système différentiel, seule
C cette subroutine sera modifiee
SUBROUTINE Derivee(x,y,fx,fy)
REAL x,y,fx,fy
COMMON /Systeme/l
REAL g, l
PARAMETER (g = 9.81)
fx = y
fy = -(g/l)*SIN(x)
END
C Subroutine Algorithme d'Euler
SUBROUTINE Euler(x,y,f1,f2,i,dt)
REAL x(*), y(*), f1, f2, dt
INTEGER i
x(i+1) = x(i)+f1*dt ! on retrouve ici les formules du schéma d'Euler
y(i+1) = y(i)+f2*dt
END
Quelques remarques à propos de ce dernier programme:
J'utilise GnuPlot 3.7, que vous pouvez très facilement charger, que vous soyez sous Windows ou Linux (www.gnuplot.info). Je n'ai pas encore essayé la dernière version 4!
GnuPlot est un logiciel destiné au tracé et à l'impression des données en 2D et 3D. Ce n'est sans doute pas le meilleur, mais c'est le plus répandu. Il possède la majeure partie des fonctions utilisées quotidiennement pour faire du tracé de courbes.
Il n'est pas dans mon intention de faire de grands développements sur GnuPlot ici. Celui-ci possède une fonction d'aide (tapez help après le prompt gnuplot) qui vous donnera toutes les informations que vous voulez connaître. Nous allon voir ici les principales commandes pour tracer une courbe.
Nous partirons du principe que nos données sont stockées dans un fichier texte, par exemple Euler.dat, produit par un programme de calcul.
Nous voulons tracer la courbe correspondante sur une feuille propre avec des légendes, un titre, etc. Voyons les principales commandes nécessaires. Je ne fais pas apparaître le prompt de GnuPlot..:
Voyons comment tracer la courbe correspondante à l'exécution de notre programme Euler. Comme je suis un peu paresseux, je vais stocker toutes mes commandes gnuplot dans un fichier de commandes. Je créé donc un fichier avec un éditeur quelconque (Notepad par exemple), que j'appelle Euler.p. Voici son contenu:
# Fichier de tracé de la courbe theta = f(t)selon la methode d'Euler
set title "Variation de l'angle en fonction du temps"
set xlabel "Temps"
set ylabel "Angle"
plot 'Euler.dat'
Comme vous l'avez deviné, le # désigne une ligne de commentaires
Je place ce fichier dans le répertoire qui contient le fichier Euler.dat, qui doit être le répertoire de travail de GnuPlot.
Maintenant, au prompt de GnuPlot, je tape load 'Euler.p' et j'obtiens une magnifique courbe!
Je vous donne rendez-vous maintenant sur la page des schémas numériques, où nous pourrons aborder des schémas numériques plus élaborés, en particulier la méthode Runge-Kutta, clé de voute de la résolution des équations différentielles ordinaires (EDO) et donc de la simulation.
Si vous souhaitez plus de renseignements, voici quelques références:
Contenu et design par Dominique Lefebvre - tangenteX.com mars 2006 -- Vous pouvez me joindre par ou sur PhysiqueX
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France .