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