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 (en 2006). 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 - www.tangenteX.com mars 2006
Contact : PhysiqueX ou
Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.