TangenteX.com

Les solitons

Ondes et physique linéaire

Le concept d'onde est familier aux élèves du lycée et des classes préparatoires. Cependant, une onde, solution d'une équation de propagation comme l'équation de D'Alembert, présente quelques caractéristiques sur lesquelles on ne revient peut-être pas suffisamment. Une onde, telle qu'elle est modélisée mathématiquement, ne possède ni extension spatiale, ni extension temporelle : elle se propage indéfiniment, ce qui est physiquement assez ennuyeux. Pire, elle transporte une énergie infinie, ce qui est encore plus ennuyeux. Dans une page consacrée aux paquets d'ondes, j'ai abordé ces problèmes.

Les équations de propagation, D'Alembert ou Schrödinger, les équations de Maxwell ou celles de la physique quantique, sont des filles de la physique linéaire. Elles ont été établies en laissant de coté tous les aspects tenus pour négligeables, qui auraient pu introduire une non-linéarité, un sinus ou un terme quadratique par exemple. Sauf lorsque c'était impossible : les équations de Navier-Stokes sont non linéaires et en cela elles furent et sont toujours d'ailleurs, problématiques.

J'ai déjà abordé la non-linéarité dans plusieurs pages de TangenteX, par exemple dans l'étude des oscillateurs, et découvert une nouvelle branche de la physique née de cette non-linéarité : le chaos.

Il existe des objets physiques dont la physique est non linéaire, mais qui ne sont pas chaotiques. Ces objets sont observables dans la nature : ce sont les solitons. La première observation d'une drôle d'onde sur un canal, onde qui se propage à vitesse constante, sans se déformer et qui possède une extension spatiale finie, fut rapportée par John Scott Russel, un architecte naval, en 1834. Puis d'autres formes de ces ondes exotiques furent observées dans d'autres domaines de la physique, en particulier en astrophysique et en physique du solide. L'étude théorique de ces objets donna lieu à la construction de deux grandes familles de solitons : les solitons non-topologiques, ceux de l'hydrodynamique pour simplifier, et les solitons topologiques, ceux de la physique du solide, toujours pour simplifier. Les premiers sont décrits par l'équation de Korteweg-de Vries et les seconds par l'équation de sine-Gordon (SG) que j'ai déjà abordé sur TangenteX.

Dans cette page, j'aborderai l'équation de Korteweg - de Vries encore nommée KdV et les solitons non-topologiques ou hydrodynamiques, "non-topologiques" parce qu'ils modifient la topologie de leur milieu de propagation lors de cette propagation.

Tous les scripts Python cités dans cette page sont téléchargeables dans la bibliothèque de codes de TangenteX.

Qu'est-ce qu'un soliton

Un soliton est une onde particulière, d'amplitude constante et de faible extension spatiale, qui naît et se propage sur des distances qui peuvent être grandes dans un milieu dispersif.

Le terme "soliton", du à Zabusky et Kruskal en 1965, vient de la contraction des mots "onde solitaire", puisqu'il s'agit d'une onde qui se propage seule et sans déformation, aux effets dissipatifs près. Le suffixe "on" rappelle que l'on peut assimiler le soliton à une particule, comme un électron ou un proton : il possède une extension spatiale restreinte, il transporte une énergie finie et possède un analogue d'une masse.

Etablir l'équation de Korteweg-de Vries

D'où vient l'équation KdV

L'équation KdV est une équation de mécanique des fluides, plus précisément d'hydrodynamique, qui modélise la propagation d'ondes de grande longueur d'onde et de faible amplitude relative, à la surface d'un liquide.

Comme je l'ai écris plus haut, c'est un architecte naval écossais qui observa en 1834 un curieux phénomène sur un canal devenu célébre : une onde solitaire qui se propageait sur une grande distance (plusieurs kilomètres) sans que son amplitude et sa forme ne varient. Ce phénomène fut modélisé bien plus tard (en 1895) par deux physiciens hollandais Diederik Korteweg et Gustav de Vries.

Ils sont partis de l'équation fondamentale de l'hydrodynamique, l'équation d'Euler pour un fluide parfait (incompressible, non visqueux) :

\( \rho  \dfrac{\partial \overrightarrow{v}}{\partial t} + \rho (\overrightarrow{v}.\overrightarrow{\nabla}).\overrightarrow{v} = \rho \overrightarrow{g} -  \overrightarrow{\nabla} p \)

avec \( \overrightarrow{\nabla}.\overrightarrow{v} = 0 \)

Leur modèle s'articule sur ce système et sur certaines autres équations comme l'équation de conservation de la masse, la condition à la limite au fond et la condition à la limite cinématique. Ils posèrent également des hypothèses : négliger la tension superficielle à la surface du liquide,  supposer le liquide irrotationnel (il n'y a pas de tourbillon). Enfin, ils supposèrent que la solution de l'équation devait être non divergente puisque cette solution existait, ayant été observée expérimentalement.

Pour simplifier les choses, fort complexes comme d'habitude en mécanique des fluides, ils travaillaient dans le plan vertical, sur un modèle bidimensionnel.

Ils obtinrent une équation assez compliquée, qu'ils s'empressèrent d'adimensionner en choissant des variables caractéristiques : longueur, amplitude, célérité et temps. Ils ont également fait quelques hypothèses d'échelle pour pouvoir négliger certains paramètres : non-linéarité pas trop forte, propagation en eau peu profonde (la profondeur est faible par rapport à l'extension spatiale, ce qui est très relatif dans le cas d'un tsunami qui se propage sur des milliers de km dans un océan profond en moyenne de  4 km).

Bref, après toutes ces manipulations, ils produisirent l'équation KdV que nous allons examiner.

L'expression de l'équation KdV

Dans sa forme la plus épurée l'équation KdV s'écrit :

\(  \dfrac{\partial u}{\partial t} + 6u\dfrac{\partial u}{\partial x} + \dfrac{\partial^3 u}{\partial x^3} = 0 \)           (1)

où u est une fonction solution recherchée, x et t les coordonnées resp. spatiale et temporelle adimensionnées.

Vous la verrez le plus souvent écrite avec la notation simplifiée :

\( u_t + 6uu_x + u_{xxx} = 0 \)               (2)

Une brève analyse de l'équation KdV

Examinons la forme de cette équation. Nous observons trois termes :

  • Le premier terme \( \dfrac{\partial u}{\partial t} \) est le terme de propagation.
  • Le seconde terme  \( u\dfrac{\partial u}{\partial x} \) est le terme non-linéaire de l'équation.
  • Le troisième terme \( \dfrac{\partial^3 u}{\partial x^3} \) est son terme dispersif.

La forme d'un soliton non topologique résulte de l'équilibre entre le terme non-linéaire et le terme dispersif dans l'équation. Nous verrons tout à l'heure comment observer cet équilibre.

On démontre des caractéristiques intéressantes de cette équation. Elle conserve la masse. Elle conserve aussi l'énergie, propriété mise ne pratique en optique et laser, domaines dans lesquels les solitons ont une grande importance. Enfin, l'équation KdV n'est pas invariante sous la transformée de Lorentz...

Solution analytique de l'équation de KdV

L'équation (1) admet une infinité de solutions analytiques de la forme :

\( u = \dfrac{c}{2}sech^2 \left( \dfrac{\sqrt{c}}{2}(x - ct - x_0) \right) \)          (3)

où sech désigne la sécante hyperbolique \( sech(x) = \dfrac{1}{cosh(x)} \), c la célérité de l'onde et x0 la position initiale de son maximum.

Le script KdV_UnSolitonAnalytique.py permet de visualiser la propagation d'un soliton non topologique déterminée par l'équation KdV. Sur le graphe ci-dessous, il s'agit d'un soliton dont la célérité initiale c = 1.0 et la position initiale de son maximum est fixée à x0 = 8.

Soliton non topologique - Analytique c = 1.0

Un examen rapide de l'équation (3) nous permet d'affirmer que l'amplitude du soliton dépend de sa célérité. Nous pouvons d'ailleurs le confirmer en modifiant la valeur de c dans le script, en conservant tous les autres paramètres inchangés. Nous obtenons le graphe suivant :

Soliton non topologique - Analytique c = 0.5

L'amplitude du soliton a été divisée par deux. Sur l'animation vous constaterez également que la vitesse de propagation est inférieure. A propos de vitesse de propagation, vous pourrez aussi vérifier que la vitesse de propagation est supérieure à la vitesse c paramétrée. On dit que le soliton KdV est "supersonique".

Sur le graphe ci-dessus, vous pouvez observer un autre phénomène : la largeur du soliton de célérité c = 0.5 est plus grande que celle du soliton de célérité c = 1.0. La largeur d'un soliton croît inversement proportionnellement à sa célérité, selon la formule \( l = \sqrt{\dfrac{2}{c}}\).

Résoudre KdV par la méthode pseudo-spectrale

Quelques mots sur les méthodes pseudo-spectrales

Remarquons d'abord que dans une équation aux dérivées partielles (EDP), comme notre équation KdV, nous observons la présence d'opérateurs, en l'espèce les opérateurs de dérivation dans le temps, dans le terme \( u_t \) et dans l'espace, dans les termes \( 6uu_x \) et \( u_{xxx} \).

Fondamentalement, la méthode spectrale consiste à transformer la recherche des solutions d'une EDP en un calcul des valeurs propres et des vecteurs propres de l'opérateur spatial. Ce dernier calcul consiste à déterminer le "spectre" de l'opérateur, d'où le nom de ces méthodes. On part du principe que le calcul du spectre de l'opérateur sera plus précis que la résolution de l'équation par d'autres méthodes. Dans beaucoup de cas, ce principe est avéré, mais ce n'est pas toujours le cas, en particulier lorsque le calcul du spectre de l'opérateur est inabordable analytiquement.

On applique alors les méthodes pseudo-spectrales. Comme pour les méthode spectrales, il s'agit de trouver un ensemble de fonctions orthogonales qui formeront la base des solutions de notre équation et sur lesquelles notre opérateur spatial sera appliqué. Mais on procédera différemment en ne calculant pas directement le spectre de l'opérateur.

Le principe des méthodes pseudo-spectrales est basé sur la possibilité de séparer l'opérateur temporel et l'opérateur spatial dans l'équation KdV, en posant :

\( u_t  = -6uu_x - u_{xxx} \)

On peut transformer cette EDP en un système de deux EDO dans le domaine spatial en appliquant une transformée de Fourier. Je vous renvoie au livre de Claude Aslangul "Des mathématiques pour les sciences" §8.2.2 pour découvrir ou redécouvir la relation très utile entre dérivée et transformée de Fourier.

Sur le plan numérique, il y a grand intérêt à utiliser une transformée de Fourier, pour laquelle il existe des algorihtmes de calcul très efficaces (voir la FFT), pour calculer les dérivées, plutôt que les méthodes aux différences finies avec lesquelles on perd rapidement de la précision.

Il y a cependant une condition stricte à l'utilisation des méthodes pseudo-spectrales : le domaine de définition doit être périodique pour pouvoir utiliser les transformées de Fourier, et la fonction u(x,t) recherchée doit être de carré sommable, pour que les intégrales ne divergent pas. On sera donc amener à "périodiser" le domaine d'intégration.

A ceux que l'aspect mathématique de la méthode intéresse, je recommande la lecture du §8 de l'ouvrage "Méthodes mathématiques pour les sciences physiques" de Jean-Michel Bony.

Le script

Pour calculer les dérivées spatiales, je vais utiliser la fonction diff() du fftpack de scipy, à ne pas confondre avec la fonction diff() de scipy. Pour marquer cette différence, je la nommerai psdiff(), comme c'est l'habitude en calcul numérique avec scipy. Vous trouverez la documentation de cette fonction ici.

Le script KdV_UnSolitonPseudoSpectral.py se présente ainsi :

from scipy import cosh, sqrt, linspace
from scipy.integrate import odeint
from scipy.fftpack import diff as psdiff
from matplotlib.pyplot import figure,xlim,ylim,grid,xlabel,ylabel,tight_layout,show, title,plot
import matplotlib.animation as animation

# Fonction sécante hyperbolique
def sech(x):
    y = 1./cosh(x)
    return y

# Solution analytique de l'équation KdV
def KdVAnalytique(x,t,c,x0):
    u = (c/2.)*sech((sqrt(c)/2.)*(x - c*t - x0))**2
    return u
   
# Définition de l'EDP en employant la méthode pseudo-spectrale
def SystemeKdV(u,t,L):
    # calcul des dérivées en x d'ordre 1 et 3
    ux = psdiff(u, period=L)
    uxxx = psdiff(u, period=L,order=3)
    # calcul de la dérivée temporelle d'ordre 1 (définition de l'EDP)
    ut = -6*u*ux - uxxx
    return ut
   
# Définition du domaine spatial
L = 80.0                # période spatiale
Nx = 256                # nombre de pas de discrétisation du domaine spatial
x0 = 0
x = linspace(x0,L,Nx)

# parametres temporels
tmax = 200
t0 = 0
Nt = 225
t = linspace(t0,tmax,Nt)
  
# Paramètres du soliton initial
c1 = 0.5; shift1 = 0.15*L
u0 = KdVAnalytique(x,t0,c1,shift1)

# intégration numérique de l'équation KdV par la méthode pseudo-spectrale
y = odeint(SystemeKdV, u0, t, args=(L,), mxstep=5000)
 
# fonction de tracé de l'animation
def KdV(i,x,data,line):
    line.set_data(x,data[i,:])
    return line,

# tracé de l'évolution
fig1 = figure(figsize=(8,6))
title("Soliton de Korteweg de Vries - Méthode pseudo-spectrale")
xlim(0,L)
ylim(0,(c1/2.0) + 0.1)
grid(True)
xlabel("X")
ylabel("U(x,t)")

# affichage du soliton initial
lineplot, = plot(x,u0)
# lancement animation
ani = animation.FuncAnimation(fig1,KdV,frames=len(t),fargs=(x,y,lineplot), repeat=False)
tight_layout()
show()

Le coeur du script est la définition du système KdV qui sera intégré par odeint() :

def SystemeKdV(u,t,L):
    # calcul des dérivées partielles en x d'ordre 1 et 3
    ux = psdiff(u, period=L)
    uxxx = psdiff(u, period=L,order=3)
    # calcul de la dérivée temporelle d'ordre 1 (définition de l'EDP)
    ut = -6*u*ux - uxxx
    return ut

Chaque dérivée constituant l'équation (2) est calculée par la fonction psdiff(). En paramétres de cette fonction, on passe la fonction u à dériver et le domaine de dérivation L, considéré comme périodique.

Pour information, cette partie du script est fournie sur le site de scipy.

Le résultat

En lançant le script, nous obtenons l'animation suivante :

KdV Pseudo-spectrale

Comme vous le notez, le résultat obtenu est pratiquement identique à la solution analytique.

Résoudre KdV avec le schéma de Zabusky et Kruskal

Zabusky et Kruskal furent les premiers, en 1965, à proposer un schéma de résolution de l'équation KdV basé sur la méthode des différences finies. C'est un schéma de type leapfrog à trois niveaux, de premier ordre dans le domaine spatial et de second ordre explicite dans le domaine temporel. Il a été améloré par Taha et Ablowitz en 1984, pour donner ceci en partant de l'équation (2) :

\( u_i^{j+1} = u_i^{j-1} - 2\dfrac{\Delta t}{\Delta x}[u_{i+1}^{j}  + u_i^j + u_{i-1}^{j}][u_{i+1}^{j} - u_{i-1}^{j}] - \dfrac{\Delta t}{(\Delta x)^3} [u_{i+2}^{j} +2u_{i-1}^{j} -2u_{i+1}^{j} - u_{i-2}^{j}] \).

C'est un schéma très instable numériquement, qui nécessite un pas de temps \( \Delta t\) très petit. Sa condition de stabilité numérique, son CFL (voir la page sur les méthodes de DF) est  : \( \dfrac{\Delta t}{\Delta x} \le \dfrac{2}{\sqrt{3}} \dfrac{1}{6u_0 - \dfrac{3}{(\Delta x)^2}}\) où \( u_0\) est l'amplitude maximum du soliton.

Pour démarrer le calcul, il faut initialiser \(u_i^1 \) avec un soliton analytique aux conditions initiales et calculer \(u_i^2 \) avec l'expression :

\( u_i^2 = u_i^1 - 2\dfrac{\Delta t}{\Delta x}[u_{i+1}^1  + u_i^1 + u_{i-1}^1][u_{i+1}^1 - u_{i-1}^1] - \dfrac{\Delta t}{(\Delta x)^3} [u_{i+2}^1 +2u_{i-1}^1 -2u_{i+1}^1 - u_{i-2}^1] \)

Cet algorithme nécessite donc le calcul de deux conditions initiales et un CFL très couteux en temps. Si vous faites des expériences ou un TIPE sur les solitons, je vous recommande d'utiliser plutôt une méthode pseudo-spectrale que ce schéma qui est assez difficile à stabiliser.

A la rigueur, si vous tenez vraiment à utiliser une méthode DF, alors utilisez la méthode de Crank-Nicholson implicite de deuxième ordre. Elle est un peu plus simple à mettre en oeuvre car plus stable.

Néanmoins, pour ceux que l'aventure tenterait, je propose une implémentation du schéma ZK dans le script KdV_UnSolitonZK.py. En voici les parties essentielles :

  • la construction du soliton initial :
U[:,0] = KdVAnalytique(x,0,c1,shift1,mu,nu)
U[0,1] = 0. ; U[0,2] = 0.
U[Nx-1,1] = 0. ; U[Nx-1,2] = 0.
  • le calcul du terme U[:,1], nécessaire au démarrage du schéma :
for  i in range (1, Nx-1):                             
    if i > 1 and  i < Nx-2:
       a2 = U[i+2,0] + 2*U[i-1,0] - 2*U[i+1,0] - U[i-2,0]
    else:
       a2 = U[i-1, 0] - U[i+1, 0]
    U[i, 1] = U[i, 0] - \
              (mu/6.)*(dt/dx)*(U[i + 1, 0] + U[i, 0] + U[i - 1, 0])*(U[i+1, 0] - U[i-1, 0]) - \
              (nu/2.)*(dt/(dx**3))*a2
  • la fonction de calcul pour chaque pas :
def KdV(num):
    for  i in range (1, Nx-2):
         if i > 1 and i < Nx-2:
            a2 = U[i+2,1] + 2*U[i-1,1] - 2*U[i+1,1] - U[i-2,1]
         else:
            a2 = U[i-1, 1] - U[i+1, 1]
         U[i, 2] = U[i, 0] - \
                   A*(U[i + 1, 1] + U[i, 1] + U[i - 1, 1])*(U[i+1, 1] - U[i-1, 1]) - \
                   B*a2
    # sauvegarde des résultats
    line.set_data(x,U[:,2])
    # permutation des buffers    
    U[:, 0] = U[:, 1]
    U[:, 1] = U[:, 2]
    return line,

Voici le résultat de la simulation :

Simulation schéma ZK

Attention, comme je l'ai indiqué plus haut, le schéma est très instable. Aussi, si vous modifiez les paramètres des domaines temporel et spatial, vous risquez fort de constater des courbes très bizarres, voire de plantages du script. Le schéma ZK est difficile à "tuner" !

Collision de solitons

La collision de deux solitons est un phénomène très courant dans la nature, s'agissant par exemple de vagues dans un estuaire à marée montante. La simulation nous permet de reproduire très facilement ce phénomène.

Le script KdV_DeuxSolitonsAnalytique.py propose l'animation de la collision de deux solitons d'amplitude, et donc de vitesse, différente. Avant la collision, les deux solitons se présentent ainsi :

Collision solitons-1

Nous pouvons observer là aussi que le soliton d'amplitude la plus faible est plus étalé et moins rapide que le soliton d'amplitude la plus forte.

Voilà un aperçu de leur aspect pendant la collision :

Collision solitons-2

En observant la collision des deux solitons, vous noterez l'évolution de l'amplitude et de la largeur des deux solitons. L'équation KdV n'est pas linéaire, donc l'amplitude maximum constatée lors de la collision n'est pas égale à la somme des amplitudes des deux solitons. Il ne s'agit donc pas de la simple superposition de deux ondes, mais d'une interaction particulière.

Et après la collision :

Collision solitons-3

On observe qu'après la collision, les deux solitons retrouvent leur forme et leur célérité initiale, on parle de collision élastique, mais qu'ils ont subi un déphasage, que l'on voit mieux sur l'animation.

Influence des termes non-linéaire et dispersif

Une forme générale de l'équation KdV

Considérons la forme générale de l'équation (2) :

\( u_t + \mu uu_x + \nu u_{xxx} = 0 \)

dans laquelle \( \mu \) est le coefficient de non-linéarité, qui vaut 6 dans l'équation (2) et \( \nu \) le coefficient de dispersion, qui vaut 1 dans l'équation (2).

Dans la suite, j'utilise le script KdV_Etude.py en faisant varier les coefficients mu et nu. Ce script met en oeuvre la solution analytique de l'équation générale KdV.

L'influence du terme non-linéaire

Lançons le script pour les valeurs mu = 6 et nu = 1, cas nominal de l'équation. Nous obtenons :

mu=6, nu=1

Fixons maintenant la valeur de mu à 4. Nous obtenons :

mu=4, nu=1

Nous observons que l'amplitude du soliton est plus élevèe. Si mu = 2, nous obtenons :

mu=2, nu=1

L'augmentation d'amplitude se confirme.

Inversement, si vous augmentez le coefficient de non-linéarité, vous constaterez que l'amplitude du soliton diminue, sa largeur restant constante.

L'influence du terme dispersif

Conservons maintenant le coefficient non-linéaire mu = 6 et faisons varier le coefficient nu du terme dispersif. Fixons le à nu = 2, nous obtenons :

mu = 6 nu = 2

Nous observons que le soliton est plus aplati que lorsque le coefficient de dispersion nu est égal à 1. Fixons maintenant la valeur de nu = 3 :

mu = 6 nu = 3

L'aplatissement du soliton se confirme.Si maintenant, nous réduisons la valeur du coefficient de dispersion à 0.1, par exemple , nous obtenons :

mu = 6 nu = 0.1

Le soliton se resserre, la dispersion de la courbe est moindre.

Le soliton : un équilibre entre non-linéarité et dispersion

La propagation d'un soliton résulte d'un équilibre stable entre le poids respectif des coefficients de non-linéarité et de dispersion. Lorsqu'il s'agit de solitons réels, une vague ou un tsunami, ces coefficients dépendent de grandeurs physiques, en particulier du rapport entre la longueur d'onde de la vague et la profondeur du bassin.

Physiquement, quand le terme non-linéaire l'emporte dans l'équilibre, l'amplitude du soliton, de la vague, est grande. La partie haute de la vague va plus vite que la partie basse et la vague déferle. C'est ce qui se passe à l'approche de la plage, lorsque la profondeur devient très petite par rapport à la longueur d'onde de la vague.

Lorsque le terme dispersif l'emporte, terme qui est linéaire, c'est à dire quand la profondeur est grande, alors nous observons des petites vaguelettes, bien plates.

A l'équilibre, nous observons un soliton !

Contenu et design par Dominique Lefebvre - www.tangenteX.com novembre 2020    Licence Creative Commons    Contact :