2- On construit alors notre algorithme nous permettant de calculer g1(i,j) en fonction de g(i,j) :

do i=2,imax-1             do j=2,jmax-1 c On calcule les dérivées partielles :                 gxx=(g(i+1,j)-2*g(i,j)+g(i-1,j))/(dx**2)               gyy=(g(i,j+1)-2*g(i,j)+g(i,j-1))/(dy**2)               t1=g(i+1,j+1)-g(i+1,j-1)-g(i-1,j+1)+g(i-1,j-1)               gxy=(t1)/(4*dx*dy)   c On peut alors calculer g1(i,j) en fonction de g(i,j) :              g1(i,j)=g(i,j)+dt*(f(i,j)-gxx*gyy+(gxy)**2)             enddo          enddo

3- Il faut maintenant modéliser la condition d'arrêt :

 Pour cela On doit se fixer une norme pour , on prendra la norme discrétisée suivante :

 Dans ces conditions , on décidera d'arrêter notre algorithme lorsque , où sera fixé au début de notre programme .

  Nous pouvons alors programmer le calcul de cette norme ainsi que la condition d'arrêt :

C Calcul du résidu :          temp=0          do i=1,imax             do j=1,jmax             temp=temp+abs(g(i,j)-g1(i,j))             enddo          enddo C on normalise le résidu :       if (kt.eq.1) reso=temp se sera la valeur du 1° résidu       if (reso.eq.0) goto 10       res=temp/reso        c Condition d'arrêt :             if (res.lt.1e-4) goto 10

 

 Il se peut que pour des valeurs de dx , dy et dt données la valeur ne tende pas vers 0 , et même 'explose' . Dans ces conditions nous limiterons le nombre de calcul ( la variable Ktmax représentera ce nombre ) et si nous observons que le résidu explose nous pourrons aussi choisir d'arrêter le programme .

 Nous allons maintenant tester notre programme en nous fixant certaines valeurs .

3) Premiers tests du programme :

 Dans un premier temps pour vérifier si l'algorithme créée est correct nous allons nous donner une fonction solution , puis nous pouvons à partir de ses dérivées partielles construire un   vérifiant :

 Nous discrétiserons alors ce f et nous regarderons si la fonction solution donnée par notre programme est proche de celle donnée par .

Prenons pour valeur de  la fonction :

Cette fonction est bien nulle sur le bord de . Après calcul nous en déduisons notre fonction  :

Avec le maillage precedent définit nous avons le graphe de u :

 Fig 1 : Graphe de u

 Nous allons donc maintenant faire tourner notre programme . Pour cela nous allons nous donner une fonction initiale nulle sur . Nous allons d'abord partir de la solution , c'est à dire : , si notre algorithme est correct , le programme devrait s'arrêter à la première itération , et cela quel que soit le pas de temps choisit . Hors ce n'est pas le cas . Non seulement le programme ne s'arrête pas mais en plus ,au bout d'un certain nombre d'itération , le résidu 'explose' .

 Par exemple à l'itération 69181 le résidu a une valeur proche de 10000 et la fonction donnée prend cette forme :

Fig 2 : graphe obtenue

Fig 3 : Evolution du résidu

 Le résultat est éloigné de la fonction initiale . Nous voyons apparaître deux 'pics' . Si nous ne nous préoccupons pas de ces 'pics' , en les tronquant nous remarquons que notre résultat est proche de la fonction solution .

Fig 4: graphe n°3 obtenue sans les deux 'pics'

 En fait ce premier problème peut être résolu en prenant dans notre programme pour la fonction , non plus les valeurs exactes données par :

 Mais par les valeurs approchées données par :

 Nous n'avons alors plus le problème lié à l'approximation des dérivées , qui faisait que notre programme ne s'arrêtait pas à la première itération . Nous allons donc calculer directement à partir des valeurs de la fonction .

c On définit notre u(i,j) do i=1,imax          do j=1,jmax       u(i,j)=sin(pi*(-1+2*((i-1)*dx)))*sin(pi*(-1+2*((j-1)*dy)))        enddo       enddo c on calcule alors f(i,j) à partir de u :       do i=2,imax-1        do j=2,jmax-1                 uxx=(u(i+1,j)-2*u(i,j)+u(i-1,j))/(dx**2)        uyy=(u(i,j+1)-2*u(i,j)+u(i,j-1))/(dy**2)        t1=u(i+1,j+1)-u(i+1,j-1)-u(i-1,j+1)+u(i-1,j-1)        uxy=(t1)/(4*dx*dy)        f(i,j)=uxx*uyy-(uxy)**2        enddo       enddo

Nous pouvons alors relancer notre programme et nous apercevoir que cette fois-ci il s'arrête bien à la première itération et nous donne la fonction solution sans y avoir touché ( Le résidu trouvé est nul , il correspond à la variation entre la fonction au temps 0 , c'est à dire , et celle au temps 1 , c'est à dire notre solution ) .

 Nous allons maintenant essayer de partir d'une fonction proche de la solution et regarder si notre programme retrouve la solution .

 Pour cela nous allons considérer la fonction de départ :

 Nous avons donc perturbé notre fonction par une petite variation :  sur  et 0 sur :

Fig 5: Fonction perturbation

 Nous allons nous fixer un pas de temps . Si nous faisons tourner notre programme , nous-nous rendons compte que notre résidu tourne autour de la valeur 1 et puis à partir d'un certain temps explose et le résultat obtenu n'est pas en accord avec la solution souhaitée .

Fig 6 : résultat à t=15172 , résidu =

Fig 7 : Evolution du résidu

  La petite perturbation aux autres points du maillage ne s'atténue pas :

Fig 8: Différence entre la solution théorique et la solution donnée par le programme

Nous voyons aussi apparaître dans le graphe n°6 un 'pic' , correspondant à l'instabilité de notre schéma en ces points . Nous voyons donc ici  se dessiner le principal problème rencontré lors de l'utilisation d'un schéma explicite , à savoir l'instabilité du schéma .En effet dans la pratique quel que soit le pas de temps fixe pris notre résidu explosera tôt ou tard et le résultat trouvé ne sera pas correct .

 Si à l'itération 151482 le résidu avoisinera 1000646

 Si à l'itération 159 le résidu prendra une valeur proche 770882 .

Fig 9 : résultat pour

 Les points d'instabilité restent les mêmes . La seule chose qui change c'est la vitesse à laquelle notre schéma diverge .Nous allons donc devoir exhiber une condition de stabilité qui , nous l'espérons sera suffisante .

4) Condition de stabilité :

 La principale difficulté vient du fait que notre E.D.P n'est pas linéaire ! En effet nous ne pouvons directement utiliser l'analyse de Fourier . Nous allons donc pour cela linéariser notre problème :

 Le schéma explicite reviens à définir une fonction , par :

 Le problème vient du fait qu'à chaque itération notre fonction est une valeur approchée , on a des erreurs d'arrondis . En fait on trouve une fonction proche du  exact : . On voudrait s'assurer que notre schéma ne s'éloigne pas trop de notre solution . c'est à dire , si l'on note  le résultat trouvé :

 pour  proche de zéro .

 Or :

 avec

c'est à dire :

  Donc on va vouloir s'assurer que notre ( que l'on notera pour alléger les notations : ) soit stable . 

 On va calculer au point du maillage i ,j . En fait la fonction , ne dépend que des paramètres : avec .

 Dans ces conditions , on a alors :

On peut alors calculer  :

C'est à dire pour un point du maillage i,j :

                                                (3)

Grâce à la linéarité de peut maintenant utiliser l'analyse de Fourier pour dégager une condition de stabilité pour .

 On pose :

 avec et

Alors en remplaçant dans la formule (3) , et en prenant  ( on prend un maillage uniforme ) on trouve :

 Alors on doit avoir la condition nécessaire de stabilité de Von Neumman ,c'est à dire :

On en déduit , alors la condition de Stabilité :

                (4)

 Nous allons maintenant implanter cette condition de stabilité dans notre programme et voir si elle nous permet de retrouver dans le cas :   la bonne solution . On va traduire la condition (4) par :

lorsque                   (5)

 Dans le programme , cela donne les lignes de code suivantes :

     
       cond=2*(abs(gxy))+4*(abs(gxx)+abs(gyy))

       dtt(i,j)=dt                 au cas où cond soit nul on fixe dtt(i,j)

       if (cond.gt.0) dtt(i,j)=(dx*dy)/cond
 

 On peut tracer le graphe des valeurs de vérifiant la condition (5) pour .

Fig 10 : Graphe de   pour .

 On fait maintenant tourner notre programme avec la condition de stabilité . On part de la fonction perturbée : . Mais notre programme ne converge toujours pas . Le résidu explose rapidement : à la 59 itération, il prend une valeur proche de 121610 . Le graphe obtenu n'est pas non plus satisfaisant :

Fig 11 : résultat à t=15172 , résidu =121610

Fig 12 : Evolution du résidu

 Les points d'instabilité sont plus nombreux qu'en prenant un pas de temps fixe . Et lorsque l'on tronque les nombreux pics , le reste du graphe n'est guère plus ressemblant à notre fonction solution :

Fig 13 : Graphe n°10 sans les 'pics'

 En fait la condition que l'on vient de trouver n'est pas suffisante pour assurer la stabilité de notre schéma . Nous voyons donc bien que dans le cas non linéaire , exhiber une condition de stabilité est difficile . Le schéma explicite n'est donc pas le plus approprié . Nous allons néanmoins essayer de voir des cas ou notre programme ne diverge pas trop rapidement afin de voir quels types de résultat nous pourrions obtenir .

5) Un cas de convergence du programme :

 On prend un pas de temps constant : .

Nous allons partir de fonctions telles que :

. si

. sinon ( où c est une constante réelle )  

Fig 14 : Exemple de fonction pour c=2

  1° cas : c=2 :

 Le résidu tend vers 0 pendant un certain nombre d'itération , puis se met à diverger . Par exemple à l'itération 15745 le résidu est de l'ordre de 0.08 .  Mais à l'itération 169442 il explose et prend la valeur 134 , après il tend vers l'infini .

 

Fig 15 : Evolution du résidu

Il est intéressant de regarder la fonction donnée par notre programme juste avant que le résidu n'explose :

 

 Fig 16 : Graphe donnée à l'itération 15745

Fig 17 : Evolution du résidu à l'échelle logarithmique

  On voit donc que notre résidu décroît de plus en plus rapidement , et notre fonction donnée par le programme tend vers une fonction que nous pouvons qualifier de 'cake' à 4 bosses . Notons que si l'on diminue ou l'on augmente le pas de temps on ne change en rien le résultat , mais le résidu se met à diverger plus ou moins vite . Par exemple :

 -Pour à l'itération 1695 le résidu sera proche de 342

 -Pour à l'itération 169227 le résidu sera proche de 100

De plus la condition de stabilité trouvée précédemment ne nous d'aucune utilité , le schéma divergeant dés la seconde itération .

 

Fig 18 : A droite le graphe obtenue pour et le résidu proche de 100 . A gauche le graphe obtenue avec la condition de stabilité à l'itération 59 .

 Nous voyons donc apparaître des points d'instabilité pour un pas de temps fixe . Et nous noterons aussi que la condition de stabilité trouvée nous donne un résultat aberrant en faisant exploser le nombre de point d'instabilité !

 

2° cas : c=0 :

 De la même manière que pour le cas c=2 , le résidu converge vers 0 , avant d'exploser . Le résultat trouvé avant la divergence du résidu , nous montre aussi la formation de 4 bosses , on peut aussi noter l'apparition de 4 points d'instabilité :

 



Fig 19 : Graphe donnée à l'itération 3063

2° cas : c= :

On a toujours le même comportement du résidu . On remarque que la fonction donnée est de la même forme que dans le cas c=2 :

 

Fig n°20 :Graphe donnée avant divergence du résidu

 

  Il semble donc que notre programme tende vers une solution différente de . En fait il semblerait que cette fonction en forme de 'cake' à 4 bosses soit un minimum local pour la fonction d'énergie des déplacements et non le minimum global solution de notre problème .

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6) Autre approche possible :

 Nous voyons donc que la méthode explicite n'est pas appropriée à la situation . La principale difficulté est due au fait que l'E.D.P est non linéaire et qu'il est ainsi difficile d'extraire des conditions de stabilité .

 Il sera donc nécessaire pour résoudre ce problème d'explorer d'autres voies , notamment celle des méthodes de gradient . En effet ces dernières s'applique bien au cas non linéaire .

Par exemple :

  On considère la fonction énergie :

 On se donne une fonction de départ . et on calcule de manière numérique :

où i représente chaque point du maillage . Par exemple avec le maillage uniforme choisis précédemment on a .

 On calcule alors J dans la direction du gradient :

i.e : On prend pour a réel et on calcule .

 On cherche la valeur de pour laquelle est minimale .

 A titre culturel , Le professeur Y. Brenier a résolue le problème de Monge-Kantorovich par une méthode de Gradient . On remplace le problème par un nouveau problème où les inconnues dépendent d'une variable supplémentaire assimilée au temps. Le problème devient un problème de minimisation convexe sous contraintes linéaires. Il est équivalent à un problème de point selle et peut être résolu par une méthode de gradient de type Uzawa.

 7) Annexe : 

Programme :

  program MongeKantorovich 


c*****************************
c* Définition des paramètres *
c*****************************
 
      parameter(imax=30,jmax=30,ktmax=1000000)

      integer i,j,kt

      real ctrl 

      double precision g(imax,jmax),g1(imax,jmax),f(imax,jmax)

      double precision t1,temp,reso,u(imax+1,jmax+1),xa,xb,toto,nj

      double precision res,gxx,gxy,gyy,dx,dy,pi,dt,cond,epsi,aj(1000)

      double precision resl,g2(imax,jmax),dtt(imax,jmax)

     

 
c Définition du maillage et d'un pas de temps fixe :


      dx=1./(imax-1)

      dy=1./(jmax-1)

      dt=1e-8
 
      pi=3.1415926535 Valeur approchée du nombre 
      tf=0

      res=2         Valeur initiale du résidu ( pour le test résidu 1 )   

      ctrl=1000     On affiche le résidu toute les CTRL itérations 

      condition=0   condition de stabilité (0=inactive)

      condt1=0      1° test sur le résidu (0=test inactif)

      condt2=1.5    Valeur du 2° Condition de test sur le résidu

 
c***************************************   
c* Définition des fonctions de départs *
c***************************************


c On Définit le second membre f :

      do i=1,imax

         do j=1,jmax

           g(i,j)=0

         enddo

      enddo

      do i=2,imax-1

         do j=2,jmax-1

      u(i,j)=sin(pi*(-1+2*((i-1)*dx)))*sin(pi*(-1+2*((j-1)*dy))) 

         enddo

      enddo
      do i=2,imax-1

       do j=2,jmax-1

             

       gxx=(u(i+1,j)-2*u(i,j)+u(i-1,j))/(dx**2)

       gyy=(u(i,j+1)-2*u(i,j)+u(i,j-1))/(dy**2)

       t1=u(i+1,j+1)-u(i+1,j-1)-u(i-1,j+1)+u(i-1,j-1)

       gxy=(t1)/(4*dx*dy)

       f(i,j)=gxx*gyy-(gxy)**2
       dtt(i,j)=dt

       

       enddo

      enddo
 
c on donne les valeurs au Bord de  :


      do i=1,imax

         g(i,1)=0

         g1(i,1)=0

         g(i,jmax)=0

         g1(i,jmax)=0

      enddo
      do j=1,jmax

         g(1,j)=0

         g1(1,j)=0

         g(imax,j)=0

         g1(imax,j)=0

      enddo


c on définie les valeurs de la fonction initiale v dans  :
 
      do j=2,jmax-1

         do i=2,imax-1

        g(i,j)=........   On donne la valeur désirée de v ici

        g2(i,j)=g(i,j)

         enddo

      enddo


c**********
c* Boucle *
c**********
 
      open(1,file='résidu',type='unknown')
   

      do kt=1,ktmax
 14      continue
         do i=2,imax-1

            do j=2,jmax-1
 15            continue
              gxx=(g(i+1,j)-2*g(i,j)+g(i-1,j))/(dx**2)

              gyy=(g(i,j+1)-2*g(i,j)+g(i,j-1))/(dy**2)

              t1=g(i+1,j+1)-g(i+1,j-1)-g(i-1,j+1)+g(i-1,j-1)

              gxy=(t1)/(4*dx*dy)
 
c Condition de stabilité :
 
              if (condition.eq.0) goto 501
              cond=2*(abs(gxy))+4*(abs(gxx)+abs(gyy))

              dtt(i,j)=dt

              if (cond.gt.0) dtt(i,j)=(dx*dy)/cond
 501          continue

       

            g1(i,j)=g(i,j)+dtt(i,j)*(f(i,j)-gxx*gyy+(gxy)**2)

           enddo

         enddo
 
c******************
c* test du résidu *
c******************
 
c Calcul du résidu et normalisation :
 
         resl=res
         temp=0

         do i=1,imax

            do j=1,jmax

            temp=temp+abs(g(i,j)-g1(i,j)) 

            enddo

         enddo   

      res=temp

      if (kt.eq.1) reso=temp

      if (reso.eq.0) goto 10    Si le 1° résidu est nul on arrête le programme

      tf=tf+1

      res=temp/reso             Normalisation du résidu
 
C On affiche le résidu à chaque CTRL itération :
 
      if (tf.eq.ctrl) write(*,*) ' res=',res,' kt=',ktmax-kt

            write(1,*) kt,res

      if (tf.eq.ctrl) tf=0        

 
c contrôle des valeurs du résidu :
 
      if (condt1.eq.0) goto 503
 
      if (resl.lt.res) goto 16
 
 503  continue

      if (condt2.lt.res) goto 16

         do i=2,imax-1

           do j=2,jmax-1

            g(i,j)=g1(i,j)

           enddo

         enddo
      if (res.lt.1e-4) goto 10
 
      enddo
 
c Le résidu est considéré comme invalide par les tests :
 
 16   continue
c RESL correspond au résidu à l'itération précédente
      write(*,*) 'Résidu incorrect res=',res,' resl=',resl 
 
 10   continue

 
c**************  

c* sauvegarde *
c**************
 
  Cette partie concerne la sauvegarde des divers résultats , 
elle n'est donc pas retranscrite car peu intéressante .
 

 Contact et liens internet :

e-mail : [email protected]