Un dérivateur formel en (presque) 5 minutes

Un dérivateur numérique, c'est dangereux alors qu'un dérivateur formel, c'est si facile à écrire...en Haskell bien sûr !

(Par curiosité, vous pouvez allez jeter un coup d'œil sur la version Python ou la version OCAML de ce dérivateur...)

Haskell
module DerivFormel where
 
data Fonction =
    Const Float
        | Param String
        | Var String
        | Add Fonction Fonction
        | Mul Fonction Fonction
        | Ln Fonction
        | Exp Fonction
        | Puis Fonction Fonction
        | Cos Fonction
        | Sin Fonction
 
 
infix 9 @^
infix 5 @*
infix 3 @+
    
-- les fonctions arithmétiques avec leurs simplifications             
(@*) :: Fonction -> Fonction -> Fonction
(@*) (Const k) (Const k') = Const (k*k')
(@*) (Const k) (Mul (Const k') f) = (Const (k*k')) @* f
(@*) (Const k) (Mul f (Const k')) = (Const (k*k')) @* f
(@*) (Mul (Const k') f) (Const k) = (Const (k*k')) @* f
(@*) (Mul f (Const k')) (Const k) = (Const (k*k')) @* f
(@*) (Const 0) _ = Const 0
(@*) _ (Const 0) = Const 0
(@*) f (Const 1) = f
(@*) (Const 1) f = f
(@*) f g = Mul f g
 
 
(@+) :: Fonction -> Fonction -> Fonction
(@+) (Const k) (Const k') = Const (k + k')
(@+) (Add (Const k') f) (Const k) = (Const (k+k')) @+ f
(@+) (Add f (Const k')) (Const k) = (Const (k+k')) @+ f
(@+) f (Const 0) = f
(@+) (Const 0) f = f
(@+) f g = Add f g
 
(@^) :: Fonction -> Fonction -> Fonction
(@^) _ (Const 0) = Const 1
(@^) f (Const 1) = f
(@^) f n = Puis n f 
 
-- L'opérateur de dérivation formelle                  
deriv :: Fonction -> Fonction -> Fonction
deriv (Const _) _  = Const 0
deriv (Param _) _  = Const 0
deriv (Add f g) x  = (deriv f x)  @+ (deriv g x)
deriv (Mul f g) x  = ((deriv f x) @* g) @+ ((deriv g x) @* f)
deriv (Puis n f) x = (deriv f x) @* (n @* (f @^ (n @+ (Const (-1)))))
deriv (Ln f) x     = (deriv f x)  @* (f @^ (Const (-1)))
deriv (Exp f) x    = (deriv f x)  @* (Exp  f)
deriv (Cos f) x    = ((Const (-1)) @* (deriv f x)) @* (Sin f)
deriv (Sin f) x    = (deriv f x)  @* (Cos f)
deriv (Var y) (Var x) = if x == y then Const 1 else Const 0
 
-- Pour le calcul de la dérivée n-eme
derivn :: Fonction -> Int -> Fonction -> Fonction
derivn f n x = foldl (\ facc _ -> deriv facc x) f [1..n]
 
-- Pour évaluer une fonction en 1 flottant
eval :: Fonction -> Float -> Float
eval (Var _) t    = t
eval (Const k) _  = k
eval (Mul f g) t  = (eval f t) * (eval g t)
eval (Add f g) t  = (eval f t) + (eval g t)
eval (Puis n f) t = (eval f t)**(eval n t)
eval (Ln f) t     = log (eval f t)
eval (Cos f) t    = cos (eval f t)
eval (Sin f) t    = sin (eval f t)
eval (Exp f) t    = exp (eval f t)
 
-- pour un pas trop horrible affichage au format LateX
montre :: Fonction -> String
montre (Const x)            = show x
montre (Var x)              =  x
montre (Param k)       = k  
montre (Add f (Const 0))    = montre f
montre (Add (Const 0) f)    = montre f
montre (Mul (Const 0) _)    = show (0)
montre (Mul _ (Const 0))    = show (0)
montre (Mul (Const 1) f)    = montre f
montre (Mul f (Const 1))    = montre f
montre (Mul (Const ((-1))) f) = "-" ++ montre f
montre (Mul f (Const ((-1)))) = "-" ++ montre f
montre (Add f g)            = "(" ++ (montre f)  ++ ") + (" ++ (montre g) ++ ")"
montre (Mul f g)            = "(" ++ (montre f) ++  ") * (" ++ (montre g) ++ ")"
montre (Ln f)               = "\\ln(" ++ (montre f) ++ ")"
montre (Exp f)              = "\\exp(" ++ (montre f) ++ ")"
montre (Cos f)              = "\\cos(" ++ (montre f) ++ ")"
montre (Sin f)              = "\\sin(" ++ (montre f) ++ ")"
montre (Puis (Const k) f) =
    if k == 1 then montre f
    else
        if k == 0.5 then "\\sqrt{" ++ (montre f) ++ "}"
        else
            if k > 0 then "(" ++ (montre f) ++ ")^{" ++ (show k) ++"}"
            else "\\frac{1}{" ++ (montre  (Puis (Const (-k)) f)) ++ "}"
montre (Puis k f) = "(" ++ (montre f) ++ ")^{" ++ (montre k) ++ "}"
 
instance  Show Fonction where
    show = montre

Par exemple, pour dériver $x\mapsto \ln(\ln(\ln(x)))$:

Haskell
λ> let x = Var "x"
λ> let f = Ln (Ln (Ln x))
λ> f
\ln(\ln(\ln(x)))
λ> deriv f x
((\frac{1}{x}) * (\frac{1}{\ln(x)})) * (\frac{1}{\ln(\ln(x))})

Ce qui se traduit en copiant-collant cette dernière ligne dans ce navigateur : $$((\frac{1}{x}) * (\frac{1}{\ln(x)})) * (\frac{1}{\ln(\ln(x))})$$

et pour la dérivée seconde:

Haskell
λ> derivn f 2 x
((((-\frac{1}{(x)^{2.0}}) * (\frac{1}{\ln(x)})) + (((\frac{1}{x}) * (-\frac{1}{(\ln(x))^{2.0}})) * (\frac{1}{x}))) * (\frac{1}{\ln(\ln(x))})) + ((((\frac{1}{x}) * (\frac{1}{\ln(x)})) * (-\frac{1}{(\ln(\ln(x)))^{2.0}})) * ((\frac{1}{x}) * (\frac{1}{\ln(x)})))

qui se lit : $$((((-\frac{1}{(x)^{2.0}}) * (\frac{1}{\ln(x)})) + (((\frac{1}{x}) * (-\frac{1}{(\ln(x))^{2.0}})) * (\frac{1}{x}))) * (\frac{1}{\ln(\ln(x))})) + ((((\frac{1}{x}) * (\frac{1}{\ln(x)})) * (-\frac{1}{(\ln(\ln(x)))^{2.0}})) * ((\frac{1}{x}) * (\frac{1}{\ln(x)})))$$

...c'est encore perfectible mais c'est déjà pas mal pour 5 minutes de boulot...

Autre exemple: $x\mapsto \frac{1}{1+x}$.

Haskell
λ> let g = ((Const 1) @+ x) @^ (Const (-1))
λ> derivn g 5 x
(-120.0) * (\frac{1}{((1.0) + (x))^{6.0}})

c'est-à-dire $$(-120.0) * (\frac{1}{((1.0) + (x))^{6.0}})$$

Avec une racine carrée: $x\mapsto \sqrt{1+x}$.

Haskell
λ> let h = ((Const 1) @+ x) @^ (Const (-0.5))
λ> deriv h x
0.5 * (\frac{1}{\sqrt{1.0 + (x)}})

c'est-à-dire $$0.5 * (\frac{1}{\sqrt{1.0 + (x)}})$$

Avec un paramètre: $x\mapsto (x+1)^k$

Haskell
λ> let x = Var "x"
λ> let k = Param "k"
λ> let g = ((Const 1) @+ x) @^ k
λ> derivn g 3  x
((((-2.0) + (k)) * (((1.0) + (x))^{(-3.0) + (k)})) * ((k) + (-1.0))) * (k)

c'est-à-dire $$g^{(3)}(x)=((((-2.0) + (k)) * (((1.0) + (x))^{(-3.0) + (k)})) * ((k) + (-1.0))) * (k)$$

Avantage du calcul différentiel formel sur le calcul numérique

On voudrait évaluer la dérivée 5ème du cosinus en 0.

En formel, aucune surprise:

Haskell
λ> let x = Var "x"
λ> let f = Cos x
λ> eval (derivn f 5 x) 0
-0.0
λ> derivn f 5 x
-\sin(x)

En numérique, définissons rapidement les taux de variations:

Haskell
-- Dérivation numérique
deriveNum :: (Float -> Float) -> Float -> (Float -> Float)
deriveNum f h = \ t ->  (f(t + h) - f(t)) / h
 
derivenemeNum :: (Float -> Float) -> Float -> Int ->  (Float -> Float)
derivenemeNum f _ 0 = f
derivenemeNum f h n = derivenemeNum (deriveNum f h) h (n - 1) 

Alors, pour la dérivée cinquième de cos (i.e. -sin) en 0, on obtient:

Haskell
λ> let gnum = derivenemeNum cos 1e-3 5
λ> gnum 0
5.960451e8

$5,96\times 10^{8}$, c'est un peu beaucoup pour un sinus...