Projets TP méthodes numériques INFO1

Tant que vous tenterez de travailler sans papier ni crayon...

...La classe expr avec son polynôme de Taylor :

Python
from copy import copy
from math import exp, log, cos, sin
class expr:
 
    def __init__(self, valeur, gauche, droite) :
        self.val = valeur
        self.g = gauche
        self.d = droite
 
    def __add__(self, other) :
        if self.val == 0 :
            return other
        if other.val == 0 :
            return self
        if type(self.val) != str and type(other.val) != str :
            return expr(self.val + other.val, None, None)
        #if type(other.val) != str and other.val < 0 :
        #    return expr('-',self,-other)
        return expr('+', self, other)
    
 
    def __pow__(self,k) :
        if k.val == 1 :
            return self
        if k.val == 0 :
            return expr(1,None,None)
        return expr('**',self,k)
 
    def __truediv__(self,other) :
        if other.val == 1 :
            return self
        return self * (other ** expr(-1,None,None))
        
    def __mul__(self, other) :
        if type(self.val) != str and type(other.val) != str :
            return expr(self.val * other.val, None, None)
        if self.val == 0 or other.val == 0 :
            return expr(0,None,None)
        if self.val == 1. :
            return other
        if other.val == 1. :
            return self
        if self == other :
            return expr('**',self,expr(2,None,None))
        if self.val == '**' and other.val == '**' and self.g == other.g :
            return expr('**', self.g, self.d + other.d)        
        if type(self.val) != str and type(other.val) != str :
            return expr(self.val * other.val, None, None)
        return expr('*',self,other)
 
    
    def __neg__(self) :
        if type(self.val) == str :#self.val in {'*','+','**','cos','sin','ln','exp'} :
            return expr(-1,None,None) * self
        return expr(-self.val, None, None)
 
    def __sub__(self, other) :
        if self.val == 0 :
            return -other
        if other.val == 0 :
            return self
        if type(self.val) != str and type(other.val) != str :
            return expr(self.val - other.val, None, None)
        return expr('+',self,-other)
 
    def __repr__(self) :
        def par(v) :
            if v.val == None:
                return ''
            if v.g == None :
                return str(v)
            return '(' + str(v) + ')'
        if self.val in {'+', '-', '*', '/'} :
            return par(self.g) + self.val + par(self.d)
        if self.val in {'ln', 'cos', 'sin', 'exp'} :
            return self.val + '(' + str(self.g) + ')'
        if self.val == '**' :
            if self.d.val == 1 :
                return str(self.d)
            if self.d.val == 0.5 :
                return 'R(' + str(self.g) + ')'
            if type(self.d.val) == str or self.d.val > 0 :
                return par(self.g) + '^' + par(self.d)
            else :
                return '1/' + par(self.g ** (-self.d))
        return str(self.val)
 
    def eval(self, var, nb) :
        if self.val == '+' :
            return (self.g.eval(var,nb)) + (self.d.eval(var,nb))
        if self.val == '-' :
            return (self.g.eval(var,nb)) - (self.d.eval(var,nb))
        if self.val == '*' :
            return (self.g.eval(var,nb)) * (self.d.eval(var,nb))
        if self.val == '/' :
            return (self.g.eval(var,nb)) / (self.d.eval(var,nb))
        if self.val == '**' :
            return (self.g.eval(var,nb)) ** (self.d.eval(var,nb))
        if self.val == 'ln' :
            return log(self.g.eval(var,nb))
        if self.val == 'cos' :
            return cos(self.g.eval(var,nb))
        if self.val == 'sin' :
            return sin(self.g.eval(var,nb))
        if self.val == 'exp' :
            return exp(self.g.eval(var,nb))
        if self.val == var :
            return nb
        if type(self.val) == str :
            raise ValueError('Évaluation avec un paramètre impossible')
        return self.val
        
 
    def deriv(self, var) :
        if self.val == var :
            return expr(1,None,None)
        if self.val == '+' :
            return self.g.deriv(var) +  self.d.deriv(var)
        if self.val == '-' :
            return self.g.deriv(var) -  self.d.deriv(var)
        if self.val == '*' :
            return (self.g.deriv(var) * self.d) + (self.d.deriv(var) * self.g)
        if self.val == '/' :
            return ((self.g.deriv(var) * self.d) - (self.d.deriv(var) * self.g)
                    /(self.d ** expr(2,None,None)))                    
        if self.val == '**' :
            return self.d * self.g.deriv(var) * (self.g **(self.d - expr(1,None,None)))
        if self.val == 'ln' :
            return self.g.deriv(var) / self.g
        if self.val == 'exp' :
            return self.g.deriv(var) * self
        if self.val == 'cos' :
            return -self.g.deriv(var) *  expr('sin',self.g,None)
        if self.val == 'sin' :
            return self.g.deriv(var) *  expr('cos',self.g,None)
        return expr(0,None,None)
 
    def derivn(self,n,var) :
        if n == 0 :
            return self
        return (self.deriv(var)).derivn(n-1,var)
         
 
    def taylor(self,n,a) :
        f_derive_ordre_k = self # f^(k) c'est f au départ quand k = 0
        ix_moins_a  = Var('x') - Const(a) # x - a
        coeff = Const(1) # le coeff qui multiplie f^(k)(a) vaut 1 au départ
        Somme = Const(f_derive_ordre_k.eval('x',a))  # S contient f(a)
        for k in range(1, n + 1) :
            f_derive_ordre_k = f_derive_ordre_k.deriv('x') # on dérive f une fois de plus
            coeff = coeff * ix_moins_a / Const(k) #  coeff mult. par (x-a)/k
                                            # donne (x-a)^k / k!
            Somme += coeff * Const(f_derive_ordre_k.eval('x',a)) # on rajoute f^(k)(a) * coeff à la somme
        return Somme
 
 
 
    
# Des racourcis :            
 
def Ln(e) :
    return expr('ln',e,None)
 
def Cos(e) :
    return expr('cos',e,None)
 
def Sin(e) :
    return expr('sin',e,None)
 
def Exp(e) :
    return expr('exp',e,None)
 
def Const(k) :
    return expr(k, None, None)
 
def Var(x) :
    return expr(x, None, None)

Les lynx et les lapins :

Python
import matplotlib.pyplot as plt
import numpy as np
 
def euler1(f, x0, y0, xf, nb_pts):
    """
    Résolution de y' = f(y,x) par la méthode d'Euler explicite
    x0, y0 : position initiale
    xf : abscisse finale
    nb_pts : nb de points calculés
 
    renvoie un couple d'array numpy:
       - X : le tableau des abscisses xi = xi-1 + pas
       - Y : le tableau des ordonnées yi = yi-1 + pas*f(yi-1,xi-1)
 
    Exemple : résolution de y' = 3*y + x^2 avec P0=(0,1) et xf=5
        >>> euler(lambda y,x: 3*y + x^2, 0, 1, 5, 256)
        
    """
    h = (xf - x0) / (nb_pts - 1) # un intervalle de moins que le nb de piquets
    X = [x0 + k*h for k in range(nb_pts)]  # les abscisses subdivisées
    Y = [y0]*nb_pts # une liste d'ordonnées initialisée à y0
    for k in range(1, nb_pts): # on va calculer chaque yk en fonction du précédent
        Y[k] = Y[k-1] + h*f(Y[k-1], X[k-1])
    return (np.array(X),np.array(Y)) # on sort des array pour utiliser plt.plot ensuite
 
def euler2(f, x0, y0, yp0 ,xf, nb_pts):
    """
    Résolution de y'' = f(y',y,x) par la méthode d'Euler explicite
    x0, y0 : position initiale
    yp0 : dérivée initiale
    xf : abscisse finale
    nb_pts : nb de points calculés
 
    renvoie un couple d'array numpy:
       - X : le tableau des abscisses xi = xi-1 + pas
       - Y : le tableau des ordonnées yi = yi-1 + pas*f(yi-1,xi-1)
 
    Exemple : résolution de y'' = 0.05*sqrt(1 + h'**2) avec P0=(0,20), y'0=0 et xf=15
        >>> euler(lambda yp, y,x: 0.05*(1 + yp**2)**0.5, -15, 20, 0,15, 256)
        
    """
    h = (xf - x0) / (nb_pts - 1) # un intervalle de moins que le nb de piquets
    X = [x0 + k*h for k in range(nb_pts)]  # les abscisses subdivisées
    Y = [y0]*nb_pts # une liste d'ordonnées initialisée à y0
    Yp = [yp0]*nb_pts
    for k in range(1, nb_pts): # on va calculer chaque yk en fonction du précédent
        Yp[k] = Yp[k-1] + h*f(Yp[k-1], Y[k-1], X[k-1])
        Y[k]  = Y[k-1] + h*Yp[k-1]
    return (np.array(X),np.array(Y)) # on sort des array pour utiliser plt.plot ensuite
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
def trace_euler1(f,x0,y0,xf,nb_pts):
    X, Y = euler1(f,x0,y0,xf,nb_pts)
    plt.plot(X, Y)
    plt.show()
 
def trace_eulerV(f,x0,y0,xf,nb_pts) :
    T, X = euler1(f,x0,y0,xf,nb_pts)
    plt.show(T, X)
 
def LV(a,b,c,d) :
    return lambda X, t: np.array([X[0]*(a - b*X[1]), X[1]*(-c + d*X[0])])
    
def trace_euler2(f,x0,y0,yp0,xf,nb_pts):
    X, Y = euler2(f,x0,y0,yp0,xf,nb_pts)
    plt.plot(X, Y)
    plt.show()
 
    
def clouseau() :
    clock, temperature = euler(lambda temper, tailleme : -0.034*temper, 9, 30, 0, 1025)
    indice = 0
    while temperature[i] < 37 :
        indice += 1
    heure = clock[indice]
    h = int(heure)
    m = int((heure - h) * 60)
    print('Heure du crime -> ' + str(h) + ':' + str(m))
    return heure

courtesy of webmatter.de