Comment s'initier par un projet C++( suite) sur editeur de texte dans TERMINAL

Initiation au FORTRAN

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.

Pourquoi le FORTRAN

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, bien que C/C++ se rapproche. Les raisons de cette longévité sont nombreuses:

FORTRAN est un langage structuré, de syntaxe simple et facile à apprendre, qui a su évoluer.
C'est un langage utilisé depuis 40 ans : les programmes et bibliothèques de routines disponibles se chiffrent en millions de lignes. On trouve toujours le bout de code que l'on cherche dans une bibliothèque existante.
C'est un langage universel: la plupart des scientifiques mathématiciens, physiciens, chimistes, géologues, biologistes, ingénieurs, etc. connaissent ou ont côtoyé FORTRAN.

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.

Un peu d'histoire

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:

1978, la publication de la norme ANSI 77, qui donna le FORTRAN 77, encore largement utilisé aujourd'hui (les exemples ci-dessous sont écrits en 77). Le FORTRAN 77 introduit la programmation structurée dans les habitudes des programmeurs FORTRAN
1992, la publication de la norme ANSI 90, et pour la première fois la reconnaissance ISO (normalisation internationale). Le FORTRAN 90, outre quelques améliorations dans la gestion de la mémoire, amène essentiellement les outils nécessaires au calcul parallèle sur les machines multiprocesseurs. Il se décline en une nouvelle norme, qui n'est pas ANSI, et que l'on rencontre chez certains constructeurs, dont IBM, le HPF pour High Performance FORTRAN.
1997, la dernière norme en date, qui est une norme ISO: le FORTRAN 95.

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!).>

Les bases du FORTRAN

Organisation d'un programme FORTRAN

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

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)

END

type FUNCTION FunctionName (Arguments)

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)

END

type FUNCTION FunctionName (Arguments)

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 types de données

Les principaux types de donnée simples (nous verrons les tableaux plus loin) que nous utiliserons sont:

INTEGER : entier signé sur 4 octets (31 bits pour la valeur et un bit de signe). Nombre compris entre -2³¹ et 2³¹ -1
REAL : nombre "réel" codé en virgule flottante sur 4 octets. Nombre égal à +ou- 0.m*2^e , où la mantisse m est codée sur 23 bits et l'exposant e sur 8 bits signés. En gros, un REAL peut prendre en valeur absolue une valeur entre 1.4*10^-45 et 3.4*10³⁸
DOUBLE PRECISION : REAL dont la mantisse est codée sur 52 bits et l'exposant sur 11 bits. Un DOUBLE PRECISION peut prendre une valeur comprise entre 4.9*10^-324 et 1.8*10³⁰⁸
COMPLEX : 2 REAL dans un même nombre

Il en existe d'autres (CHARACTER,LOGICAL) que nous verrons plus loin.

Les constantes

Selon le type de la constante son écriture sera différente (attention, c'est un piège classique):

une constante INTEGER s'écrira simplement 1, -1, 3, 0
une constante REAL s'écrira obligatoirement avec un . ou le E de la notation en virgule flottante, par exemple: 0. , 1.24 , 31415E-4, 6.02E23
une constante DOUBLE PRECISION doit obligatoirement s'écrire en virgule flottante, avec un D en lieu et place du e : 31415D-4 , 0D0, 1.D0
une constante complex s'écrira sous forme d'un couple entre parenthèse en respectant les conventions des REAL: (0.0,1.0) , (1.1E-10, 3.5E-5)

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

Les variables

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:

Le nom d'une variable peut contenir des chiffres, des lettres, un blanc ou un underscore (_). Il doit commencer par une lettre. Le FORTRAN est indifférent aux minuscules/majuscules. J'ai pris l'habitude d'écrire mes noms de variables en minuscules pour les distinguer des mots clé. Ce n'est pas très conforme au standard FORTRAN des puristes, mais c'est personnel...
Seuls les 6 premiers caractères d'une variable sont significatifs. Attention au piège : MinimumX et MinimumY représentent la même variable! Bug assuré!!
Sauf indication contraire, en FORTRAN les variables commençant par i,j,k,l,m et n sont implicitement considérées comme des INTEGER. Les autres lettres désignent un REAL. Cela permet de ne pas déclarer explicitement les variables. Je vous déconseille formellement cette pratique. C'est la raison pour laquelle vous verrez en tête de tous mes programmes la déclaration IMPLICIT NONE, qui oblige la déclaration de chaque variable.
Chaque variable est locale par rapport à son bloc fonctionnel. Si vous déclarez un INTEGER toto dans le programme principal, il sera inconnu de toutes les subroutines et fonctions de ce programme. Pour partager des variables entre plusieurs blocs, on utilise l'instruction COMMON (voir son utilisation plus loin).

Les opérateurs et les fonctions

FORTRAN supporte les opérateurs arithmétiques suivants:

Addition +
Soustraction -
Multiplication *
Division /
Puissance **

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.

Les entrées/sorties

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.

Lire:

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).

Ecrire:

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é.

Les particularité d'un code source FORTRAN

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:

une ligne de commentaire doit comporter un C ou un * en colonne 1
les instructions doivent commencer en colonne 7 et leur longueur ne doit pas dépasser la colonne 72
les étiquettes (qu'on n'utilise plus tellement vu qu'on les trouve dans les GO TO, honnis, et les écritures formattées) doivent être écrites dans les colonnes 2 à 5
pour couper une instruction en plusieurs lignes, on positionne un caractère quelconque (je mets un !) dans la colonne 6 de chacune des lignes suivant la première ligne.

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.

Notre premier programme

Le programme VolumeSphere

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:

La saisie de la variable rayon. J'aurai pu écrire WRITE(*,'("Rayon de la sphere: ",$)'), ce qui aurait donné la même chose. Notez le $, qui empêche le retour à la ligne après l'affichage de la chaîne 'Rayon de la sphere : '.
Le calcul du volume: rien à dire tellement cela ressemble à la formule mathématique
Puis l'affichage, comme on l'a définit ci-dessus. En mode formaté, j'aurai pu écrire WRITE(*,'("Volume de la sphere = ",e15.8)') volume

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.

Installer un environnement de développement

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 2000). 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:

un environnement de développement (IDE) qui associe un éditeur, un compilateur et un linker. J'ai choisi VFort. VFort travaille en FORTRAN 77. Il comporte une aide et un guide de programmation. Il est téléchargeable ici
un logiciel de tracé de courbes, pour visualiser nos résultats. J'ai choisi GnuPlot, qui existe sous Windows et Linux. Il est téléchargeable sur le site www.gnuplot.info

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.

Saisir le programme

Lancer le programme VFort. La vue principale de l'IDE s'affiche.

La fenêtre de gauche liste les différents fichiers de votre projet.
La fenêtre de droite est notre fenêtre de travail. C'est là que nous saisirons notre code.
La fenêtre du bas est la console: c'est là que s'afficheront les différents messages de compilation (des messages d'erreur, souvent, mais aussi le résultat de la compilation) et de link (aussi des messages d'erreur...).

Pour saisir notre programme, il faut:

Créer un nouveau fichier (dans le menu FILE/New). Il porte un nom par défaut (Fort1, par exemple)
Saisir le texte du programme (ou faire un copier/coller). Attention au respect des colonnes : FORTRAN est intransigeant sur le sujet...
Sauvegarder le fichier dans votre répertoire (FILE/Save), en le nommant, VolumeSphere. On donne généralement le nom du programme au fichier source du programme principal. Le nom du fichier portera une extension .for
Créer un projet (Project/Create project) dans votre répertoire, en le nommant, VolumeSphere par exemple ou tout autre nom
Attacher le fichier que vous venez de créer (VolumeSphere dans notre exemple) au projet (Project/Add Files)
ATTENTION : sous Windows, pensez à ajouter à la variable d'environnement système PATH le chemin d'accès aux programmes de VFORT (le compilateur, le linker) sinon, vous auriez une erreur au build du projet. Dans mon cas, le chemin est C:\Program\VFort\Bin. Attention aussi en manipulant la variable PATH à ne pas effacer les chemins existants!

Le compiler - linker - exécuter

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...

Un peu plus compliqué

Dans cette partie un peu longue, on va examiner:

les conditions
les boucles
les tableaux
les fonctions et subroutines
les fichiers

Les conditions

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 :

x .EQ. y signifie x = y
x .NE. y signifie x <> y
x .GT. y signifie x > y
x .LT. y signifie x < y
x .GE. y signifie x >= y
x .LE. y signifie x <= y

Il définit aussi les opérateurs logiques unaires:

.OR. OU
.AND. ET
.XOR. OU exclusif
.NOT. NON

Munis de ces expressions logiques, on peut définir les instructions conditionnelles suivantes:

IF (condition) THEN

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

ELSE

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

ELSE IF (condition2) THEN

ELSE IF (condition3) THEN

ELSE

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

ELSE IF (x .EQ. 2) THEN

ELSE IF

(x .EQ. 3) THEN

ELSE

C cas de x différent de 1,2 et 3

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.

Les boucles

Sans doute, ce sont les instructions les plus utilisées en calcul. FORTRAN 77 comporte deux types de boucles.

La boucle DO...ENDDO

Sa syntaxe

DO var= debut, fin, pas

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:

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...

la boucle DO ... WHILE

Sa syntaxe

DO WHILE(condition)

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...

Les tableaux

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:

Premier point important, en FORTRAN, un tableau s'adresse en partant de 1, contrairement au C où l'on part de 0.
Deuxième point important: l'indice d'un tableau ne doit jamais dépasser la dimension fixée à sa déclaration. Sinon, plantage assuré.

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)

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

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 subroutines

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:

d'un calcul qui retourne une valeur simple à partir de variables. Par exemple, calculer le volume d'une sphère à partir de son rayon. Dans ce cas, on utilise une fonction. Sa principale caractéristique est de ne savoir retourner au programme qu'une seule valeur.
d'un traitement qui retourne plusieurs valeurs. Par exemple, on peut chercher les extrema et minima locaux dans un tableau. Dans ce cas on utilise une SUBROUTINE.

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).

Les fonctions

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

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é!

Les subroutines

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

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.

Les fichiers

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.

Créer et ouvrir un fichier séquentiel formaté

La syntaxe de l'instruction pour ouvrir un fichier est:

OPEN(unité,

FILE= chaine,

FORM = chaine,

STATUS = chaine

ERR = numlabel)

où

unité : numéro d'identification de l'unité du fichier. C'est un nombre entier compris entre 10 et 99, choisi arbitrairement par le programmeur. Il existe des conventions historiques, dont on ne s'encombrera pas ici. Chaque fois qu'une instruction voudra adresser le fichier ouvert, elle le désignera par ce numéro d'unité.
FILE : chaîne de caractères contenant le nom et le chemin du fichier, par exemple 'd:\NumLab\LabPhysique\data.txt'
FORM: chaîne de caractères: 'formatted' pour les fichiers formatés (valeur par défaut) ou 'unformatted' pour les fichiers binaires (dont nous ne parlerons pas)
STATUS : chaîne de caractères :
- 'new': le fichier n'existe pas, il est créé.S'il existe déjà le programme retourne une erreur
- 'old'le fichier existe déjà. S'il n'existe pas, le programme retourne une erreur
- 'unknow': ouvre le fichier, s'il n'existe pas il le créé, s'il existe il écrit à la fin de celui-ci. C'est le mode par défaut
ERR : numéro de label vers lequel le programme continuera son exécution s'il rencontre une erreur

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...

Lire un fichier séquentiel formaté

La syntaxe de l'instruction de lecture est:

READ(unité,

format,

ERR = numlabel,

END = numlabel) liste de données

où

unité : numéro d'identification de l'unité du fichier alloué à son ouverture
format : dans notre cas (fichier formaté), on indiquera *
ERR : numéro de label vers lequel le programme continuera son exécution s'il rencontre une erreur
END : numéro de label vers lequel le programme continuera son exécution lorsqu'il n'aura plus rien à lire dans le fichier (optionnel)

Exemple : READ (10,*,ERR=90, END =100) toto lit une donnée dans le fichier data.txt et la stocke dans la variable toto.

Ecrire une fichier séquentiel formaté

La syntaxe de l'instruction de lecture est:

WRITE(unité,

format,

END = numlabel) liste de données

où

unité : numéro d'identification de l'unité du fichier alloué à son ouverture
format : dans notre cas (fichier formaté), on indiquera *
END : numéro de label vers lequel le programme continuera son exécution lorsqu'il n'aura plus rien à lire dans le fichier (optionnel)

Exemple : WRITE (10,*) toto écrit la variable toto dans le fichier data.txt.

Fermer un fichier

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.

Exemple d'utilisation

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:

l'usage d'une variable CHARACTER et son assignation. On créé une variable contenant une chaîne de caractères, ici d'une longueur maximum de 80 caractères. On lui assigne simplement le nom du fichier que l'on écrit entre cotes.
la fonction REAL (ne pas confondre avec le type de variable!). La raison de l'usage de cette fonction, qui convertit un DOUBLE PRECISION en REAL est un peu curieux. Elle tient au logiciel de plot que vous utiliserez. certains, comme Matlab n'aime pas le 'D' utilisé par FORTRAN dans la notation double précision. Ils ne comprennent que le 'E' des REAL. Alors pour être sur de ne pas avoir d'ennuis, on convertit en REAL. Il risque d'y avoir quelques pertes de précision, mais pour faire de l'affichage graphique, cela a généralement peu d'importance.

Le premier programme de calcul - La méthode d'EULER

Résolution d'une équation différentielle par la méthode d'Euler

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 x_i et x_i+1 par un développement de Taylor d'ordre 1 (on dit que c'est une méthode à un pas):

y_i+1 = y_i + h*f(x_i, y_i) où h = x_i+1 - x_i

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.

Un exemple

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 :

d²(theta)/dt² = -(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:

theta_i+1 = theta_i + h*omega_i (1)

omega_i+1 = omega_i + h*a_{i =}omega_i + h*(-(g/l)sin(theta_i)) (2)

Nous retrouverons tout ce petit monde dans le programme ci-dessus:

le tableau X1 contient les valeurs d'angle calculées
le tableau X2 contient les valeurs de vitesse angulaire calculées
le tableau t contient les pas temporels (nous sommes, comme toujours en simulation, dans un monde de temps discret - A ne jamais oublier).
les conditions initiales sont : angle initial 60° (attention à l'exprimer en rd) et vitesse initiale nulle.

Le programme

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:

L'usage de COMMON entre le programme principal et les subroutines. Cette instruction permet de partager des variables entre deux ou plusieurs blocs fonctionnels (je vous rappelle que les variables en FORTRAN sont locales). Attention, COMMON ne permet pas de partager les variables déclarées en PARAMETER, ce qui explique que j'ai recopié dans la subroutine Derivee le PARAMETER g. Le nom du COMMON est quelconque. Vous pouvez nommer un COMMON comme vous le souhaitez, à condition que ce nom soit unique dans le programme!
Les commentaires en fin de ligne. Vous avez noté le signe ! qui permet d'écrire des commentaires en fin de ligne. Pratique pour commenter ses instructions...

Utiliser GnuPlot pour tracer ses résultats

GnuPlot : présentation du traceur de courbes

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.

Les principales commandes de GnuPlot

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..:

plot 'Euler.dat' : c'est la commande la plus simple. Je suppose que le répertoire de travail de gnuplot contient le fichier Euler.dat. Sinon, il faut le changer en utilisant la commande chdir du menu de Gnuplot. Vous obtenez un tracé avec les axes par défaut, sans légende, brut de fonderie.
set data style lines : pour modifier l'aspect de la courbe avant le tracé. On lui précise ici de relier les points par des segments, c'est plus joli! Un détail, cette option sera active tant que vous n'en changerait pas...
set xrange[0.0:2.0] ou set yrange[-8:8] pour fixer les bornes du tracé
set xlabel "Temps" pour indiquer la légende de l'axe des X. Evidement set ylabel "Angle" légende l'axe des Y.
set title "Variation de l'angle en fonction du temps" pour indiquer un titre à votre courbe...

Voilà pour le principal...

Comment tracer un fichier de données produit par un code FORTRAN

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!

Pour aller plus loin

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:

La programmation FORTRAN

Fortran 77 de Ficheux-Vapne
Numerical Recipes in Fortran: The Art of Scientific Computing de William H. Press, et al
Programmer en Fortran 90. Guide complet de Claude Delannoy
Fortran 90/95 for Scientists and Engineers de Stephen J. Chapman
Computer Applications in Physics With Fortran And Basic de Suresh Chandra