Polynômes de Taylor : graphiques et animations

Comment construire les polynômes de Taylor à n'importe quel ordre de fonctions usuelles sans calcul formel sur Haskell ?

Pour construire $T_{n,0}(x)= \frac{x^0}{0!}f^{(0)}(0) + \frac{x^1}{1!}f^{(1)}(0) + \frac{x^2}{2!}f^{(2)}(0) + \cdots $ on va créer deux listes infinies : $[\frac{x^0}{0!}, \frac{x^1}{1!} , \frac{x^2}{2!} ,...]$ et $[f^{(0)}(0), f^{(1)}(0), f^{(2)}(0),....]$.


La première liste est facile à créer avec scanl qui renvoie la liste des différents états de l'accumulateur d'un pliage:

Haskell
-- Renvoie le vecteur [1, x, x^2/2!, x^3/3!, ..]
vecTaylor :: Double -> [Double]
vecTaylor    x      =  scanl (\ acc k -> acc * x / k) 1.0 [1 .. ]

Attention ! C'est une liste infinie...

Haskell
*AnalyseNumerique> take 5 (vecTaylor 1)
[1.0,1.0,0.5,0.16666666666666666,4.1666666666666664e-2]

Pour la deuxième liste, sans calcul formel, on va créer des listes pour chaque fonction élémentaire en cherchant une relation de récurrence.


Par exemple, pour l'exponentielle, c'est [1,1,1,1,....], pour $x \mapsto \frac{1}{1-x}$ c'est [1,1,2!,3!,4!,...], pour le sinus c'est [0,1,0,-1,0,1,0,-1,...].


On commence par donner des alias aux fonctions:

Haskell
-- Fonctions de base
-- fonction :: Double -> Double
invPlus    = (\ x -> 1 / (1 + x))
invMoins   = (\ x -> 1 / (1 - x))
logPlus    = (\ x -> log (1 + x))
racinePlus = (\ x -> sqrt (1 + x))

Ensuite, on crée un couple (fonction, liste de ses dérivées successives en 0) :

Haskell
{- On renvoie le couple (fonctions, liste de ses dérivées successives en 0)
   On cherche le plus souvent une formule de récurrence et on utilise scanl
   Sinon, on répète des motifs
tTruc :: (Double -> Double, [Double]) -}
tExp      = (exp,        repeat 1.0)
tInvMoins = (invMoins,   scanl (\ acc k -> acc * k   )               1.0 [1 .. ])
tInvPlus  = (invPlus,    scanl (\ acc k -> acc * (-k))               1.0 [1 .. ])
tRacine   = (racinePlus, scanl (\ acc k -> acc * (-0.5) * (2*k - 1)) 1.0 [0 .. ])
tLogPlus  = (logPlus,    0.0 : (snd tInvPlus))
tSin      = (sin,        cycle [0.0,1.0,0.0,-1.0])
tCos      = (cos,        cycle [1.0,0.0,-1.0,0.0])

Ce sont encore des listes infinies...


Par exemple pour $x \mapsto \frac{1}{1-x}$:

Haskell
*AnalyseNumerique> take 10 (snd tInvMoins)
[1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0]

Il ne reste plus qu'à effectuer le produit scalaire des deux listes infinies avec sum et zipWith :

Haskell
{- Renvoie l'évaluation du poly de Taylor de degré n en x de la fonction
   dont on donne la liste des dérivées en 0 -}
polyTaylor :: [Double]   -> Int -> Double -> Double 
polyTaylor    listeDerivees n      x      =
    sum $ take (n + 1) $ zipWith (*) (vecTaylor x) listeDerivees 

Pour l'afficher, on peut utiliser la simplissime bibliothèque EasyPlot qui permet d'interfacer GNUPlot sous Haskell.

On n'oubliera pas d'importer EasyPlot

Haskell
import Graphics.EasyPlot

On peut alors utiliser plot X11 (Function2D [Title ..., Color ... ] [Range ... ..., Step ...] fonction) :

Haskell
-- On trace sur un même graphique la fonction et son poly de Taylor sur [a,b] avec un pas de 1/2^8
compareDL :: Double -> Double -> ((Double -> Double) , [Double]) -> Int -> IO Bool
compareDL    a         b         tf                                 n   = 
  plot  X11 
     [ 
    Function2D [Title ("n = " ++ (show n))   , Color Blue ] [Range a b, Step (2**(-8))]
                   (polyTaylor (snd tf) n),
    Function2D [Title "Vraie fonction", Color Red  ] [Range a b, Step (2**(-8))]
                   (fst tf)
     ] 

Ensuite, on demande l'affichage dans un terminal GNUplot avec l'option X11 (ou Windows ou Aqua selon les OS) où un export en PDF, PNG, JPG, etc.

Haskell
*AnalyseNumerique> compareDL (-0.8) 1.2 tLogPlus 4

Renvoie, pour ${\rm ln}(1+x)$:



On peut en profiter pour créer une petite animation. Pour rester libre, on évitera Flash, GIF et cie....


Voici une petite animation PNG. On exporte les images en PNG et on utilise par exemple apngasm.


On commence par modifier compareDL pour pallier au problème de l'absence de contrôle sur l'échelle des y:

Haskell
animIniDL :: Double -> Double -> ((Double -> Double) , [Double]) -> Int -> IO Bool
animIniDL    a         b         tf                                 n   = 
  plot  (PNG ("animIndi" ++ (show (n + 10)) ++ ".png")) (
      [
       Function2D [Title " "   , Color White ] [Range a b, Step (2**(-8))] -- pour pallier l'abscence de contrôle sur l'échelle des y
                      (polyTaylor (snd tf) k) | k <- [0..10]
      ] ++ [ 
    Function2D [Title ("n = " ++ (show n))   , Color Blue ] [Range a b, Step (2**(-8))]
                  (polyTaylor (snd tf) n),
    Function2D [Title "Vraie fonction", Color Red  ] [Range a b, Step (2**(-8))]
                  (fst tf)
           ] 
            )

Ensuite on crée une liste d'objets graphiques :

Haskell
animDL :: Double -> Double -> (Double -> Double, [Double]) -> Int -> [IO Bool]
animDL a b tf n = [animIniDL a b tf k | k <- [0 .. n] ]

Haskell doit communiquer avec le monde extérieur via ici la monade IO. On utilise l'intruction do qui permet d'exécuter des commandes IO séquentiellement : cela donne l'illusion de faire de la programmation impérative...

Haskell
suite :: [IO Bool] -> IO ()
suite []     = return ()
suite (x:xs) = do x
                  suite xs

Il ne reste plus qu'à créer une fonction qui crée la suite de figures de manière transparente :

Haskell
anim :: Double -> Double -> (Double -> Double, [Double]) -> Int -> IO ()
anim a b tf n = suite (animDL a b tf n)

Par exemple pour les 11 premiers polynômes de Taylor de $x\mapsto {\rm ln} (1+x)$ sur $[-0.8,1.3]$:

Haskell
anim (-0.8) 1.3 tLogPlus 10

Il reste, dans un shell, à regrouper les fichiers PNG

Bash
$ apngasm animLog.png animIndi*.png 1 1

pour créer cette animation:



courtesy of webmatter.de