Méthodes spectrales pour la simulation PDE : théorie, mise en œuvre et quand les utiliser

Les méthodes spectrales sont des techniques numériques d’ordre élevé pour résoudre les équations aux dérivées partielles. Contrairement aux méthodes de différences finies, qui se rapprochent des dérivés par le biais de pochoirs locaux, les méthodes spectrales représentent la solution avec des fonctions de base qui s’étendent sur l’ensemble du domaine de calcul.

Les choix courants incluent les modes de Fourier pour les problèmes périodiques et les polynômes de Chebyshev ou Legendre pour les domaines bornés et non périodiques. Lorsque la solution est suffisamment lisse, ces approximations globales peuvent atteindre une précision très élevée avec relativement peu de degrés de liberté.

Cet avantage s’accompagne d’importantes limitations. Les méthodes spectrales fonctionnent mieux sur des domaines simples avec des solutions fluides et des conditions aux limites bien définies. Les discontinuités provoquent des oscillations, une géométrie complexe affaiblit la commodité des bases globales et l’intégration explicite dans le temps peut devenir sévèrement restreinte à mesure que la résolution augmente.

Quelles sont les méthodes spectrales ?

Une méthode spectrale se rapproche d’une fonction comme une somme pondérée de fonctions de base lisses :

u(x) ≈ Σ cₙ φₙ(x)

Les fonctions φₙ forment la base spectrale, tandis que les coefficients cₙ décrivent la contribution de chaque mode. Au lieu de stocker uniquement des informations locales, chaque coefficient peut affecter l’approximation sur le domaine complet.

La base doit refléter la structure du problème. Les fonctions de Fourier sont naturelles pour les domaines périodiques car elles satisfont déjà à la périodicité. Les polynômes de Chebyshev et Legendre sont couramment utilisés sur des intervalles finis tels que [-1, 1].

Les méthodes spectrales ne sont pas simplement des schémas de différences finies très élevées. Ils suivent une stratégie d’approximation différente. Les méthodes locales construisent la solution à partir d’informations dans des cellules ou des éléments à proximité. Les méthodes spectrales utilisent des modes globaux qui peuvent décrire une fonction fluide avec une efficacité remarquable.

Trois principales formulations spectrales

Les méthodes spectrales peuvent appliquer la PDE gouvernante de plusieurs manières.

Méthodes de collocation

Une méthode de collocation nécessite que l’équation différentielle soit maintenue aux points de grille sélectionnés. Ces points sont souvent connectés à la règle de base et de quadrature, tels que les points de grille de Fourier ou les nœuds Chebyshev-Gauss-Lobatto.

La collocation est populaire car elle convertit les dérivés en opérations matricielles ou en transformations spectrales. Il est également relativement facile à combiner avec des termes non linéaires évalués dans l’espace physique.

Méthodes de Galerkin

Une méthode de Galerkin projette le résidu de la PDE sur les fonctions de base sélectionnées. Le résidu n’a pas besoin de disparaître à chaque point, mais il doit être orthogonal à l’espace d’approximation.

Cette formulation fournit une forme faible naturelle et peut préserver d’importantes propriétés d’énergie ou de symétrie lorsqu’elle est conçue avec soin.

Méthodes TAU

Une méthode TAU modifie les équations sélectionnées associées aux modes d’ordre le plus élevé afin que les conditions aux limites puissent être imposées. Il est étroitement lié à l’approche de Galerkin mais gère les contraintes différemment.

Les logiciels spectraux modernes peuvent masquer une grande partie de ces détails d’implémentation. Cependant, les utilisateurs doivent encore comprendre comment la formulation choisie représente les conditions aux limites et les opérateurs différentiels.

Méthodes spectrales de Fourier

Les méthodes de Fourier représentent une solution périodique comme une somme des modes trigonométriques :

u(x) ≈ Σ ûₖ exp(ikx)

La différenciation devient particulièrement simple dans l’espace spectral :

dûₖ/dx = ik ûₖ

Un solveur numérique peut donc transformer la solution en coefficients de Fourier, multiplier chaque coefficient par le nombre d’onde approprié et transformer le résultat en espace physique.

La transformée de Fourier rapide réduit le coût de ces opérations à environ O(N log N). Cela rend les méthodes de Fourier très efficaces pour les simulations périodiques impliquant des ondes, la turbulence, la dynamique des fluides et la formation de motifs.

Un dérivé de Fourier en Python

La fonction suivante calcule la première dérivée d’une fonction périodique échantillonnée sur une grille uniforme :

import numpy as np
from scipy.fft import fft, ifft, fftfreq

def fourier_derivative(values, domain_length):
    """Return the first derivative of periodic grid data."""

    number_of_points = len(values)
    spacing = domain_length / number_of_points

    wave_numbers = (
        2.0
        * np.pi
        * fftfreq(number_of_points, d=spacing)
    )

    spectral_values = fft(values)
    spectral_derivative = (
        1j
        * wave_numbers
        * spectral_values
    )

    return ifft(spectral_derivative).real


# Example
length = 2.0 * np.pi
points = 128

x = np.linspace(
    0.0,
    length,
    points,
    endpoint=False
)

values = np.sin(3.0 * x)
derivative = fourier_derivative(values, length)

exact_derivative = 3.0 * np.cos(3.0 * x)
error = np.max(np.abs(derivative - exact_derivative))

print("Maximum error:", error)

Pour une fonction périodique lisse qui est bien résolue par la grille, la dérivée peut être extrêmement précise. La méthode évite également l’erreur de troncature associée à un pochoir à différences finies courtes.

Méthodes spectrales de Chebyshev

Les modes de Fourier ne conviennent pas lorsque la solution n’est pas périodique. Les méthodes de Chebyshev offrent une alternative courante sur un intervalle fini.

Les nœuds Chebyshev-Gauss-Lobatto sont définis par :

xⱼ = cos(πj / N),  j = 0, ..., N

Ces points se regroupent à proximité des points de terminaison. Le regroupement améliore l’interpolation polynomiale et aide à contrôler les grandes oscillations qui peuvent se produire avec une interpolation d’ordre élevé équidistant.

Le même clustering crée également un défi de pas de temps. Le plus petit espacement près des limites devient beaucoup plus petit que l’espacement moyen de la grille, ce qui peut imposer des limites de stabilité restrictives aux méthodes explicites.

Construire une matrice de différenciation de Chebyshev

L’implémentation suivante crée la matrice de différenciation de Chebyshev standard du premier ordre :

import numpy as np

def chebyshev_differentiation_matrix(order):
    """Return Chebyshev nodes and first derivative matrix."""

    if order == 0:
        return (
            np.array([1.0]),
            np.array([[0.0]])
        )

    indices = np.arange(order + 1)
    nodes = np.cos(np.pi * indices / order)

    coefficients = np.ones(order + 1)
    coefficients[0] = 2.0
    coefficients[-1] = 2.0

    coefficients *= (-1.0) ** indices

    node_matrix = np.tile(
        nodes,
        (order + 1, 1)
    )

    differences = (
        node_matrix.T
        - node_matrix
    )

    ratio_matrix = np.outer(
        coefficients,
        1.0 / coefficients
    )

    derivative_matrix = (
        ratio_matrix
        / (
            differences
            + np.eye(order + 1)
        )
    )

    derivative_matrix -= np.diag(
        np.sum(
            derivative_matrix,
            axis=1
        )
    )

    return nodes, derivative_matrix


# Example
order = 32
x, derivative_matrix = (
    chebyshev_differentiation_matrix(order)
)

values = np.exp(x)
numerical_derivative = derivative_matrix @ values
exact_derivative = np.exp(x)

error = np.max(
    np.abs(
        numerical_derivative
        - exact_derivative
    )
)

print("Maximum error:", error)

La matrice de différenciation est dense car chaque fonction de base globale influence l’intervalle entier. Une multiplication directe matrice-vecteur a un coût d’environ O(N²).

Pour des problèmes unidimensionnels modérés, cela peut toujours être pratique. Les simulations plus grandes peuvent utiliser des méthodes basées sur des transformations, des reformulations clairsemées, une décomposition de domaine ou des bibliothèques spécialisées.

Pourquoi la convergence spectrale peut être si rapide

Le principal avantage des méthodes spectrales est leur taux de convergence pour des solutions fluides. Une différence finie d’ordre faible ou une méthode d’éléments finis converge normalement de manière algébrique :

Error ≈ C N⁻ᵖ

La valeur de p dépend de l’ordre de la méthode. Par exemple, doubler le nombre de points dans une méthode de second ordre peut réduire l’erreur d’environ un facteur de quatre lorsque la solution est dans la plage de convergence asymptotique.

Pour une solution analytique, une approximation spectrale peut converger géométriquement ou exponentiellement :

Error ≈ C exp(-αN)

Cela signifie que l’augmentation du nombre de modes peut réduire l’erreur beaucoup plus rapidement que l’augmentation de la résolution d’une méthode locale d’ordre inférieur.

Une introduction pratique avec des exemples numériques est disponible dans Méthodes spectrales dans Matlab.

L’exigence de douceur

La convergence exponentielle ne se produit pas pour toutes les fonctions. Cela dépend de la régularité de la solution exacte.

• Une solution analytique peut produire une convergence géométrique ou exponentielle.
• Une solution infiniment différentiable mais non analytique peut produire une convergence plus rapide qu’algébrique sans un taux exponentiel fixe.
• Une solution avec seulement un nombre fini de dérivés produit normalement une convergence algébrique.
• Une solution discontinue crée des oscillations de Gibbs et supprime le principal avantage d’une base lisse globale.

Les coefficients lisses ne garantissent pas une solution lisse. Les coins, les conditions aux limites incompatibles, le forçage discontinu, les interfaces de matériaux et les données initiales singulières peuvent réduire la régularité.

Avant de sélectionner une méthode spectrale, les chercheurs doivent examiner la douceur attendue de la solution plutôt que uniquement l’apparence de l’équation gouvernante.

Le phénomène Gibbs

Une expansion globale de Fourier ou polynomiale ne peut représenter un saut sans osciller à proximité. Ce comportement est connu comme le phénomène Gibbs.

Au fur et à mesure que le nombre de modes augmente, la région oscillatoire devient plus étroite, mais le dépassement maximum près de la discontinuité ne disparaît pas de la même manière que l’erreur de région lisse ordinaire.

Ces oscillations peuvent créer des concentrations négatives, des valeurs de pression non physiques ou des calculs non linéaires instables. Le filtrage peut les réduire, mais le filtrage supprime également les informations à haute fréquence et introduit la dissipation.

Les méthodes spectrales globales pures sont donc rarement le premier choix pour les lois de conservation dominées par les chocs.

Méthodes spectrales multidomaines

Une façon de préserver la précision spectrale est de diviser le domaine en sous-domaines. Chaque sous-domaine reçoit sa propre expansion spectrale en douceur.

Si une interface de discontinuité ou de matériau se situe exactement sur une limite de sous-domaine, l’approximation à l’intérieur de chaque région peut rester lisse. Conditions d’interface puis connectez les solutions de sous-domaine.

Cette approche est courante en astrophysique et en relativité numérique. Un examen détaillé est disponible dans le Méthodes spectrales pour la relativité numérique.

Les formulations multidomaines créent également un pont entre les méthodes spectrales globales et les techniques d’éléments spectraux.

Approches de capture des chocs

Plusieurs techniques ont été développées pour stabiliser les approximations spectrales à proximité des chocs. Ils comprennent le filtrage spectral, la viscosité de la disparition spectrale, la relaxation et la suppression périodique des modes à haute fréquence non résolus.

Des travaux récents sur la relaxation spectrale et la purge spectrale examinent comment des noyaux soigneusement conçus peuvent contrôler les oscillations tout en conservant des informations utiles à petite échelle. Un exemple est l’étude Méthodes spectrales nouvelles pour la capture des chocs et la suppression de Tygers dans Dynamique des fluides informatique.

Ces techniques peuvent améliorer une simulation spectrale, mais elles ne font pas un problème discontinu équivalent à un problème lisse. La méthode, la force du filtre, la résolution et les propriétés de conservation nécessitent toujours une validation minutieuse.

Restrictions de pas de temps

Une précision spatiale élevée ne supprime pas les limites de stabilité temporelle. En fait, les discrétisations spectrales peuvent produire de grandes valeurs propres qui rendent l’intégration explicite du temps restrictive.

Pour la discrétisation de Fourier de l’advection du premier ordre, le plus grand nombre d’ondes augmente proportionnellement à N. Une limite de stabilité explicite est donc souvent approximative :

Δt ∝ N⁻¹

Pour la discrétisation de Fourier de la diffusion, les valeurs propres augmentent comme le carré du nombre d’ondes :

Δt ∝ N⁻²

Le clustering de points Chebyshev rend les limites explicites plus restrictives. Pour les problèmes de première dérivée, la limite pratique peut évoluer à peu près comme N⁻². Pour les opérateurs de diffusion de la deuxième dérivée, cela peut devenir encore plus grave.

La condition exacte dépend de la PDE, du traitement des limites, de la formulation et de l’intégrateur de temps. Il ne doit pas être réduit à un seul exposant universel.

Intégration implicite et IMEX

Les méthodes implicites peuvent éviter les restrictions de stabilité les plus fortes associées à la diffusion linéaire ou à d’autres termes rigides. Les formules de Crank-Nicolson et de différenciation vers l’arrière sont des choix courants.

Une méthode IMEX traite implicitement les termes linéaires rigides et évalue explicitement les termes non linéaires ou moins restrictifs :

∂u/∂t = L(u) + N(u)

L’opérateur linéaire L peut représenter la diffusion, tandis que N contient une advection ou une réaction non linéaire. Cette structure est largement utilisée dans les logiciels Spectral PDE.

Une comparaison plus large est disponible dans le guide Méthodes d’intégration temporelle pour les solveurs PDE : schémas explicites ou implicites.

Termes non linéaires et alias

Les produits non linéaires créent des modes avec des fréquences supérieures à la résolution d’origine. Lorsque ces modes sont échantillonnés sur la grille existante, ils peuvent apparaître de manière incorrecte en tant que composants à basse fréquence. C’est ce qu’on appelle le crénelage.

Les solveurs pseudospectraux calculent généralement des dérivés dans l’espace spectral et des produits non linéaires dans l’espace physique. Avant de transformer le produit, ils peuvent appliquer la distribution.

La règle commune des deux tiers supprime les modes de Fourier les plus élevés après une multiplication non linéaire. Une autre approche renforce la représentation spectrale vers une grille plus grande, effectue la multiplication et tronque le résultat.

Sans négocier, une simulation peut devenir inexacte ou instable même lorsque la grille spatiale semble suffisamment fine.

Utiliser Dedalus

La construction d’un solveur spectral multidimensionnel complet nécessite une gestion de base, des transformations, des équations de limites, une distribution parallèle et une intégration temporelle. Le DeDalus Documentation décrit un framework Python conçu spécifiquement pour les simulations PDE spectrales.

Dedalus prend en charge les bases de Fourier et polynomiales, les problèmes de valeurs initiales, les problèmes de valeurs limites, les problèmes de valeurs propres et l’exécution parallèle. Il fournit également des outils basés sur Tau pour imposer des contraintes dans des domaines non périodiques.

Les utilisateurs doivent suivre la syntaxe de la version de Dedalus installée, car son API a changé entre les principales versions. Conceptuellement, le flux de travail reste cohérent :

1. Sélectionnez les coordonnées et les bases spectrales.
2. Créez des champs pour les variables dépendantes.
3. Définissez les équations et les contraintes de limites.
4. Sélectionnez un intégrateur de temps ou un solveur linéaire.
5. Définissez les tâches de résolution, de négociation et de sortie.
6. Exécutez les contrôles de convergence et de stabilité.

Méthodes spectrales et conditions aux limites

Les conditions aux limites périodiques sont construites naturellement dans une base de Fourier. Les conditions non périodiques nécessitent plus de travail.

Les conditions de Dirichlet ou de Neumann peuvent être imposées en remplaçant les équations de collocation, en construisant des fonctions de base qui satisfont déjà aux conditions ou en ajoutant des variables et des contraintes TAU.

L’approche sélectionnée affecte le conditionnement de la matrice et la structure du système final. Les conditions aux limites doivent donc être prises en compte lors de la sélection de la base, et non seulement une fois la discrétisation spatiale terminée.

La recherche sur les bases orthogonales des EDP dépendant du temps fournit des moyens supplémentaires de classer les systèmes de base et le comportement aux limites. Une discussion mathématique récente est disponible dans Foundations mathématiques des méthodes spectrales pour les PDE dépendantes du temps .

Méthodes spectrales vs Galerkin discontinu

Les méthodes de Galerkin spectrale et discontinue utilisent toutes deux une approximation polynomiale, mais elles répartissent différemment la base.

Une méthode spectrale traditionnelle utilise une base globale dans l’ensemble du domaine. Une méthode DG attribue une base polynomiale distincte à chaque élément et permet des sauts entre les éléments voisins.

Les méthodes d’élément spectral et de différence spectrale combinent la décomposition du domaine local avec une approximation d’ordre élevé à l’intérieur de chaque élément. Un exemple orienté python de la connexion entre une approximation locale d’ordre élevé et DG est décrit dans Quail : un code Galerkin discontinu léger open-source en python.

Quand les méthodes spectrales fonctionnent mieux

Une méthode spectrale globale est un choix fort lorsque :

• La solution attendue est lisse ou analytique.
• Le domaine est périodique, rectangulaire ou unidimensionnel.
• Les conditions aux limites correspondent à la base sélectionnée.
• Une précision spatiale élevée est plus importante que la flexibilité géométrique.
• Le problème peut utiliser des transformations basées sur FFT ou des matrices denses modérées.
• Les discontinuités et les interfaces matérielles nettes sont absentes.

Les applications typiques comprennent la propagation des ondes lisses, l’analyse de stabilité, les flux incompressibles sur des domaines simples, les modèles quantiques, la formation de motifs et les problèmes sélectionnés en géophysique et en astrophysique.

Quand une autre méthode est meilleure

Le volume fini ou les méthodes DG sont généralement plus naturels lorsque les chocs, les discontinuités de contact ou la conservation locale stricte dominent le problème.

Les méthodes des éléments finis et des éléments spectraux peuvent être plus appropriées pour la géométrie compliquée, le raffinement local et les limites irrégulières.

Les méthodes de différences finies d’ordre inférieur peuvent rester préférables lorsque la facilité de mise en œuvre, l’algèbre linéaire clairsemée et le comportement local prévisible sont plus importants qu’une précision extrême.

La décision doit être basée sur la régularité de la solution, la géométrie du domaine, les conditions aux limites, les restrictions d’échelle de temps et la sortie qui doit être prédite.

Un tableau de sélection pratique

erreurs de mise en œuvre courantes

Une erreur courante consiste à sélectionner une méthode spectrale uniquement parce qu’une grande précision est attendue. La solution doit d’abord être vérifiée pour la douceur.

Parmi les autres problèmes fréquents, citons :

• Utilisation des modes de Fourier pour des données non périodiques sans extension appropriée
• Ignorer l’alias dans les équations non linéaires
• Utilisation d’un pas de temps explicite qui viole la limite de stabilité spectrale
• Application des conditions aux limites de manière incohérente
• En supposant que toutes les fonctions fluides produisent le même taux exponentiel
• Utiliser trop de modes sans surveiller le conditionnement
• Interpréter les oscillations de Gibbs comme un comportement physique
• Sauter les comparaisons avec des méthodes d’ordre inférieur ou local

Comment valider un solveur spectral

Commencez par une fonction lisse dont la solution dérivée ou PDE est connue analytiquement. Augmentez le nombre de modes et mesurez l’erreur.

Pour un problème analytique, l’erreur devrait diminuer rapidement jusqu’à ce qu’elle atteigne les limites causées par la précision à virgule flottante, le conditionnement, l’erreur d’intégration temporelle ou une solution de référence insuffisamment précise.

Pour les problèmes non linéaires, répétez l’expérience avec et sans traitement. Vérifiez les quantités conservées, les résidus de limites et la décroissance des coefficients spectraux.

Une solution spectrale utile montre normalement des coefficients diminuant vers les modes résolus les plus élevés. Si les coefficients finaux restent élevés, la simulation peut être sous-résolue.

Guides connexes

• Qu’est-ce que la simulation scientifique et pourquoi cela est important ?
• Méthodes d’intégration de temps pour les solveurs PDE : schémas explicites ou implicites

Lectures complémentaires

• Trefethen, L. N. — Méthodes spectrales dans Matlab
• Matériau SIAM supplémentaire sur les approximations spectrales
• Iserles, A. — Fondements mathématiques des méthodes spectrales pour les PDE dépendant du temps
• Recherche sur la relaxation spectrale, la purge et la capture des chocs
• Grandclément et Novak — Méthodes spectrales pour la relativité numérique
• Documentation du projet Dedalus
• Caille : un code Galerkin discontinu léger open-source en Python

Conclusion

Les méthodes spectrales se rapprochent des solutions PDE avec des bases globales de Fourier ou polynomiales. Pour les solutions analytiques sur des domaines appropriés, elles peuvent atteindre une très grande précision avec beaucoup moins de degrés de liberté que les méthodes locales de faible ordre.

Leurs performances dépendent fortement de la douceur. Les discontinuités provoquent des oscillations de Gibbs, la géométrie irrégulière affaiblit la commodité des bases globales et les discrétisations de Chebyshev à haute résolution peuvent imposer de sévères restrictions de pas de temps explicites.

Les méthodes de Fourier sont particulièrement efficaces pour les problèmes périodiques, tandis que les techniques de Chebyshev et Legendre prennent en charge des domaines non périodiques délimités. L’intégration implicite ou IMEX du temps, la négociation et le traitement rigoureux des limites sont souvent nécessaires dans les simulations pratiques.

Lorsque le domaine est complexe ou que la solution contient des interfaces nettes, des éléments spectraux, des méthodes discontinues de Galerkin ou de volume fini peuvent fournir un meilleur équilibre. La méthode correcte est déterminée non seulement par la précision souhaitée, mais également par la régularité, la géométrie, les exigences de conservation et le coût de calcul.