QUELQUES
FONCTIONS APL SIMPLES POUR
CLAUDE
CHACHATY
claude.chachaty@wanadoo.fr.
Cet article m'a été suggéré par la lecture de
l'excellent ouvrage de P. Hellings sur l'initiation aux calculs numériques en
astrophysique (1). Le premier chapitre
de cet ouvrage traite en grande partie des méthodes classiques de résolution
des équations et systèmes différentiels, également décrites avec leurs variantes
dans la réf. (2). Il m'a semblé utile de montrer que ces méthodes
permettant de traiter une multitude de problèmes parfois compliqués, peuvent se traduire en quelques
lignes d'APL comme on le
voit sur la fonction suivante :
EQUADIF
'Equation dy/dx = F(x,y) : ' ª DIMF„½FAB„FXY„
©
© - Substitution de x et y par A et B dans FAB:
©
FAB[(FAB='x')/¼DIMF]„'A' ª FAB[(FAB='y')/¼DIMF]„'B'
'Methodes: Euler [1], point-milieu [2]'
'prévision-correction [3], Runge-Kutta [4]' ª METH„Œ
'Valeurs initiales x0 et y0 : ' ª (x y)„(x0 y0)„Œ
'Nombre d''itérations et incrément dx: '
(NIT dx)„Œ ª X„Y„¼I„0
©
L1: I„I+1
–(METH=1)/'Y„Y,y„y+dx×–FXY ª x„x+dx'
–(METH=2)/'A„x+0.5×dx ª
B„y+0.5×dx×–FXY'
–(METH=2)/'Y„Y,y„y+dx×–FAB ª x„x+dx'
–(METH=3)/'A„x+dx ª
B„y+dx×fxy„–FXY'
–(METH=3)/'Y„Y,y„y+0.5×dx×fxy+–FAB ª x„A'
…(METH¬4)/L2 ª A„x+0.5×dx ª
K1„dx×–FXY
B„y+0.5×K1 ª K2„dx×–FAB ª B„y+0.5×K2 ª K3„dx×–FAB
x„A„x+dx ª B„y+K3 ª K4„dx×–FAB
Y„Y,y„y+(K1+K4+2×K2+K3)÷6
©
L2:X„X,x ª …(I<NIT)/L1
La méthode
d'Euler (1) est une simple itération
de la formule des accroissements finis. Les méthodes du point-milieu (2)
et de prévision-correction (3)
, représentent deux modes
d'interpolation des valeurs de f(x,y) obtenues selon la méthode 1. La méthode de Runge-Kutta (4) correspond à l'approximation de la
fonction f(x,y) par un développement limité
en série de Taylor. La précision de ces méthodes peut être évaluée par comparaison avec les
valeurs obtenues pour une équation différentielle dont la solution analytique
est connue; elle est en général dans l'ordre
1 << 2 ~ 3 <
4.
La méthode d'Euler, peu précise, n'est citée ici que pour mémoire. La
méthode 4, la plus précise, est deux
à trois fois plus lente que 2 et 3,
inconvénient mineur si l’on
considère les performances actuelles des microordinateurs.
La fonction SYSEQUA2, de structure
analogue à EQUADIF, s'applique à deux équations différentielles du premier
ordre couplées ou à une équation différentielle du deuxième ordre écrite sous la forme 
SYSEQUA2
'Equation dy/dx = F(x,y,z) :' ª DIMF„½AC„FXYZ„,
sysequa2 ª FABC„AC
''
'Equation dz/dx = G(x,y,z) :' ª DIMF„½AC„GXYZ„,
sysequa2 ª GABC„AC ©sysequa2 substitue A, B et C
à x, y et z dans FABC et GABC
'Méthodes : Point-milieu [1], prévision-correction [2], 'Runge-Kutta [3]'
METH„Œ
''
'Valeurs initiales x0, y0, z0 : '
X„Y„Z„¼I„0 ª (x y z)„(x0 y0 z0)„Œ
''
'Nombre d''itérations et incrément dx: '
(NIT dx)„Œ
L1: I„I+1
…(METH¬1)/L2 ª A„x+0.5×dx ª
B„y+0.5×dx×–FXYZ
C„z+0.5×dx×–GXYZ ª Y„Y,y„y+dx×–FABC
Z„Z,z„z+dx×–GABC ª x„x+dx
©
L2:…(METH¬2)/L3 ª 'A„x+dx ª B„y+dx×fxyz„–FXYZ
C„z+dx×gxyz„–GXYZ ª Y„Y,y„y+0.5×dx×fxyz+–FABC
Z„Z,z„z+0.5×dx×gxyz+–GABC ª x„A
©
L3:…(METH¬3)/L4 ª
A„x+0.5×dx ª K1„dx×–FXYZ
B„y+0.5×K1 ª M1„dx×–GXYZ
C„z+0.5×M1 ª K2„dx×–FABC ª
M2„dx×–GABC
B„y+0.5×K2 ª C„z+0.5×M2 ª K3„dx×–FABC
M3„dx×–GABC ª x„A„x+dx ª B„y+K3
C„z+M3 ª K4„dx×–FABC ª M4„dx×–GABC
Y„Y,y„y+(K1+K4+2×K2+K3)÷6
Z„Z,z„z+(M1+M4+2×M2+M3)÷6
©
L4:X„X,x ª …(I<NIT)/L1
.
Le même
procédé de calcul peut être
facilement étendu à un plus grand nombre
d'équations.
EXEMPLES D'APPLICATIONS
Evolution thermique d'une réaction en chaîne.
La cinétique de la plupart des réactions en
chaîne comme la combustion, la décomposition thermique, la polymérisation, dépend d'une manière
critique de la différence entre les
vitesses de dégagement et d'évacuation de la chaleur produite en milieu
confiné. Ce problème a fait l'objet d'un
nombre considérable de travaux dont l'un
des plus importants est la réf. (3).
Prenons l'exemple d'une polymérisation où
la terminaison des chaînes ne s’effectue que par consommation du monomère. La
vitesse de formation du polymère est :
[1]
R0 :
concentration du catalyseur, C0 : concentration initiale de
monomère.
Ap et Ep :
facteur préexponentiel et énergie d'activation de propagation.
L'équation d’évolution de la température
est la somme de deux termes correspondant respectivement au dégagement et à
l'évacuation de la chaleur produite par la réaction :
[2]
Qp, cp
et : enthalpie de réaction, chaleur spécifique et densité de la solution.
Te : température externe du réacteur de
surface S et de volume V.
: coefficient d'échange thermique.
La figure 1 représente l'évolution de la
température calculée pour un réacteur
contenant une solution de styrène en cours de polymérisation. Le système
d'équations [1, 2] qui n'a pas de
solution analytique, est résolu avec la fonction SYSEQUA2 par la méthode de
Runge-Kutta. La concentration initiale
des réactifs, le rapport surface/volume du réacteur et la température externe
de celui-ci ont une une influence cruciale sur le régime thermique de la
réaction. Cette figure montre par
exemple que l' augmentation d'une dizaine de degré de la température externe
provoque la transition d'un régime thermique quasi-stationnaire ou à variation
lente, à une régime transitoire qui se manifeste par une brusque élévation de
la température interne suivie d’une décroissance rapide due à la consommation
complète du monomère.
Les équations [1] et [2] sont d'un intérêt tout à fait général; elles s'appliquent avec quelques variantes à une réaction en chaîne comme la polymérisation, à l'inflammation d'une allumette, à la détonation d'un explosif , à la surchauffe due à la fermentation des graines dans un silo etc...
|
|
|
Figure
1. Différence entre les températures
interne et externe du réacteur.
Paramètres
: |
,
, 
s
Trajectoire d'un objet dans un champ gravitationnel.
L'exemple suivant concerne le problème à
trois corps dont une version simplifiée est donnée dans la réf. (1). On
considère deux astres en interaction gravitationnelle. L'astre le moins massif
(II) est supposé parcourir une orbite quasi-circulaire à la vitesse angulaire 
autour de l'astre I .
La trajectoire du troisième corps (III) dont la masse est négligeable devant
celle de I et II, est exprimée dans un référentiel tournant à la vitesse autour du centre de gravité du système I-II pris comme origine.
L'unité de temps est la période T, l'unité de longueur est la distance
I-II et l'abscisse de l'astre II est
choisie négative.
La trajectoire de l'objet III est obtenue par
résolution de deux équations différentielles du second ordre couplées, que l'on
peut remplacer par le système de quatre équations du premier ordre :
[3a]
[3b]
[3c]
sont les masses
réduites des astres I et II 
dont les coordonnées
sont 
. Les deux premiers termes des équations 3b et 3c
correspondent à la loi de gravitation de Newton, les autres sont introduits
pour tenir compte de la rotation du référentiel XY (voir réf. (1)).
Le système [3a,b,c] a été résolu par la
méthode de Runge-Kutta au moyen de la fonction SYSEQUA4, analogue à SYSEQUA2,
valable pour quatre équations du premier ordre. Les trajectoires représentées
sur la figure 2 ont été calculées sur 1200 points avec un pas de 0.01 T,
soit douze révolutions de
La fonction NEWTON de Chaitin (4) permet la
résolution numérique du problème à N corps dans toute sa généralité et ne
nécessite donc pas les simplifications
mentionnées plus haut. Que l'on utilise SYSEQUA4 ou NEWTON, on constate que les
formes des trajectoires sont extraordinairement dépendantes des conditions
initiales et donc difficilement prévisibles.
|
|
|
|
|
|
Figure 2. Trajectoire d'un objet de faible masse dans le système Terre (m1=0.98786, x1= 0.01214) - Lune (m2 = 0.01214,
x2= - 0.98786) .
Positions initiales : y0 = 0,
x0 variable. Vitesses initiales : u0 =
0, v0 = - 1. Durée totale 12
T.
Les fonctions APL décrites dans cet
article, comme la plupart de celles que j'utilise, existent en APL*PLUS II
version 5.1 pour DOS et en deux applications Windows: APL*PLUS III version
1.2 et APL2win (IBM APL2 beta level
7). Les trajectoires de la figure 2 par
exemple, ont été calculées avec ces logiciels
en 10.9, 9.3 et 11.3 s respectivement, sur un PC Compaq 133 MHz avec
processeur Pentium. Ces logiciels ont donc des performances à peu près
équivalentes dans le cas d' une fonction opérant en boucle avec un grand nombre
d'itérations. Pour ce qui est des fonctions
ne comportant que des opérations sur des vecteurs et matrices réels, APL2win est environ deux fois
plus rapide que les APL*PLUS II et III. Dans les opérations sur les complexes,
en particulier pour les inversions de matrices, les vitesses relatives de
calcul sont en moyenne dans l'ordre 3, 1.5 et 1 respectivement pour APL2win ,
APL*PLUS III et
APL*PLUS II. Les opérations sur
les variables complexes sont plus simples et plus rapides en APL2 IBM qu'en
APL*PLUS du fait qu'elles sont considérées comme de même rang que les variables réelles équivalentes et
font appel aux mêmes opérateurs, à l’exception du domino, remplacé par une
fonction. Ainsi un vecteur à N composantes complexes est un vecteur de dimension N en APL2 IBM et un tableau de dimension 2,N en
APL*PLUS; de même le produit scalaire de deux vecteurs complexes s'écrit V+.W en APL2 IBM et V CXMATTIMES W en APL*PLUS etc...
Un autre avantage d'APL2win
sur les deux autres APL dont il est question ici, est l’existence de
la zone (alias workspace) GRAPHPAK qui
contient diverses fonctions d'utilisation facile, pour graphiques bi et
tridimensionnels. Pour une impression de bonne qualité des graphiques, il est
cependant nécessaire de transférer les tableaux numériques d'APLwin à EXCEL.
Pour le moment je n'ai réussi ce transfert que par la séquence APL2win fichier ASCII sous DOS reformater en APL*PLUS III EXCEL. Je n'ai pas trouvé non plus en APL2win, l'équivalent des
commandes )CMD ou CMD qui en APL*PLUS permettent de communiquer rapidement avec DOS en
cours de session APL.
En APL*PLUS II seulement, CMD [instruction DOS] dans une fonction APL, permet d’exécuter un
fichier *.exe et de récupérer automatiquement les tableaux calculés. APL*PLUS II a également l'avantage de
permettre de disposer de zones plus
grandes qu'en APL sous Windows. Ainsi
pour 16 Mo de RAM, on dispose de zones de 14 Mo en APL*PLUS II alors qu'en
APL*PLUS III et APL2win , il est
préférable de limiter leur taille à 5 Mo pour
ne pas compromettre la rapidité des calculs.
Il est donc avantageux d'utiliser APL2win pour les calculs longs que l'on rencontre souvent dans les
problèmes d'optimisation, APL*PLUS III si l'on souhaite l'édition immédiate et
de bonne qualité de données numériques et de graphiques sous EXCEL, APL*PLUS II
pour des fonctions générant des variables de grande taille ou faisant appel à
des fichiers DOS en cours d'exécution.
Références
1 - P. Hellings, Astrophysics with
a PC, an Introduction to Computational Astrophysics. Willmann-Bell,
Richmond, Virginia (USA) 1994. ISBN 943396-43-3.
2 - W.H. Press, S. A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran, the Art of
Scientific Computing, seconde édition, chap. 16, Cambridge University
Press, Cambridge 1992. ISBN 0 - 521 -
43064 - X.
3 - D.A. Frank-Kamenetzky, Diffusion
and Heat Exchange in Chemical Kinetics, chap. 6 et 7, Princeton University
Press, Princeton 1955.
4 - G.J. Chaitin, An APL2 Gallery
of Mathematical Physics. A course Outline. Proceedings Japan 85 APL
Symposium, pp 1-56, IBM Japan 1985.