mise à jour : 31-03-2010

cinet.chim

complément
Intégration numérique

logo Sa

complément

Intégration numérique

L'intégration numérique est un très vaste et complexe sujet. Notre but n'est pas de le traiter, ni de façon exhaustive, ni dans le détail, mais, aussi simplement que possible, d'en dégager les principes utiles pour comprendre comment les intégrateurs fonctionnent et comment les utiliser correctement.

Nous n'évoquerons ici que l'intégration numérique des équations différentielles avec conditions initiales, les problèmes avec conditions aux limites, plus difficiles, ne concernant guère la cinétique chimique.

Vous pourrez trouver un exposé très complet, orienté vers la programmation, et relativement abordable dans Numerical Recipes. Les trois articles suivants sont plus techniques et concernent directement l'algorithme SRK utilisé dans Sa.

Bibliographie

- Numerical Recipes in C.  Chapitre 16 : Integration of Ordinary Differential Equations http://www.nr.com

- Generalized Runge-Kutta Methods of Order Four with Stepsize Control for Stiff Ordinary Differential Equations. P. Kaps, P. Rentrop ; Numerische Mathematik, 33, 55-68, 1979

- A Reliable Rosenbrock Integrator for Stiff Differential Equations.  B. A. Gottwald, G. Wanner ; Computing, 26, 355-360, 1981

- Application of a variable-order semi-implicit Runge-Kutta method to chemical models.  P. Kaps, P. Rentrop ; Computers & Chemical Engineering, 8, 393-396, 1984


1. Le problème à résoudre

Etant donné un ensemble de N équations différentielles du premier ordre couplées, de la forme générale

dyi(x)/dx  =  fi (x, y1, ..., yN)      i = 1 à N     (1)

où les fonctions fi sont connues,

et des conditions initiales connues

x0 , y1,0 , ... , yN,0

on veut trouver numériquement les valeurs de

yi, j  =  G(xj)

pour une plages de valeurs de x : x0, ..., xj, ..., xM.

L'équation (1) peut s'écrire aussi sous la forme vectorielle :

y' = f(x, y)

Les équations différentielles d'ordre supérieur peuvent toujours se ramener à un ensemble d'équations du premier ordre.


2. La démarche de base

Pour simplifier nous ne considérerons d'abord que le problème à une seule variable (N = 1). L'équation (1) se réduit alors à :

dy(x)/dx  =  f (x, y)  

L'idée de base est de calculer une valeur correspondant à xn+1 = xn+ h à partir des valeurs connues pour xn, en commençant par x0. Cette valeur devient alors la valeur initiale pour un nouveau pas, et de proche en proche on construit la table des valeurs désirée.

L'équation différentielle écrite sous forme discrète devient :

Δy  =  Δx f (x, y) = h f (x, y)

équivalente à (1) si h → 0, et on peut donc calculer le point suivant par la formule d'Euler :

yn+1  =  yn + h f (xn, yn)     (2)

qui est en fait le début d'un développement de Taylor :

yn+1  =  yn + h f (xn, yn)+ O(h2)     (3)

méthode d'Euler

Une telle méthode est dite du premier ordre, puisqu'elle néglige les termes d'ordre supérieur. L'erreur de troncature (due à la méthode) est proportionnelle à h2.

D'autre part, elle extrapole la valeur de la dérivée au point initial, Mn, à tout l'intervalle à parcourir. Le nouveau point obtenu, M'n+1, est donc différent du point théorique attendu, Mn+1. Le saut appliqué à y est :

k1 = h f (xn,yn)

Ainsi, même si le pas d'intégration h est très petit, ce qui suppose un temps de calcul très important, les erreurs commises à chaque pas, aussi petites soient elles, non seulement s'accumulent mais ont tendance à s'amplifier : l'intégration diverge.

Cette méthode, insatisfaisante telle quelle, est cependant à la base de toutes les méthodes plus sophistiquées, et montre que le pas d'intégration reste le paramètre fondamental de l'intégration numérique.

3. Méthode de Runge-Kutta

Pour remédier, du moins en partie, aux défauts de la méthode d'Euler, l'idée de la méthode de Runge-Kutta est de remplacer l'évaluation de la dérivée au point Mn par une moyenne des dérivées sur l'intervalle h.

Pour cela, par exemple, évaluons la dérivée au point milieu (flèche rouge) de l'intervalle précédent. Cela nous donne un nouveau terme de saut pour y :

k2 = h f (xn + h/2, yn + k1/2)

Runge-Kutta deuxième ordre

qu'on applique à yn. L'équation (2) devient alors :

yn+1  =  yn + k2     (4)

Et on peut démontrer que cette opération annule l'erreur du premier ordre, de sorte qu'on peut écrire plus rigoureusement :

yn+1  =  yn + k2 + O(h3)     (5)

d'où le nom de Runge-Kutta du second ordre.

Le calcul de la moyenne des dérivées sur l'intervalle h peut être encore amélioré en calculant celles-ci sur quatre points, comme indiqué sur la figure ci-dessous.

Runge-Kutta quatrième ordre

Les valeurs à moyenner sont alors :

k1 = h f (xn,yn)

k2 = h f (xn+h/2, yn+k1/2)

k3 = h f (xn+h/2, yn+k2/2)

k4 = h f (xn+h, yn+k3)

et on obtient une approximation de yn+1 :

yn+1  =  yn + k1/6 + k2/3 + k3/3 + k4/6 + O(h5)     (6)

Cette méthode, dite de Runge-Kutta du quatrième ordre (RK4), est une des plus utilisées dans des problèmes relativement simples, se traduisant par des courbes "douces". Elle ne convient pas, par contre, pour des problèmes présentant de très fortes variations des dérivées, comme l'autocatalyse chimique en est un exemple.

Ses principaux avantages sont la facilité de programmation et de modification du pas d'intégration h, et sa stabilité (en général).

La question de sa précision est complexe, car elle n'est pas exclusivement liée à l'ordre de la méthode, et hors du champ de cet exposé. Théoriquement, il est possible d'améliorer celle-ci en diminuant le pas d'intégration. Mais c'est au prix d'un temps de calcul qui devient rapidement prohibitif. En effet, pour exécuter un pas d'intégration, il faut calculer 4 fois la dérivée. Si ce calcul est peu coûteux, ce n'est guère un problème, mais ça le devient dans le cas, par exemple, où il y a beaucoup de variables (voir l'équation du problème général, occultée dans cette présentation schématique)

4. Pas d'intégration auto-adaptatif

Il est clair que si l'équation à intégrer résultait en une droite, un seul pas d'intégration couvrant tout le domaine serait suffisant, et assurerait une précision limitée seulement par les erreurs d'arrondi du calculateur. A l'inverse, plus la courbure augmente, plus le pas d'intégration doit être choisi petit pour assurer un minimum de précision. Or, la méthode précédente imposerait d'appliquer ce pas minimum tout au long de l'intégration, y compris dans des régions quasi-linéaires.

D'autre part, il lui manque un test interne qui permettrait d'affirmer que tel pas d'intégration a été correctement traité, ou non, selon des critères objectifs donnés. Devant un premier résultat négatif de ce test, le pas d'intégration pourrait être automatiquement réduit, jusqu'à ce que le test redevienne positif. Tandis que devant un premier résultat positif, le pas pourrait être automatiquement augmenté tant que celui-ci reste positif.

Bien entendu, la mise en œuvre d'une telle procédure aura également un coût en temps de calcul, qu'il faudra minimiser de façon à être globalement gagnant.

Une technique simple consiste à effectuer 2 pas d'intégration de 2 manières différentes :

- une fois avec un pas égal à 2h ; résultat : y1

- une deuxième fois en faisant 2 pas égaux à h ; résultat : y2

Nous admettons que la différence

Δ  =  y2y1

est une estimation de l'erreur de troncature commise.

Supposons maintenant que nous considérions comme acceptable une erreur égale ou inférieure à Δ0 (valeur qui dépend du paramètre habituellement appelé tolérance de l'intégrateur), alors, on peut démontrer qu'une formule du type

h'  =  h | Δ0 / Δ | α

α est un exposant positif (par exemple α = 0.2),

permet d'obtenir une valeur du pas h' qui donnera directement une erreur de l'ordre de Δ0. Ainsi, si Δ < Δ0 (pas d'intégration réussi), le pas pourra être augmenté (h' > h) pour l'étape suivante. Au contraire, si Δ > Δ0 , l'intégration est erronée et doit être recommencée avec un pas plus petit (h' < h).

Nous nous en tenons ici au principe du pas auto-adaptatif. Il existe d'autres façons de faire, d'autres estimations de l'erreur de troncature et d'autres formules de correction du pas.

Noter que Δ et Δ0 sont en réalité des grandeurs vectorielles : la tolérance doit évidemment s'appliquer à toutes les variables du système à intégrer. Ainsi, au cours d'une intégration d'un système multivariables, c'est la variable la plus "fautive" qui imposera le pas d'intégration.

5. Algorithmes explicite, implicite et semi-implicite

Pour illustrer ces notions, revenons à la formule d'Euler :

yn+1  =  yn + h f (xn, yn)     (2)

cette formule est explicite parce qu'elle donne directement yn+1 en fonction de la valeur connue yn. Et nous avons vu que sa stabilité n'est pas garantie : l'intégration effectuée de la sorte, ou d'une façon analogue, peut diverger. Auquel cas, même l'état final du système (t ∞ pour une cinétique) peut s'avérer complètement faux.

Si nous modifions la formule d'Euler en évaluant la dérivée non en xn,yn mais en xn+1,yn+1, nous obtenons la nouvelle formule

yn+1  =  yn + h f (xn+1, yn+1)     (7) 

qui est implicite, car faisant appel au calcul de la dérivée au point qu'on se propose justement de calculer.

Supposons maintenant que l'équation (ou les équations) à intégrer soit linéaire, c'est-à-dire de la forme (qui doit vous rappeler quelque chose) :

dy/dx = −ky     (k étant une constante > 0)

l'équation (7) s'écrit alors

yn+1  =  yn + h(−kyn+1)

qui peut être rendue explicite (parce que l'équation différentielle est linéaire) :

yn+1  =  yn / (1 + hk)     (8)

Or l'application de cette formule, indépendamment de sa précision, est stable et converge vers l'état final correct.

Ce résultat se généralise bien sûr dans le cas d'un système d'équations linéaires. Le coût supplémentaire pour utiliser une méthode implicite est alors l'inversion d'une matrice à chaque pas.

Si les équations ne sont pas linéaires (c'est alors qu'on a le plus besoin d'un algorithme stable), les choses sont un peu plus compliquées. On doit procéder à une linéarisation locale. Nous n'entrerons pas dans le détail de cette procédure et ne mentionnerons que le fait qu'elle s'obtient également par l'inversion d'une matrice construite à partir de la matrice jacobienne du système. Ainsi, partant de l'équation générale (1), vectorielle, la formule d'intégration implicite analogue à (8) devient :

yn+1  =  yn + h [ 1hf/∂y|yn ]−1. f(yn)     (9)

y et f désignent des vecteurs, 1 la matrice unité et ∂f/∂y la matrice jacobienne.

La matrice jacobienne d'une fonction vectorielle f d'un vecteur y est la matrice de ses dérivées partielles par rapport aux composantes de ce vecteur.

Cette résolution de l'équation implicite par linéarisation est appelée méthode semi-implicite. La formule (9) correspond donc à la méthode d'Euler semi-implicite. Il n'est pas totalement garanti qu'elle soit stable, mais elle l'est en général car localement le comportement est en général semblable à celui des systèmes linéaires.

L'algorithme SRK (Semi-implicite Runge-Kutta) ou méthode de Kaps-Rentrop, combine les méthodes de Runge-Kutta du quatrième ordre, du pas d'intégration auto-adaptatif et de la linéarisation semi-implicite. Il allie ainsi une excellente stabilité et précision, ainsi qu'un temps de calcul très réduit. Il est de plus relativement léger à programmer (6 ou 7 pages en langage C) et simple à paramétrer.

Il ne faut pas oublier cependant qu'on peut toujours mettre un intégrateur numérique en défaut... et le corollaire : il faut toujours contrôler ce qu'on fait !

6. L'intégrateur de Sa

L'intégrateur de Sa, srkvi, est un intégrateur SRK. Quelques améliorations mineures, mais très utiles, lui ont été apportées.

Tout d'abord, dans Sa, l'intégrateur est le plus souvent utilisé couplé à l'optimisateur. Or les données expérimentales, cinétiques pour ce qui nous concerne, ne se présentent pas toujours comme une table de valeurs prises à intervalle régulier : d'une part il peut y avoir des données absentes, parce que la mesure n'a pas fonctionné correctement, d'autre part l'incrémentation de la variable indépendante n'est pas forcément constante, elle peut être logarithmique par exemple. Dans ces conditions, ce qu'on attend de l'intégrateur, c'est qu'il nous donne les valeurs correspondant le plus exactement possible aux points expérimentaux, de façon à pouvoir comparer données calculées et données expérimentales et effectuer l'optimisation des paramètres du modèle.

La version srkvi (SRK Variable Incrément) permet de passer à l'intégrateur une table quelconque de valeurs de la variable indépendante (la variable ind[ ]). Cette table est générée par Sa en fonction du problème à traiter. Et l'intégrateur interprète en interne chaque intervalle de cette table comme un sous-domaine à intégrer.

Cette façon de faire a permis en outre d'ajouter la possibilité d'ajuster de façon automatique le pas maximum d'intégration à l'incrément en cours de la variable indépendante (mode Auto).

Il peut être utile de préciser le vocabulaire employé :

Le pas d'intégration est le pas réellement utilisé en interne par l'intégrateur. Le pas maximum d'intégration est la valeur maximum de ce dernier autorisée par l'utilisateur (ou Sa en mode Auto). L'incrément de la variable indépendante est l'intervalle entre deux valeurs successives pour lesquelles on attend que l'intégrateur fournisse un résultat (cela est quelquefois appelé le pas de sortie de l'intégrateur).

Le mode Auto simplifie le choix du paramètre pas max et améliore la rapidité du calcul dans la plupart des cas. Toutefois, il n'est pas indiqué dans les problèmes difficiles.

Le test d'erreur de troncature de la version originale de SRK pouvait poser des problèmes lorsque certaines variables devenaient très petite, car il comporte une division par ces variables. Pour y remédier, nous avons introduit un nouveau paramètre, c_min, qui se substitue aux variables qui lui deviennent inférieures (uniquement au niveau du test bien entendu). Le résultat est une atténuation de l'erreur de troncature commise sur ces variables. Cela ne perturbe en général pas l'intégration... si les variables en question ne jouent pas un rôle crucial dans la dynamique du système, ce qu'il est parfois difficile de savoir à priori. Par exemple, la valeur c_min = 10−8 que nous préconisons pour la plupart des problèmes de cinétique en solution, aux concentrations courantes, a due être réduite à 10−15 pour traiter des réactions radicalaires où des radicaux demeuraient à des concentrations extrêmement faibles mais jouaient un rôle important.

L'appel de srkvi se fait par l'instruction suivante :

srkvi(n_diff, &ca[first_var], ind, npt, h0, tol, iset, jacob, h_compt, c_min);

Les quatre premiers arguments concernent la communication avec les données du problème :

n_diff : nombre d'équations différentielles

first_var : numéro de la première variable à intégrer

ind : table des valeurs de la variable indépendante pour lesquelles on attend un résultat

npt : nombre de ces valeurs

A ces quatre arguments, s'en ajoute un cinquième qui demeure masqué pour l'utilisation courante :

eq_system : c'est le nom générique de la fonction qui calcule les dérivées. Par défaut, c'est la fonction eqdiff que vous avez déjà utilisé, mais il est possible de lui affecter le nom d'une autre fonction. Cela permet par exemple d'appeler successivement plusieurs fois l'intégrateur avec des systèmes différentiels différents, ou même de changer de système différentiel au cours de l'intégration. Voir le Manuel de Sa, chapitre XII, page 71 et les Exemples commentés, exemple XI, page 83.

Les six derniers arguments sont les paramètres d'intégration :

h0 : pas maximum d'intégration autorisé (son effet est annulé dans le mode Auto et le pas maximum est alors égal à l'incrément en cours).

Il est crucial de régler le pas maximum d'intégration en accord avec l'échelle de temps du problème traité. S'il est évident qu'il ne doit pas être trop grand et en tout cas ne pas dépasser l'incrément de la variable indépendante, il ne doit pas non plus être trop petit. En effet, outre un temps de calcul inutilement allongé, un pas trop petit peut provoquer des problèmes, notament au démarrage de l'intégration, alors que le pas interne de l'intégrateur ne s'est pas encore adapté.  Cela peut induire des essais d'intégration vers les temps négatifs, et se traduire par une modification des conditions initiales !

tol : tolérance du test d'erreur de troncature, destiné à régler le pas réel d'intégration.

iset : peut être égal seulement à 1 ou 2 et détermine lequel des deux jeux de paramètres internes de srkvi doit être utilisé. Cela n'a le plus souvent aucun effet, mais dans certaines situations difficiles, il peut arriver que le changement de jeu suffise à surmonter un problème d'intégration.

jacob : égal à 1 par défaut, la matrice jacobienne est alors calculée numériquement par une méthode de différences finies. Si l'on sait la calculer analytiquement, on peut mettre ce paramètre à zéro. Il faut alors fournir la fonction correspondante. Nous conseillons de s'en tenir au mode par défaut.

h_compt : paramètre obsolète, conservé uniquement pour des raisons de compatibilité avec des travaux antérieurs aux améliorations de srkvi.

c_min : voir ci-dessus.

D'autre part, dans la version srksta, nous avons adapté l'intégrateur srkvi pour obtenir directement l'état stationnaire d'un système. Nous avons donc ajouté un test interne permettant de décider si un état stationnaire est atteint, et l'intégration se poursuit automatiquement jusque là. Seul le résultat correspondant à l'état stationnaire est délivré. Cette version est très utile pour étudier la variation de l'état stationnaire d'un système en fonction d'un paramètre donné. Par exemple, on peut étudier ainsi l'évolution de l'état stationnaire d'un réacteur CSTR en fonction du débit. Une autre application plutôt inattendue est son utilisation pour résoudre les systèmes d'équations algébriques non linéaires, par exemple dans les problèmes d'équilibres chimiques.

L'appel de srksta est analogue à celui de srkvi :

srksta(n_diff, ca1, ind1, 0, 0, 0, h0, tol, iset, jacob, h_compt, c_min);

où l'on reconnaît essentiellement les mêmes arguments. La seule différence concerne les deuxième et troisième arguments, qui doivent être préalablement déclarés et initialisés. Voir les Templates de l'éditeur Sa et les Exemples commentés, exemple 3 page 15, et exemple 9.2 page 72.

Signalons pour terminer deux petits outils de Sa parfois utiles :

- le transfert des valeurs finales d'un calcul dans les valeurs initiales (onglet Variables, bouton portant une flèche bleue). Cette possibilité n'est pas spécifique à l'intégration numérique, mais permet, par exemple, de poursuivre celle-ci si l'état stationnaire n'est pas atteint à la fin de l'intervalle initialement choisi. Cela peut s'avérer très utile si le temps de calcul est très long.

- un outil de suivi de l'intégration numérique (fenêtre principale de Sa), permettant d'afficher les résultats pendant l'intégration, d'intégrer "pas à pas", ou de mettre un point d'arrêt. Sans être un "débogueur", cela peut être très utile pour détecter les causes d'une intégration défectueuse.