Discrétisation d'équations différentielles

Avec Haskell

Voici de petites illustrations des méthodes d'Euler et du point milieu évoquées dans l'exercice 2.8 du poly

Haskell
--
-- MÉTHODES D'EULER
--
 
-- y' = f(x,y) avec un pliage
listeEuler f x0 y0 pas n =
    scanl (\ (x,y) k ->  (x +  pas, y + pas * (f x y) )) (x0,y0)  [0 .. n]
 
 
-- y'' = f(x,y,y') avec pliage 
listeEuler2 f x0 y0 y0' pas n  =
  map fst $ scanl
    (\ ((x,y),y') k -> ((x + pas, y + pas * y'),  y' + pas * (f x y y'))) ((x0,y0),y0') [0 .. n]
 
-- spés et signatures 
 ---- f : fonction définissant l'équa diff : Float -> Float -> Float
 ---- x0         : abs ini :: Float
 ---- y0         : ordonnée ini :: Float 
 ---- y0'        : dérivée ini :: Float
 ---- pas        : incrément des abscisses :: Float
 ---- n          : nombre de points créés :: Int(eger)
 ---- listeEuler : renvoie le liste des points de la courbe :: [(Float,Float)]
 
-- Remarque on peut créer un type synonyme par souci de clarté :
 ---- type Point = (Float,Float) alors listeEuler renvoie une liste de points : [Point]
 
-- tracé ordre 1
euler f x0 y0 xf n = 
  let pas = (xf - x0) / n in
    plot X11 --(PDF ("eulerExpo.pdf"))
             ( 
    Data2D [Title "Méthode d'Euler d'ordre 1", Color Blue, Style Lines] [] (listeEuler f x0 y0 pas n) 
    )
 
-- tracé ordre 2
euler2 f x0 y0 y0' xf n = 
  let pas = (xf - x0) / n in
    plot (PDF ("eulerExpoDeux.pdf"))
             ( 
    Data2D [Title "Méthode d'Euler d'ordre 2", Color Blue, Style Lines] [] (listeEuler2 f x0 y0 y0' pas n) 
    )

Par exemple, pour $y'=y \cos (4x)$:

Haskell
*AnalyseNumerique> euler1 (\ x y -> y * (cos (4*x)) ) 1 1 100 (2**10)

et pour $y''= - y$:

Haskell
*AnalyseNumerique> euler2 (\ x y y' -> - y) 0 0 1 10 (2**10)

Avec Python

Pour un rendu graphique d'une grande qualité, on utilise le module pygal de Python.

Ensuite, depuis ipython, on peut visualiser sur firefox par exemple, avec la commande !firefox ./courbeEulerPython.svg

Python
from pygal import *
 
def listeEuler(f,x0,y0,pas,n):
    x,y,L = x0,y0,[]
    for k in range(n):
        L += [(x,y)]
        x += pas
        y += pas * f(x,y)
    return L
 
def euler(f,x0,y0,xf,n):
    pas = (xf - x0) / n
    courbe  = XY()
    courbe.title = "Methode d Euler"
    courbe.add("Solution approchee",listeEuler(f,x0,y0,pas,n))
    courbe.render_to_file("courbeEulerPython.svg")
    !firefox ./courbeEulerPython.svg

Voici le magnifique rendu SVG pour $y'=y$, $y(0)=1$.



courtesy of webmatter.de