FR 2947071 A1 20101224 FR 0954120 A 20090618 -1- � Procédé et système de traitement de données numériques pour la simulation numérique ou la résolution numérique d'un système physique à équations différentielles � La présente invention concerne un procédé de traitement de données numériques pour le calcul de la valeur d'une ou plusieurs grandeurs physiques, par exemple les valeurs de la température dans différents points de l'épaisseur d'un mur soumis à deux environnements de températures différentes sur ses deux faces. Ce procédé s'applique dans le cas où ces grandeurs physiques, dans cet exemple la température, dépendent d'une loi de fonctionnement représentée sous la forme d'équations différentielles du second ordre dont on obtient ainsi une résolution numérique. Domaine technique L'ordinateur est très utilisé pour obtenir par le calcul une prédiction ou une simulation numérique représentant le fonctionnement et l'état d'un système ou d'un environnement, sans avoir besoin de réaliser concrètement ce système ou cet environnement. De tels calculs sont utilisés par exemple pour prédire les contraintes dans une construction ou une pièce mécanique, et ainsi dimensionner ses éléments avant la fabrication. Ils sont utilisés aussi pour prédire le fonctionnement et les performances, par exemple pour réaliser des essais � simulés � d'un système thermique, ou de résistance à l'avancement dans l'air pour un véhicule réalisant ainsi une � soufflerie numérique �. Dans de tels cas, on utilise souvent des méthodes de résolution dites � numériques �, qui consistent à utiliser l'informatique pour trouver un résultat numérique, en reproduisant par le calcul le comportement d'un système à partir de valeurs initiales choisies de façon possiblement empirique, et comprenant plus particulièrement les paramètres de modélisation ainsi que les conditions de frontières (dites aussi conditions aux limites). On pourrait dire qu'il s'agit de partir des causes pour calculer le résultat, plutôt que de partir du résultat souhaité pour calculer les besoins. -2- Par exemple lorsque les équations du fonctionnement d'un système sont trop complexes à calculer, ou que l'on ne connaît pas de méthode mathématique permettant leur résolution analytique, on utilise une méthode de résolution numérique pour trouver des solutions numériques pour ces équations. Pour cela, on définit un maillage qui divise le système étudié en un certain nombre de petits éléments ou � mailles �. On utilise alors l'ordinateur pour obtenir des valeurs numériques au niveau de ces différentes mailles, et tester par le calcul un grand nombre de choix possibles pour ces valeurs initiales, jusqu'à obtenir une valeur plus ou moins approximative de la solution ou des solutions de l'équation. Selon la nature du système à modéliser, le maillage peut aller de quelques dizaines de points à plusieurs millions ou plus, en particulier dans les systèmes non-finis tels que par exemple des modélisations acoustiques en extérieur. Dans de tels calculs, on cherche à minimiser les erreurs de prédiction ou simulation tout en limitant le volume de calculs à réaliser. Ce volume de calcul est un point particulièrement critique dans des domaines de plus en plus nombreux, qui vont par exemple de la simulation numérique d'un écoulement de gaz de combustion d'un moteur jusqu'au calcul de prévisions météorologiques à grande échelle. L'invention s'applique plus particulièrement aux procédés et systèmes de résolution numérique appliquée à des systèmes ou des environnements dont le fonctionnement peut être représenté par des équations différentielles ou intégrales, souvent en dérivées ou intégrales secondes. Par exemple, on sait de nos jours simuler numériquement les écoulements, couplés ou non, de fluide, de chaleur, de masse ..., en résolvant numériquement, sur ordinateur, les équations non linéaires aux dérivées partielles qui régissent ces écoulements. Ces équations sont de type Navier-Stokes. Etat de la technique Il existe de nombreux programmes informatiques industriels, ou � codes � au sens d'un code programme, dédiés à des applications spécifiques (aéronautique, élaboration des matériaux (métallurgie, croissance cristalline, plasturgie, thermique des systèmes électroniques, -3- etc.) qui permettent, à moindre coût comparé aux coûts expérimentaux, d'explorer et de dimensionner des configurations de réalisation de machines, ou des procédés industriels. De tels logiciels � industriels � sont vendus par exemple par les sociétés Flomerics ou Dassault Systèmes pour le domaine des fluides, ou les sociétés Progistik ou Winbeam en résistance des matériaux, ou les sociétés Trim ou Silvaco pour le domaine du transport ionique. Le développement ou l'acquisition de ces programmes industriels est extrêmement complexe et onéreux. De plus, les solutions qu'ils fournissent sont approximatives, alors qu'elles résultent d'une dynamique non linéaire, dont on sait bien aujourd'hui qu'elle est sensible aux conditions initiales, et donc aux erreurs avec lesquelles se calculent ces solutions. Par exemple, l'identification correctement prédictive de bifurcations entre régimes d'écoulements (par exemple le passage d'un écoulement laminaire à un écoulement turbulent) ne peut que bénéficier d'une augmentation de la précision numérique. La résolution numérique de telles équations peut se faire actuellement par des méthodes dites � d'approximation locale �, telles que les méthodes par Différences Finies (� DF �), ou les méthodes par Eléments Finis (� EF �) ou par Volumes Finis (� VF �). Ces méthodes consistent, pour obtenir une valeur en un point spécifique du maillage, à utiliser quelques points voisins pour y évaluer numériquement les valeurs de la dérivée seconde. Si l'on utilise les deux points immédiatement voisins de chaque côté du point concerné, on parle d'un � schéma 3 points �. La précision de ces méthodes augmente lorsque l'on augmente le nombre de points, mais la quantité de calcul aussi. Pour l'ensemble des éléments du maillage, les différentes équations sont représentées sous la forme d'une matrice comprenant des valeurs numériques correspondant aux coefficients des inconnues dans chacune de ces équations. Concrètement, chaque calcul d'un point comprend un travail numérique de l'ordinateur sur toutes les valeurs de cette matrice, pour calculer sa matrice � inverse �. Dans ces méthodes, tous les schémas connus à ce jour (Différences Finies, Eléments Finis ou Volumes Finis) demandent l'inversion d'une matrice tridiagonale (en ID, par discrétisation à 3 points de cette dérivée -4- seconde), ou tridiagonale par blocs tridiagonaux (en 2D/3D), pour obtenir une solution de basse précision du problème différentiel posé. A titre d'exemple, la simulation numérique des écoulements (de chaleur, et/ou de masse, et/ou de fluide) se fait très couramment aujourd'hui, sur ordinateur, par résolution numérique d'équations non linéaires, aux dérivées partielles en espace-temps. Le coeur des méthodes de résolution numériques se situe dans la résolution du terme de diffusion moléculaire, un terme de dérivée seconde, qui peut être mono-dimensionnel, bi-dimensionnel, ou tri-dimensionnel d'espace, selon la situation considérée (c 'est à dire un opérateur � Laplacien � pour employer un terme mathématique). Les méthodes industrielles pratiquées à ce jour sont presque toutes basées sur l'une des deux méthodes approximatives de discrétisation des dérivées, soit les Eléments Finis (EF), soit les Volumes Finis (VF). Les méthodes industrielles utilisent en général l'approximation la moins onéreuse (coût calcul et espace mémoire), définie par une discrétisation à 3 points de la dérivée seconde dans chaque direction de l'espace. On connaît aussi des méthodes dites � d'approximation globale � ou � spectrales �, qui se basent sur tous les points du maillage pour réaliser une approximation de la solution ou des solutions. Elles sont plus précises mais nécessitent aussi des calculs beaucoup plus lourds, ce qui fait qu'elles ne sont quasiment pas utilisées dans un cadre industriel. Un but de l'invention est de fournir un procédé et un système de résolution numérique par ordinateur permettant une meilleure précision allant possiblement jusqu'à la précision spectrale, tout en limitant le volume de calcul nécessaire, en particulier pour des traitements industriels basés sur des logiciels qui n'utilisent pas de supercalculateurs de très haute puissance. Un autre objectif est de permettre l'implémentation d'un tel procédé de façon efficace, simple rapide, souple et robuste, et possiblement en réutilisant tout ou partie des applications ou systèmes existants. -5 Exposé de l'invention L'invention propose un nouveau procédé de calcul électronique de simulation ou prédiction, utilisant une nouvelle méthode pour la résolution numérique précise (spectrale) du � Laplacien � une, deux dimensions et/ou trois dimensions (2D/3D). Cette méthode fournit un schéma de résolution numérique de type spectral, donc très précise, de la dérivée seconde, et, par extension, du Laplacien 2D/3D. Tous les schémas par approximation locale connus à ce jour (Différences Finies, Eléments Finis ou Volumes Finis) demandent l'inversion d'une matrice tridiagonale (en ID, par discrétisation à 3 points de cette dérivée seconde), ou tridiagonale par blocs tridiagonaux (en 2D/3D), pour obtenir une solution de basse précision du problème différentiel posé. Dans le cas monodimensionnel, l'invention permet ainsi de réaliser ce calcul de simulation en résolvant un opérateur � d2/dx2 � (un � Laplacien � monodimensionnel) par un calcul différent du second membre du système numérique, et en inversant uniquement une matrice simple � Ll �, qui peut être la même que celle utilisée dans les méthodes d'Eléments Finis (EF) ou de Volumes Finis (VF). Cette matrice � Ll �, dont l'inversion doit être calculée à chaque nouveau calcul, est une matrice relativement légère qui est elle aussi basée sur une discrétisation à 3 points de la dérivée seconde. Cette inversion se fait avec un faible surcoût, et une occupation mémoire à peine plus grande. Mais le résultat est extrêmement précis, l'erreur pouvant être le zéro machine (c'est à dire en dessous de la précision de l'ordinateur) pour les solutions de nature polynomiale. Ainsi, selon l'invention, il est proposé un procédé pour simuler, par calcul électronique réalisant une résolution numérique, les valeurs � U � d'une grandeur physique � u � ou chimique en un ou plusieurs points d'un maillage dans un espace � x �. Cette résolution est faite d'une part à partir d'une modélisation de ladite grandeur � u � par une équation différentielle du second ordre � f(x)=d2u/dx2 � construisant une première matrice � F �, et d'autre part à partir d'une modélisation � So � des conditions aux limites. -6- Ce procédé est du type comprenant un calcul d'inversion d'au moins une deuxième matrice � L �, dépendante du maillage, déterminée à partir de ce maillage et de l'équation différentielle pour former une équation sous la forme : � L . U = F + So �. L'inversion de cette matrice permet alors d'obtenir des valeurs de résultat sous la forme : � U = L-1 . ( F + So) �. Un aspect de l'invention comprend un procédé pour simuler, par calcul électronique et/ou informatique réalisant une résolution numérique, les valeurs � U � d'une grandeur physique � u � ou chimique en un ou plusieurs points d'un maillage dans un espace � x �. Cette simulation ou prédiction se fait à partir d'une modélisation de ladite grandeur � u � par une équation différentielle du second ordre � f(x)=d2u/dx2 � construisant une matrice de modélisation � F �, et à partir d'une matrice de frontières � So � modélisant les conditions aux limites. Ce procédé est du type comprenant un calcul d'inversion d'au moins une première matrice � L �, dépendante du maillage, pour obtenir des valeurs de résultat sous la forme : � U = L-1 . ( F + So) �. Selon l'invention, ce procédé utilise en outre une deuxième matrice � L1 � et une troisième matrice � A � qui sont issues d'une factorisation de cette première matrice � L � en un produit � Ao.L1 � d'une quatrième matrice � Ao � avec cette deuxième matrice � L1 �. Ces matrices sont agencées de façon telle que : - la deuxième matrice � L1 � est creuse et ne dépend que d'un nombre de points inférieur au nombre total de points du maillage, et la quatrième matrice � Ao � est inversible pour fournir la troisième matrice � A �, laquelle est définie analytiquement en fonction d'un nombre de points du maillage supérieur à celui de ladite deuxième matrice � L1 �. Une équation de modélisation sous la forme � L . U = F + So � équivaut alors à l'équation : � (Ao . L1) . U = F + So � soit : � L1 . U = (Ao)-1 . F + (Ao)-1 . So � et donc : � L1 . U = A . F + A . So � puisque : A = (Ao)-1 ou bien: �L1.U=A.F+S� Si l'on note :S=A.So -7- Dans le procédé selon l'invention, un calcul d'inversion de � L1 � fournit alors les valeurs de résultat sous la forme : �U=(L1)-1.(A.F+A.So)�. ou bien : � U = (Ll)-1 . (A . F + S) � Si l'on note : S = A . So Cette dernière équation fournit alors directement par un simple calcul analytique les valeurs numériques � u(xp) � de la grandeur à simuler � u �, pour chacun des points � xp � du maillage. L'invention permet ainsi d'obtenir des résultats d'une grande précision numérique, en particulier du fait que beaucoup de points du maillage sont utilisé dans la troisième matrice � A �. Pour cette raison, la précision obtenue est souvent de l'ordre du zéro machine des ordinateurs actuels, et est du même ordre que les méthodes spectrales connues dans l'état de la technique, et considérablement meilleure que les méthodes locales connues. En outre, l'invention permet aussi d'obtenir cette précision par des calculs très légers donc très rapides, en particulier du fait que seuls certains des points du maillage sont utilisés dans la deuxième matrice � L1 �. En effet, cela permet d'obtenir que cette deuxième matrice � L1 � comprenne un grand nombre de zéros, ce qui simplifie considérablement les calculs numériques que l'ordinateur devra réaliser pour en faire la transformation. En effet, ce sont ces calculs d'inversion pour transformer une matrice en sa matrice inverse qui représentent une charge de travail très importante dès que le nombre de points dans le maillage grandit. Or, la troisième matrice � A �, de son côté, est définie de façon analytique ce qui permet à l'ordinateur de la calculer très simplement, surtout en comparaison des calculs d'inversion de la deuxième matrice � L1 �. Ainsi, le grand nombre de points utilisé dans la troisième matrice � A � n'augmente donc pas ou peu la charge de calcul, même lorsque l'on utilise un très grand nombre de points dans le maillage. La charge de calcul, et donc les performances de vitesse du système réalisant le calcul, reste donc similaire ou peu augmentée par rapport aux méthodes locales connues dans l'état de la technique et en tout cas nettement moins que les méthodes spectrales connues. -8- Le compromis entre précision et performance est donc grandement amélioré. En particulier, nombre d'applications industrielles qui utilisent actuellement des méthodes locales à trois points pourraient gagner considérablement en précision tout en conservant les mêmes performances de vitesse à machine constante. En outre, ces mêmes applications devraient ainsi pouvoir mieux profiter des améliorations régulières de performances du matériel pour fournir des résultats plus rapides combinés avec de grandes améliorations de précision. En outre, il est à noter que la deuxième matrice � L1 � est la même que celle qui est déjà utilisée dans les méthodes locales connues. Cette identité, combinée à la simplicité des calculs complémentaires nécessités par l'invention, permet en outre une mise en oeuvre simplifiée au sein des logiciels ou des existants voire des systèmes existants, ce qui devrait en permettre une implémentation et une maintenance robuste et à moindre coût. Selon d'autres particularités de l'invention pouvant être combinées entre elles : - la deuxième matrice � L1 � est structurée en une matrice d'approximation locale portant, pour chaque point � xp � du maillage, sur trois points du maillage, incluant typiquement ce même point de calcul � xp � et deux points � xp_1 � et � xp+1 � situés de chaque côté de ce point de calcul. - la troisième matrice � A � est structurée en une matrice d'approximation spectrale portant sur tous les points du maillage. - le procédé comprend le remplacement de la modélisation analytique de la fonction � f(x) � représentant les variations de la grandeur � u � à calculer, dans la matrice de modélisation � F �, par une approximation polynomiale prenant les mêmes valeurs que cette fonction de modélisation aux différents points du maillage. - le procédé est appliqué au calcul électronique ou informatique d'une simulation par résolution numérique au sein d'un espace à deux ou trois dimensions. -9- le procédé est appliqué à simulation par résolution numérique au sein d'un espace présentant sur au moins un côté une limite où la valeur de la grandeur à calculer � u � est indépendante de ses évolutions au sein dudit espace, configuration dite de Dirichlet. - le procédé est appliqué à simulation par résolution numérique au sein d'un espace présentant sur au moins un côté une limite où la valeur � u(x_1) � et � u(x+l) � de la grandeur à calculer � u � dépend des évolutions de cette grandeur � u � au sein dudit espace, configuration dite de Robin et possiblement de Neumann. Selon un autre aspect de l'invention, il est proposé un procédé de conception d'un système électronique ou informatique de calcul mettant en oeuvre le procédé de simulation par calcul électronique ici décrit. Ce procédé de conception comprend une conception ou une adaptation de la configuration matérielle d'au moins un composant électronique de ce système pour exécuter de façon matérielle (par opposition à logicielle) au moins une partie du calcul, par exemple par des circuits dits de � logique câblée �. Une telle architecture matérielle de calcul est particulièrement bien adaptée au calcul de construction de l'inverse de la deuxième matrice � L1 �. Une telle architecture est aussi très avantageuse pour le calcul de construction de la troisième matrice � A � (aussi notée � L3sP �), qui est presque pleine et peut être d'une taille très importante et donc consommatrice de puissance de calcul, même pour un calcul de type analytique. Selon encore un autre aspect de l'invention, il est proposé un système électronique ou informatique mettant en oeuvre le procédé exposé ici. Selon une particularité, un tel système comprend une architecture de calcul en parallèle agencée pour exécuter en parallèle, en tout ou partie, le calcul d'une pluralité de différentes valeurs utilisées pour le calcul d'inversion de la deuxième matrice � L1 � et/ou pour le calcul de construction de la troisième matrice � A � (aussi notée � L3sP �). Une telle architecture parallèle est en effet particulièrement adaptée au types - 10 - d'opérations réalisées pour ce calcul, que ce soit pour l'inversion de la deuxième matrice � Ll � mais aussi pour la construction de la troisième matrice � A � ou � L3SP � en raison de sa taille et du fait qu'elle est presque pleine. D'autres particularités et avantages de l'invention ressortiront de la description détaillée d'un mode de mise en oeuvre nullement limitatif, et des dessins annexés sur lesquels : la FIGURE 1 illustre un exemple de simulation d'un espace à simuler numériquement, dans un mode de réalisation de l'invention en une dimension pour simuler ou prédire la température en différents points à l'intérieur d'un mur entre deux milieux de températures différentes ; la FIGURE 2 illustre schématiquement un exemple d'organigramme de fonctionnement d'un système informatique ou électronique de calcul mettant en oeuvre le procédé de simulation numérique selon l'invention ; la FIGURE 3 illustre schématiquement la structure et les interactions de données mises en oeuvre dans le mode de réalisation de l'invention de la FIGURE 1. Dans le domaine de la physique au sens large, on essaie autant que possible d'identifier des lois mathématiques formant un � modèle mathématique � qui représente le fonctionnement, l'état et les variations d'un système donné. Il peut s'agir par exemple du mouvement d'un corps pour prédire ou simuler son comportement, ou de l'évolution de la température pour déterminer les flux de chaleur dans un milieu tel que par exemple un mur d'une maison qui est soumis à une température extérieure différente de la température intérieure. Ces lois mathématiques se présentent sous la forme d'équations 30 mathématiques de différents types, exprimées entre les différentes grandeurs physiques concernées. 2947071 - 11 - Notion d'Equations différentielles Dans un certain nombre de cas, il s'agit d'équations exprimant une relation entre une grandeur physique et une ou plusieurs dérivées de cette grandeur par rapport à une autre grandeur. 5 Dans le cas d'un corps en déplacement vertical, la dérivée de la hauteur � z � par rapport au temps � t � est la vitesse verticale � vZ �, c'est à dire la façon dont � z � varie au fur et à mesure du temps. Cette dérivée est notée � z' � ou dz et vaudra ici : v = dz dt z dt De la même façon, la dérivée de cette dérivée sera appelée � dérivée 10 seconde � et sera ici la variation de la vitesse verticale � vZ �, c'est à dire l'accélération � aZ �. Cette dérivée seconde sera notée � z" � ou � z(2) � ou d2z et vaudra ici : a = d2z dt2 Z dt2 Selon les lois de l'inertie, l'accélération d'un objet est déterminée par sa masse � m � ainsi que par la somme � f � des forces qui lui sont 15 appliquées, soit : Pour une chute libre dans le vide, la seule force appliquée est la pesanteur qui présente une valeur : fp =m.g, où � g � est constante dans un endroit donné, soit environ 9,81 m/s2 à la surface de la Terre. L'équation du mouvement est alors une équation différentielle en dérivée seconde par m.a = f 20 rapport au temps d'un type très simple : d 2 z m. = m.g dt2 d Pour une chute avec frottement fluide, une autre forme encore simple d'équation différentielle du second ordre apparaitra si l'on adopte l'hypothèse que les frottements créent une force � ff � proportionnelle à la vitesse � vZ � par un coefficient � b �, soit selon la formule : ff=b.vZ 25 La somme vectorielle des forces appliquées à l'objet est donc : � fp + ff � Ces forces s'appliquant en sens opposés, le mouvement est donc régi par l'équation : m.az = fp - ff soit : m.az + fz = fp A partir des formulations ci-dessus pour � fp � et � ff �, l'équation devient donc : m.dt z + b. dt = m.g - 12 - Cette équation est aussi une équation différentielle en dérivée seconde, dite aussi du second ordre. Méthodes de Résolution analytique Dans le cas de la chute d'un corps soumis à la pesanteur mais dans le vide, l'équation ci-dessus représentant la variation de l'altitude de ce corps est du type : d 2 z m. = m.g dt2 La science mathématique nous permet de résoudre cette équation différentielle seconde de façon � symbolique � ou � analytique �, en modifiant la forme de l'équation pour obtenir une nouvelle forme nous permettant de calculer précisément et facilement la valeur numérique de � z �, soit l'équation donnant analytiquement le résultat : z = 2-.g.t2 A partir de cette équation � résolue �, lorsque l'on connaît l'environnement (ici la valeur de � g �), on peut obtenir directement la valeur de � z � à n'importe quel moment c'est à dire pour n'importe quelle valeur de � t �, par un calcul numérique simple et direct. Méthodes de résolution numérique. Des équations différentielles en dérivée seconde sont fréquentes par exemple dans les systèmes comportant des flux, tels que des flux de chaleur, des flux magnétiques ou des flux de fluides liquides ou gazeux, causés par des variations d'environnement physique. Il s'agit par exemple des équations de Navier-Stokes, ou des équations de Poisson. Or, dans de nombreux cas et contrairement aux exemples cités ci-dessus, la science mathématique concernée ne nous permet pas de résoudre complètement les équations mises en jeu, ou seulement sous une forme nécessitant des calculs très lourds même pour des moyens informatiques. Dans de tels cas, on utilise alors des méthodes de résolution dites numériques pour obtenir une valeur numérique approchée, mais dont la 1 - 13 - précision peut parfois être aussi bonne que celle du calcul numérique direct issu d'une résolution analytique, en particulier lorsque l'erreur théorique est trop faible pour être visible par rapport à la précision de l'ordinateur utilisé. Les méthodes numériques peuvent s'appliquer à des systèmes comprenant une ou plusieurs équations, et portant sur une ou plusieurs variables. Ainsi qu'illustré en FIGURE 1, l'exemple d'une grandeur � u � correspondant à la température d'un mur (1) entre deux milieux à des températures � u_ � et � u+ � différentes pourra être résolu par une méthode numérique dite � 1D � c'est à dire � une dimension �, portant sur la température à différentes profondeurs dans l'épaisseur du mur. Cette dimension unique � x � sera représentée sur un axe (10) perpendiculaire à la surface du mur (1), et variant par exemple entre une valeur � -1 � sur une face du mur, � 0 � au milieu de l'épaisseur, et � +1 � sur l'autre face du mur. Des données représentant les différents paramètres de ce systèmes sont entrées dans un ordinateur (19) agencé pour exécuter une simulation numérique, qui fournit en sortie des données (199) représentant les valeurs calculées pour la grandeur � u � aux différents points (100 à 105) de l'épaisseur � x � du mur (1). D'autres systèmes à simuler, plus complexes, peuvent comporter plusieurs dimensions, et seront résolus par des méthodes de résolution numériques de même nature, mais appliquées à plusieurs � dimensions �, par exemple en 2D ou 3D. Ainsi que précisé plus loin, l'invention peut aussi s'appliquer à la résolution numérique de systèmes à plusieurs dimensions, pour en augmenter la précision ou en alléger la charge de calcul. Ces méthodes de résolutions numériques utilisent une � discrétisation � du système modélisé, c'est à dire que le système est représenté par un certain nombre de points distincts formant un maillage. Ainsi, l'épaisseur du mur (1) de l'exemple est un espace unidimensionnel continu, qui sera découpé en un certain nombre de segments formant les différents points (100 à 105) d'un maillage, plus ou moins fin en fonction de la précision recherchée et des capacités de calcul informatique disponibles. - 14 - Dans le cas d'un espace à une dimension notée � x �, une équation différentielle du second ordre à résoudre pourra porter sur une grandeur � u �, par exemple la température, en fonction de la valeur de � x �, par exemple la position dans l'épaisseur du mur. Une telle équation pourra s'écrire sous la forme d'un opérateur différentiel dxu égal à une fonction f(x), c'est à dire : dxu = f(x). Cette équation devra en outre être complétée par des conditions aux frontières, par exemple la température qui est considérée comme constante dans cet exemple (constituant une configuration de Dirichlet), et connue de chaque côté du mur, qui sont notées : u(x=�1) ou � u+ � Si la valeur de � x � vaut -1 et +1 aux deux extrémités de l'espace étudié, par exemple respectivement sur chacune des deux faces du mur. Les mêmes calculs et méthodes de l'art antérieur ou de l'invention sont extrapolables à d'autres types de conditions de frontières. Ils sont aussi extrapolables à un espace à � d � dimensions, dans lequel l'opérateur différentiel considéré sera l'opérateur appelé � Laplacien � du type : d a2 Lz-1ax, avec � d � le nombre de dimensions et � x; � la coordonnée dans Z la direction � i �. Système à résoudre Le système d'équations à résoudre dépend de la discrétisation, c'est à dire du nombre de points dans le maillage choisi et de leur position. Un exemple de maillage peut être défini sous la forme d'un nombre (N+1) de points de coordonnées � xp �, avec l'indice � p � variant de 0 à N. Les positions de ces points peuvent définis par exemple en choisissant une l o i de répartition du type : xp = ûcos(p. N) , avec xo=-1 et xN=+l. La résolution cherche à trouver la valeur de la grandeur � u � pour chacun de ces N+1 points. Dans l'équation initiale d2u ù= f(x) chacun des dx2 deux membres de l'équation est discrétisé sous la forme d'un vecteur colonne, chacun composé d'une suite de valeurs correspondant aux 2947071 - 15 - différents points du maillage, ce qui donne une équation reliant � U � et � F � entre eux, si l'on note : • vecteur � U � pour l'inconnue dans le membre de gauche de l'équation (le � premier � membre) : U = (uN (xp), p = O..N) 5 • vecteur � F � pour les données du problème (la fonction connue et les conditions aux frontières) dans le membre de droite de l'équation (le � deuxième � membre) : F(� = (u_, f(xp), p =1..Nù1,u+) Par exemple pour une valeur de N=5 illustré en FIGURE 1, l'inconnue U � s'écrit sous la forme du vecteur colonne : u5(xo) u5(xl) 10 U= u5(x2) u5(x3) u5(x4) u5(x5) Pour le même maillage, les données du problème s'écrivent dans le second membre sous la forme des vecteurs colonnes : et F = , ou bien réunis en : Far = o .f (xi) f (x2 ) f (x3 ) f(x4) o u .f (xi) .f (x2 ) f(X3) f(X4) u+ 15 Pour un tel système à résoudre, les méthodes actuellement connues peuvent être décrites comme suit. Méthodes connues 20 Les méthodes numériques de résolution d'un opérateur différentiel se basent le pus souvent sur une discrétisation de cet opérateur, associée à un maillage du domaine continu sur lequel il agit. Les méthodes � par approximation locale � sont celles utilisées le plus souvent dans les codes industriels (c'est à dire systèmes ou logiciels de 2947071 - 16 - simulation utilisés dans l'industrie). Elles font appel à une discrétisation de bas niveau de précision. Elles sont de type � Elément Fini � (EF) ou de type � Volume Fini � (VF), ou parfois de type � Différence Finie � (DF). Elles s'appuient alors sur trois points voisins pour effectuer l'approximation de 5 l'opérateur dx2 . Cette approche s'est imposée pour minimiser les coûts calcul, ainsi que l'espace mémoire requis pour résoudre ce problème différentiel, qu'il faut parfois étendre à deux ou trois dimensions d'espace. Des méthodes de type � spectral � existent aussi, qui mettent en d oeuvre une approche de discrétisation de l'opérateur 2 qui fasse dx 2 10 intervenir tous les points du maillage. La solution numérique ainsi obtenue est de grande précision, la meilleure que l'on puisse obtenir avec le maillage choisi. Cette option présente l'inconvénient d'engendrer un coût de calcul beaucoup plus important, et n'est actuellement pas vraiment utilisée dans les codes industriels. 15 De telles méthodes connues sont exposées ici dans le cas où la valeur de � u � reste constante aux frontières (conditions de Dirichlet). Elles peuvent cependant aussi être utilisées avec les mêmes inconvénients lorsque la situation correspond à un cas plus général (conditions de Robin), par exemple lorsque les conditions aux limites dépendent des variations de 20 la grandeur à l'intérieur de l'espace. Dans le cas de la FIGURE 1, cela serait le cas si l'on considère que la pièce est suffisamment petite pour que la température � u � de son air (soit � u(xo) �) dépende de celle à l'intérieur du mur 1 qui chauffe ou refroidit la pièce. 2947071 -17- Méthodes par approximation locale - DF, EF ou VF (Différence Finie, Elément fini ou Volume Fini) Les méthodes d'approximation locale ou ;4. quantité est évaluée en chaque point du intérieur ait domaine on .r.ï, pour p = 1. - , 1., par une combinais son linémea. re des valeurs de . échantillonné sur un nombre réduit de: points On parlera alors d'tin schéma 3--point; s si lapproximation de ity~s'écri sous la t�-,>tin('. .t, itN cp (.rp+f d2itN.- (5'les ap, b.p et cp étant des ixieffieients constants établis par un chemin mathématique d.e2 prédable. La sohftidn numérique s'obtient en résolvant le système tit(éirire représenté par les équations (.6) et (). .11 y a (N 1) équat.k)n pc.tn autant d'inconnus, avec c:PaliKed lus deux conditions de frontière. ä etx (.. -- a+. (h) puis l'équation différuutielle discétsék.. n, or\ �J' p ( Ï) - r• d'f a A- x [) so u s tif péi 1\ s'c . ire dé deux. mile Fe, d il fén f (a) Pour la inéthodé dite Différence Finie (DE). Sp .... f (:e: ) p = 1 - 1, Le systdmc .matricief (1))-(7) à. is--i,.soudre sécrit donc , 'p = , .N , le vecteuf�'.ol)une ayant été défini pal: i3). En .Annexi 1 se trotivfi'.une iflus- trat:011 IMilléidqUe de ce système linaire, (b) Pour les méthodes dites dlément Fini (EF) et de Volume Fini (\TF . les sources 'Obtiennent par un traitement préaabk,. du problè.me de départ. Elles s'écrivent ta, + f i:2l,i ^t il- f (xp+1 ) p 1. 1, (101 les coefficients o(pä lp et ,yi., derxindant du choix KU ou \,.T. Le ?l tCi'iile' matriciel a, résoudre s'écrit alors (11) (LOù1 5 les vecteur (iolimiries tl+) Ct, F ayant été déit-na par �in note .eplé lé niembré (U f cie ce système est le 11v4ne pour EF et On voit donc que le calcul des valeurs de l'inconnue U, c'est à dire de la matrice des valeurs de la grandeur � u � aux différents points � xp � de - 18 - l'espace � x �, nécessite un calcul de transformation, par inversion, de la matrice située dans le premier membre de l'équation finale ci-dessus, soit l'équation (11) pour les méthodes locales connues. A chaque calcul d'une nouvelle configuration, par exemple en changeant le maillage (211) ou la modélisation (212) ou les conditions aux limites (213), un système de calcul doit donc réaliser le calcul d'inversion de cette matrice L1. Pour l'exemple numérique ici décrit, et dans le cas où les conditions aux limites répondent aux conditions dites � de Dirichlet � c'est à dire présentent une valeur constante, ces matrices ont la structure et les valeurs suivantes : [J)]. EF et. \'F 11() av()11S S13t(t >tii e)11(')It 5.2-1 -7.21 2 2 -3.62 L1 2 1.)l2 -3.62 2 9 -7.21 5.21 l 1- i) >.1,` 23. i2 lO ) 5.>3 37.3 f . : r 1É}.3 37.3 8.33 8,33 23. 3 o et ou )) 5,.2S 26, 2 4. '. 7,7 43.3 t.73 7.73 43.3 .7 4.89 26.2 5.2,S On voit que la matrice de gauche (L(°F) ou L1), qui est inversée pour chaque nouveau calcul, est assez � creuse �. C'est à dire qu'elle présente systématiquement de nombreux emplacements vides (valant zéro) qui allègent de façon importante le calcul d'inversion. La matrice de droite (L3(EF) ou L3("F)) présente elle aussi une telle structure creuse. Ces structures creuses sont dues au fait que l'approximation est locale, et ne concerne que quelques points (ici trois points). Le calcul est relativement léger grâce à cet aspect local, mais la précision en est limitée d'autant. 1. loo L>~. - 19 - La quantité Si. eainsi év;nuée "Jar (5 avec uiae .certaine precision, Itre .xp d( l'ordre de max, (.r.p 21_1 i et la solution numérique, l'-.�)urnie par (9) ou (Il), sera une approximation échntillom.A), 4.lé 41ÿ.i`!. Pouf améliorer la précisieai du résultat. 11 faut soit rapprocher les poh-its voisins, donc augmenter IV et réduire ainsi la quart- tite :Lnax.p ù soit passer à. un schéma 3-ponts, en aoutant. les contributions 21 (a,ip+2) dans (5). Dans .i-.ie dernier cas erreur de-vie .ut de Fordreie rnay.p (x en ceffitrepartie, on augm ente le nomIn'e d'ééments di' la matrice a. hrvfg-ser.ä donc le coii.t de It.i. réOUtIO11 dit système, comme le voit ci-après. Méthodes par approximation globale (spectrales) Les méthodes d'approximation globale,. ou "méthodes spectrales". Ces méthode :sont très précises, car .elle sont construif?f: lit tes pour qiu toutes las iiaicurs le�:eppnaxienation drl.-Ilée Sceur1f:1e: en toiat point. conduit à lexpression survante de l'approximation dé 4.yittc dérivée seconde, 2) a Lec,) p ==== ä , i,q 02� Les (N +1)2 quantités (E12,; sont les 4l'41rnimt mat ce D2, la matrice D représentarà hl matlIce dn déiv( prerrilén. On eaf(natre que la udatnni ,fem.r'iW la dari:fiée su:onde de tout peduiiOnté u(,r� dc degré nuélanfum N x, d2 �D2 ) ) zzzii :, da-q=p 0, N u(x.) (13) 'ai .-rsénéiral, dans toute situation réeUe (industrielle) la solution n'est pas 13olynomiale. nlai.8 tÈ !tait_ réglilarité pOltr (-11tlii la pi'enSi.Cin dit re..'tallat S'aalelOre exponentiellement 4-..eff.-.c 1M-vir Iable 2�. Le ..stèille linéaire ;4éerit L1,8PI .... ti où 'a matràx! COnentIt à partir de D2 (Tieffliaile Pour l'exemple ici décrit, et dans le cas où les conditions aux limites répondent aux conditions dites � de Dirichlet � c'est à dire présentent une valeur constante, les matrices concernées ont la structure et les valeurs suivantes : Pour ,approximation spectrale, uoi:ui avons tas'a.t~aai'd .a nlatrice de data. 3ô1 ~;_ilI SeVel:)lide )Job' 11.6 -68,1 -23,6 17,6 -8 2L3 -3:1.:? 12.7 -3.69 2.21 ù(i.95:1 -1.85 7.32 -10.1 3.79 -1.91 0.714 -1.91 5.79 -10.1 7.32 -1.85 6,95 2.21. .169 12.7 -31.5 21.3 8 17,f 23 6 1� ).S 6.1 II 6 2947071 -20- L 'à. çzr Le calcul de l'inversion de cette matrice � LSP � est donc considérablement plus conséquent que pour les méthodes locales. C'est cette charge de calcul qui représente une durée bien plus importante, et/ou nécessite des moyens de calcul nettement plus performants. Méthode selon l'invention - Exemple d'un mode de réalisation La méthode selon l'invention sera ici principalement décrite dans l'hypothèse d'un espace à une seule dimension, ainsi qu'illustré en FIGURE 1, mais peut-être appliquée aussi à plusieurs dimensions comme expliqué plus bas. La FIGURE 2 illustre un exemple d'organigramme de fonctionnement d'un système informatique ou électronique de calcul mettant en oeuvre le procédé de simulation numérique selon l'invention, avec une phase de test selon que les conditions aux limites sont ou non du type Dirichlet. - 21 - La FIGURE 3 présente un schéma représentant la structure et les interactions de données mises en oeuvre dans l'exemple de mode de réalisation de l'invention ici décrit. Application Mono-dimensionnelle Nous considérons ici 1331 espace à. Iule seule direction, étant la coordonriiise représentant 1.a position Lf diiiiei1s1+:i11irelle) dans cette direction, et un probléme différentiel à. résoudre, d.}2 ?'�,r) ] ---1. I [, la solution cherchée devant -vérifier deux conditions de frontic'l'e, luit. 12 cnaque extrémité = I r•t 1, dont la fi:wnnitition la plis générale est celle dite de type Robin, Elle s'écrit = ~.. t~= f E d+ oi3 1t s et ^r'� sont (les ('r_IIistt I�tes. Les o.+ et ,"i�� sont sllpposs ~ "liier le c itère dit d`ellipticitc', iii savoir A Chaque situation physique correspondra 1311 jeu. de - a?. eurs des a+, + et 'f, On peut cependant distinguer deux configurations i:> .. tex soit de pur flux idite de type NeLlInanli'.� Si 1eS rit ÿ'annidlelit, Selon la situation physique à traiter, on peut aussi avoir une iondition de Dirichlet sur l'une des extrémités �l'un des ?f est nul). et mixte, ou encore de Neumann, sur l'autre d'autre des f est non. nul)ä Notre présentation se fera en considérant que 1' in a. soit de conditions de .Dirichlet, ! 3 f Soit des conditions mixtes, Imposées aux deux extrémités.. Elle permettra de voh�.:on-mient traiter lès autre cas. Application Multi-dimensionnelle Tout ce qui suit est extrapolable à un espace bi- ou tri-dimensionnel dans d lequel l'opérateur considéré sera le Laplacien d2 , avec d=2 à 3 et -1 dx~ 2 rata � x; � est la coordonnée dans la direction � i �. -22-Présentation mathématique Soit A le nombre d'intervalles que nous niions prendre pour ptirtitionnor le domaine [---1, 1], et . = --cos: (pli p = les coordonnees des points corieituant le .niaillage adopté. Voici ce que sont ces positions pour .N = 5, = -1, ;ri = 0809017, := ---0.309017, v.;+ = 0.309017, = 0,,109017, ;r; = 1.. On note donc ;-eo = --el et xjv = 1 les positions extrêmes. Désignons par aN (.;O le polynôme de degré ni.kyliintun N (t]1 :représente la s(;i.ution. nninér-iqu 11-ip1�KUe permettant (Dt,ceeder à rr(y). Ce polyitôtite:'•'erd déterminé. sin ce nuli;luge de (N 1) points, ( 4- eIN:Xp), polir p • • • , A , 1.m uod uisc>ns maintenant N t 1) Fonet;ous dessa t,(r), j 0, , fnu ( ti^,)ns linéaires pal mou eaux, defmies sur le maillage de la façon ù,.u:ant,c, D'aboi:d (N. ù t) fonctions sont construites sur (.r le-s noeuds int,erné,,,. X [Xi , X:i+i 'éj(:i)) +++- 0. Cenc qui. ctuiceruc le noeud en .-eij est 1�:outrée en Fg. ces ]:onctiuLs sont de I : \.pe Elles Figure 1: i+a. fon(lTft, n v.:e(j.) t-te (5à,, où rapproxftnation FF. aussi, • • N �1) : Puis Fonctctti,s sont construites sur les noeuds en frontière, f ;fr. , XII O, to,,.)) et ;cj'v , x 1 :7,0, vN(x) := 1. 7), Ces deux fore hotrs sont. ;es "detnt-dntpaux'' représentés en. Fig. 2 Elles \--d:ri:lent Elles sont donc nulles sur tous _es points iot,e leurs di]. domaine. Eigure 2; Les _fonctions d'ess :. a: +'{11;1' 1 et vN(:r) dc l'approximation -23- Invention - Exemple en conditions de Dirichlet AVeC COndäitiOnS Cri f'reilkiere 1-10.915 aVO.RS û 1) inconnues f déternriner, à savoir eust'F1i?>a' des Tl Y.,) Pour N. 1, Elles sont déliinies corn:me s>lution du Mt( le Traitement .clu terme O "'.r (à gauche dans la parenthèse de l'équation (8)) Le terme f ib. e transforme par 11te.gra,to,,(terme de gauche dans la parenthèse) p .N du. up(:r) ci"; .)) .ep(x .)) -' (lx e=l. f (12"p ftri.N x= -1 iLr. ^'.I:r ' La quantité prise dans laeeoiade est nulle, rat' toutes les fonctions vp(x) sont nulles en et. il reste pour CL terme (hi- d-uN Et r.- (ix (irace à la deiimtion n? ionetions cpi, eakule di(ectument. il devient. ie Ulerllbre deirdite de Cette eqUatiOII tie )-r air,,y dl, 'xP+ 'du A.- ä rp d :t' ce qui donne quelle que soit lapproximation �.iN(.r) choisie pour d(.r',i. ? l:1 a� vN ) 4- hi, ïiN (rp) 4- , d t' (x?,+1 , qi o û (n. *-e d ) = 1, - û 1 (IO) Eau, :rr, Traitement du terme vve4f (;e) de (terme de droite dans la parenthèse de l'équation (8)) La fonction i t est ra.rement connue anayti(1uement, et, de toutes fitçons eifectiier une int:égrab.E-. est touours (oui:eux. De plus, cette fonction est: I)rokibiemeiu. obtenue elle .nienee par voie numérique. Elle est donc ()Il tout état de cerise connue par ses (N -r 1) .. lettes nodales, à saxon_;' les .,fï, �::ep) , p= 0, ù N, Un peut alors introduire l'approximation polynorniale3 f:y(:r) de j'x), li existe en effet un polynôme riniqné, de degré N, etni prend les -valeurs fp aux points de coordonnées r:p. Il s'agit ici d'un changement intéressant par rapport à la méthode N connue des Eléments Finis, dans laquelle on écrivait : fN(x)=Efjvj(x) , c'est à dire une approximation linéaire par morceaux de la fonction � f(x) �. - 24 - Dans la méthode selon l'invention, cette approximation polynomiale s'écrit comme une somme pondérée de Polynômes de Lagrange � /p(x) � : fx' 'ri p , . propriété fp( aisé Lient (123, On (on tatf. qi i a oïl u de x) in e f(i) alt point- xp, in:, avoc 1 icrue j (�)u'.r devient fi cr,r)f(x)dx ,rr tT.iir) 1i(iv) p __.1 l' ù:t. ï..y t "ùo Système discret résultant de (8) Regroupant (9) E4 (1,1) nous. pouvons écrire léquation (S) sous la formé finale suivante, itN ) (rp -t irp+ , ri.,(.� d. Y=0 ics up, .n E; par tO),5 - 25 - Ecriture matricielle Pour écrire la version matricielle de l'équation (15), on introduit le vecteur colonne des conditions aux limites : .- . (?0 iD>rl l~ C.> 'wu lui: ecr rf! 15) sf3u, ,.6 iwt133E'i117J.itiol_E tous)ie, ,st ( J U [ t t s] }la=. 3:t:,•s mat): mat): s' Lirj'' ' t 1011 [ gel )1 1N;fY [t`cJflfOIO', if-y.1, édL)ti: (tL)i' 3 VE#ii'U , IFOd ,, t~B.S1[;i .tY vt:e1)4•)U,, C3;EUR';S wltiva Ut:)), Iyc)ts . i CE ve°tac ~t 'c -t'IaC d(C la natal i011 ,fly :)IltlOdtlii pal' E. Cette équation (21) correspond ainsi à l'équation générale s'écrivant : L1. U = S+A. F (22) et comprenant : - une matrice (ou vecteur colonne) de modélisation � F � représentant la fonction de modélisation � f(x) � combinée aux différents points xp du maillage ; une deuxième matrice � L1 � définie sous une forme � L1Dir � ne dépendant que de trois points du maillage ; et - une troisième matrice � A � définie sous une forme � L3sP,Dir � dépendant de tous les points du maillage La matrice � S � de l'équation (22) correspond au produit � A . So � énoncé dans l'exposé général présenté précédemment.20 - 26 - Deuxième matrice � Ll � La matrice � L1Dir � (c'est à dire la deuxième matrice � Ll � du cas général) située dans le premier membre (à gauche) de cette équation (21), est celle qui est à inverser pour le calcul de chaque nouvelle configuration. Les valeurs de cette matrice Ll°ir sont définies analytiquement comme suit : 1 dv dv. (L1')1 = f dx dx dx pour p=1 à N-1 et j=0 à N (23) Dir (L)OJ = 0 pour j=1 à N (L Dir),j =o pour j=0 à N-1 Dir (L1 )11 ù et (L1Dir )NN ù Avec les positions dans la matrice notées comme suit : col.0 col.1 col.2 ... col. p-1 col.p 0,0 0,1 1,0 1,2 2,1 2,2 col.p+1 col.N N,N li.p-1 p-1,p-1 P-11P li.p P,P-1 p, p p,p+1 Ii.p+1 p+ 1, p p+ 1, p+ 1 Ii.N Le calcul analytique de cette formulation (23) donne, pour chaque valeur � p � des indices de lignes autres que 0 et N : Dir _ 1 Dir = -1 Dir _ 1 (L1 )P,P-1 et (L1 )P,P et (L1 )P,P+1 xp ùx 1 xp+1 ùx 1 xp+1 ùx P-1 P-1 P - 27 - Dans le cas de l'exemple de la FIGURE 1, pour N=5, la matrice L1Dir présente les valeurs suivantes : %3. 2i -7,24 2 2 -:$,{32 1 .1>2 :L(52 -3.G2 2 2 724 5.24 Les � 1 � mis dans ses première et dernière lignes traduisent les deux conditions de Dirichlet à gauche et à droite de l'espace � x �. Cette matrice � L1(Dir) � est exactement celle (� L1EF,' �) que l'on obtenait avec les schémas locaux connus, par Elément Fini et Volume Fini, dans l'approximation linéaire par morceau de u(x), qui donnait un plus bas niveau de précision. Pour chaque valeur de l'indice de ligne � p �, on voit que les seules valeurs non nulles ne portent que sur les indices � p-1 � et � p � et � p+1 �. Pour chaque point � xp �, Les valeurs non nulles de cette matrice dépendent donc uniquement des points � xp_1 � et � xp � et � xp+1 � situés immédiatement de chaque côté. Pour chacun des points du maillage, la matrice � L � (ou � L1 �, ici notée � L1Dir �) ne dépend donc que de trois des points du maillage. Du fait de sa structure creuse, la méthode selon l'invention nécessite une faible charge de calcul, du même ordre de grandeur que celle des méthodes locales actuellement utilisées dans beaucoup de logiciels industriels. - 28 - Troisième matrice � A � La matrice L3(sP,Dir) située dans le deuxième membre de cette équation (21) correspond à la troisième matrice � A � de l'exposé général présenté précédemment. Cette matrice n'est à calculer qu'une seule fois, au moins pour un maillage déterminé. Pour un même maillage, la matrice � L3(sP,Dir) � du membre de droite n'a pas besoin d'être recalculée et est obtenue analytiquement comme suit, pour chacune des positions (L3sP,Dir)p1 correspondant aux indices respectivement � p � et � j � de lignes et de colonnes : L3 cLr p \' û 1. , NT.t On complète avec les deux lignes nulles (la première et la dernière) des conditions de frontière. Ceci donne, pour N=5: (1 O o O O O ùL2:1 27.8 9.6 2.:1 1.67 . -0.708 - 0.654 5.01 4 37 0.53 -2 . 66 1.03 1.03 -2,66 9.53 137 5t_ 1 -0.65 t - 00708, 1..61. -2.54 0.0 27.8 ùL2-1 f) t_t ii 0o o Avantages Avec les schémas connus par approximation locale, par exemple Elément Fini et Volume Fini, dans l'approximation linéaire par morceau de la fonction de modélisation � f(x) �, nous avions des matrices L3(EF,Dir) et L3(vF,Dir) tridiagonales sur leurs lignes internes, c'est à dire ne dépendant que de trois points. Avec l'invention, on voit que cette matrice L3(sP,Dir) est presque pleine, et porte donc sur la totalité des points du maillage, ce qui concoure à la précision spectrale du schéma obtenu. Cette précision est obtenue en particulier grâce au choix de l'approximation polynomiale fait pour l'équation (11), ainsi que grâce au choix et à la construction de la factorisation permettant d'arriver à ces deux matrices : _ � L1Dir � qui est creuse, et _ � L3sP,Dir � qui est presque pleine tout en étant connue analytiquement, et est en outre constante pour un même maillage. tiF'.L?sr 1 10 i 2947071 - 29 - Pour le calcul de chaque nouvelle configuration, le calcul des valeurs de l'inconnue � U �, c'est à dire de la matrice des valeurs de la grandeur � u � aux différents points � xp � de l'espace, ne nécessite qu'un calcul de 5 transformation, par inversion, de la matrice � Ll°'r � (la deuxième matrice � Ll �) située dans le premier membre de l'équation ci-dessus (21). Cette inversion conduit alors à l'équation suivante, donnant directement les valeurs de la grandeur � u � à calculer : 10 La matrice � L1Dir � , qui est inversée pour chaque nouveau calcul, est de structure � creuse �, de la même façon que pour les méthodes d'approximation locale de l'état de la technique. C'est à dire qu'elle présente systématiquement de nombreux emplacements vides (valant zéro) qui allègent de façon importante le calcul d'inversion. 15 Plus précisément, cette matrice est tridiagonale, c'est à dire que seules les valeurs situées sur la diagonale et sur une position de part et d'autre de cette diagonale présentent des valeurs non nulles (encadrées en traits gras et respectivement en traits fins dans le tableau ci-dessus). Cette structure creuse de la matrice à inverser � Lloir � concoure à la 20 rapidité de calcul de l'invention, puisque les nombreuses valeurs zéro réduisent grandement la charge de calcul à réaliser, à l'instar des méthodes locales connues, alors que la précision obtenue est de niveau spectral. Les performances de vitesse obtenues par l'invention sont similaires ou légèrement augmentées par rapport à ces méthodes locales, et reste du 25 même ordre de grandeur. Par rapport aux méthodes spectrales, qui sont d'une précision comparable à celle de l'invention, les performances de vitesse de l'invention sont donc considérablement améliorées. En outre, il est à noter que la définition analytique de cette deuxième matrice � Ll � selon l'invention (ici � L1Dir �) est la même que pour les méthodes d'éléments finis de l'état de la technique. L'implémentation de son calcul, ainsi que son inversion, peuvent donc être réalisés de la même - 30 - façon que pour les applications existantes, voire utiliser tout ou partie des applications existantes. Précisions obtenues La méthode selon l'invention permet ainsi de mettre en oeuvre une approximation de type spectral, et donc d'obtenir une bien meilleure précision que les méthodes locales, tout en n'effectuant la majorité du calcul que sur trois points, ce qui permet de garder une charge de calcul du même ordre que ces méthodes locales. Pour ces raisons, on pourrait qualifier cette méthode de � schéma spectral 3-points �. ;`iislS t`cl1El[}a:'tf131S #'.S numériques t1 N,de : ~~,,1 telles i:{ll ,#'h sont dormes par es. schémas s d'a.ppr,C;11 atim21 locale. i4 par le sellerait 3-ptaiu t~s d l'invention, en calculant l'erreur maximale commise, _L" :: Ii1Fta tl. i c i t! P i . )eu types de solutions (.) i\ .sa -,4 . D'abord u.uei solution est G11� ,.°:e_: I obino ri 3 e. de ` degré 2l.{ a' ,i ï) pour laquelle on s'attend, par a re at.~c}la (1:3), Méthode selon Invention - extension aux conditions de Robin Cas des cieux conditions de Robin (4) Id. les \'aieurs nodales YIN (.ro) et 'uA (xN ) deviennent des inconnues supplémentaires 11 .'.ii .va. falloir :11iterininer en ajoutant d(,D.N: equations (11_13 tiennent compte des den.x. conditions de Robin (4).. Partons de la relation ()'7i) que nous prc)longeons usqu'en p ,,,: 0 et p :::: iloi ù 0 , p ù 0. - , N. 2uN ctr'2 Id. interviennent les deux t'Onctions de borci, v)(r; et v N(x definies par (7), et t.)résent ées eu Fig. 2. Ou. traite ces tenues éotunie décrit en Section 3.1. Ceei conduit I� un ensemble de relations discrètes. On. retrou:ve d'ab(ed les relations (ionnees 1Ja 115) ä à savoir c .u:y (.'u-].! bv .u,..y t,, r té,j\ ( x1)+1 ) ===, 1, , N - 1. (17') .._Eh 0111pléli'111 (IS�-'( deltX I�d'atiolls de bOl (1111 arri Vol jt'S �Ondli Ions; (-].i, / vof:ii) f(x) , . .. ft, 1 (.rN- )N CN 7+f .r)f th:. 19'1 , .. -1 Les. coefficients de (17) sont déjà cOl.11111S, titsltilt'ti p;tA' (10). Pour les atitl'eS CA effiCi.elt.S, Co:::: = = 3? 1 - 3:'O et UN, 3333 333, 333: I. ù xxù.i Les intégr.tes 'des seconds membres de (18) eli (19 se cacnient comme c(3iies :de (17), en y mettant la déicoituposition poynorniaide 111. (le f (Iev3ennent 21 (':Co + il.,.N. Ti TN�..ci A,. - f, ù .L do: , (,.iv- (À.Lf.Nt) J. = J=0 -32-Ecriture matricielle Rangeons les inconnue_ tENt_2 pi. a.)E1si que les. 1 0 0 tl f+ Le tiy ti`.11l(? (17)-(19b et SOI13tioii 1It3Y11(?11(131(' ti i i'IIS`t 1It .`:do h C le la façon tiàiiv�3Yttt', ',R rats r iR.r:t,~3 + (",1'6 '1 OÙ k '4 lIl La partie intérieure de cette matrice � L1R°b �, c'est à dire sauf la première et la dernière ligne, est la même que celle de la matrice � L1Dir � correspondant aux conditions de Dirichlet. Seules les deux positions extrêmes de la diagonale sont différentes. Le calcul analytique des valeurs de cette matrice � L1R°b � donne alors les formules suivantes : Pour chaque valeur � p � des indices de lignes autres que 0 et N : 1 Rob _ -1 Rob _ 1 x _ x et C~1 ~P~P X p+1 x et (L1 ~P~P+1 X p+1 X P P-1 P+1 P-1 P+1 p L1Rob C ) p,p _ -1 -33-et pour les valeurs 0 et N : Rob -1 Rob = -1 (Li )00 _ xl + 1 et (L1 )NN 1 _ xN 1 NOUS 0.VOD,S, pour = 5, 1_ lWCI.lHnt deh Membres (1( gauche 5."24+ 5.24 5.24 7.24 2 2 --3.62 1.112 L62 --3.62 2 2 -- 7.24 5.24 5,24 524+ >ID Cli,`.int c1E 111[_'E1i~:~ï{ i: s_ droite c c' (17)-(.1.9.y, 5.' 6 4.62 -0,764 1.382 -0,272 0,122 1.24 2 r .8 9,67 2.54 1.61 11.08 1 .h> 1 -0.654 5.111 4,37 E1.53 -2.66 1.03 28 1 1.(}3 ---'266 9.5: 4.37 5.01 ---11.054 -0.70S 1.61 • 2.54 11.61 27.8 ---1.24 0 1 22! -0."2 f 2 11.382 -0.764 4.62 5.46 '1l?a.is !E:, E?lén1 ilt? L'invention peut donc s'appliquer pour programmer ou concevoir des systèmes de calcul présentant une meilleure précision que les systèmes actuels fonctionnant par approximation locale, ou une meilleure rapidité que les systèmes actuels fonctionnant par approximation globale, voire cumuler ces deux avantages. Ces avantages peuvent être obtenus pour des solutions entièrement logicielles fonctionnant sur des ordinateurs de type polyvalent monoprocesseur, mais aussi apporter des avantages pour une mise en - 35 - oeuvre sur des systèmes fonctionnant en tout ou partie en mode de calcul parallèle. L'invention peut aussi être intéressante à implémenter dans des systèmes comportant des composants électroniques de calcul programmés sous forme matérielle, pour réaliser tout ou partie des calculs de la méthode selon l'invention. Bien sûr, l'invention n'est pas limitée aux exemples qui viennent d'être décrits et de nombreux aménagements peuvent être apportés à ces exemples sans sortir du cadre de l'invention. REVENDICATIONS1. Procédé pour simuler, par calcul électronique réalisant une résolution numérique, les valeurs � U � d'une grandeur physique � u � ou chimique en un ou plusieurs points d'un maillage dans un espace � x �, à partir d'une modélisation de ladite grandeur � u � par une équation différentielle du second ordre � f(x)=d2u/dx2 � construisant une matrice de modélisation � F �, et à partir d'une matrice de frontières � So � modélisant les conditions aux limites, ledit procédé comprenant un traitement de données réalisant un calcul d'inversion d'au moins une première matrice � L �, dépendante du maillage, pour obtenir des valeurs de résultat sous la forme : � U = L-1 . ( F + So) �, caractérisé en ce qu'il utilise une deuxième matrice � L1 � et une troisième matrice � A � qui sont issues d'une factorisation de ladite première matrice � L � en un produit �Ao.L1� d'une quatrième matrice � Ao � avec ladite deuxième matrice � L1 �, lesquelles sont agencées de façon telle que : - la deuxième matrice � L1 � est creuse et ne dépend que d'un nombre de points inférieur au nombre total de points du maillage, et _ la quatrième matrice � Ao � est inversible pour fournir ladite troisième matrice � A �, laquelle est définie analytiquement en fonction d'un nombre de points du maillage supérieur à celui de ladite deuxième matrice � L1 � le calcul d'inversion de � L1 � fournissant alors les valeurs de résultat sous la forme : � U = (L1)-1 . (A . F + A . So) �. 2. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que : - la deuxième matrice (L1) est structurée en une matrice d'approximation locale portant, pour chaque point � xp � du maillage, sur trois points du maillage incluant ledit point de calcul et deux points � xp_1 � et � xp+1 � situés de chaque côté dudit point de calcul; et - la troisième matrice (A, L3SP) est structurée en une matrice d'approximation spectrale portant sur tous les points du maillage. 2947071 - 37 - 3. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce qu'il comprend le remplacement de la modélisation analytique de la fonction (� f(x) �) représentant les variations de la 5 grandeur (u) à calculer, dans la matrice de modélisation (F), par une approximation polynomiale prenant les mêmes valeurs que ladite fonction de modélisation aux différents points du maillage. 4. Procédé selon l'une quelconque des revendications précédentes, 10 caractérisé en ce qu'il est appliqué au calcul électronique ou informatique d'une simulation par résolution numérique au sein d'un espace à deux ou trois dimensions. 5. Procédé selon l'une quelconque des revendications précédentes, 15 caractérisé en ce qu'il est appliqué à la simulation par résolution numérique au sein d'un espace présentant sur au moins un côté une limite où la valeur de la grandeur à calculer est indépendante de ses évolutions au sein dudit espace, configuration dite de Dirichlet. 20 6. Procédé selon l'une quelconque des revendications 1 à 5, caractérisé en ce qu'il est appliqué à la simulation par résolution numérique au sein d'un espace présentant sur au moins un côté une limite où la valeur (u(x=-1) ou u(x=1)) de la grandeur à calculer (u) dépend des évolutions de ladite grandeur (u) au sein dudit espace, configuration dite de Robin et 25 possiblement de Neumann. 7. Système électronique ou informatique mettant en oeuvre le procédé selon l'une quelconque des revendications 1 à 6. 30 8. Système selon la revendication précédente, comprenant une architecture de calcul en parallèle agencée pour exécuter en parallèle, en tout ou partie, le calcul d'une pluralité de différentes valeurs utilisées notamment pour le calcul de construction de la troisième matrice (A, L3sP)- 38 - 9. Procédé de conception d'un système électronique de calcul mettant en oeuvre le procédé selon l'une quelconque des revendications 1 à 8, caractérisé en ce qu'il comprend une conception ou une adaptation de la configuration matérielle d'au moins un composant électronique dudit système pour exécuter de façon matérielle au moins une partie du calcul, notamment de construction de la troisième matrice (A, L3sP)