trapezesSoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for DUMMIES 19-20
# Problème 3
#
# Solution détaillée d'un problème plutot simple, quoique !
#  Vincent Legat
#
# -------------------------------------------------------------------------
# 

from numpy import *

#
# -1- Implémentation de la méthode des trapèzes !
#     Ca, c'est vraiment super simple, les gaminets :-)
#     Fallait vraiment les prendre les points associés à cela : y aura pas une seconde chance !
#

def trapezeEasy(f,a,b,n):
  X = linspace(a,b,n+1)
  U = f(X)
  h = (b-a)/n
  return (sum(U) - U[0]/2 - U[-1]/2)*h

#
#     Certains petits futés ont découvert qu'il existait une fonction numpy
#     qui faisait le job pour vous : numpy.trapz(f(X),X) : il était donc possible
#     d'écrire deux uniques lignes de code pour implémenter cela et de prendre la tangente !
#     Oui, oui : c'est une solution parfaitement admise mais cela ne rapporte rien de plus :-)
#
#     def trapezeEasy(f,a,b,n):
#       X = linspace(a,b,n+1)
#       return trapz(f(X),X)
# 

#
# -2- Implémentation de la division successive des intervalles
#     Si vous avez juste appliqué les consignes, il suffit d'écrire un code comme suit.
#     C'était pas bien compliqué : on pourrait admettre comme test d'arrêt n < nmax ou 2*n < nmax :-)
#     Comme expliqué mercredi, je reconnais une certaine ambiguité dans l'énoncé à ce sujet.
#
 
def trapezeEasyFun(f,a,b,n,nmax,tol):
  error = tol+1
  I = trapezeEasy(f,a,b,n)
  while (error > tol and n*2 < nmax):
    Iold = I
    n = n*2
    I = trapezeEasy(f,a,b,n)
    error = abs(Iold-I)
    print("    I = %14.7e (error estimation = %14.7e n = %d)" % (I,error,n))
  return I,n,error

#
#  -3- Mais, il y avait une astuce.....
#      Il faut éviter de calculer plusieurs fois f(x) pour des abscisses identiques !
#      Obtenir I(2n) peut être fait efficacement en ne recalculant que les nouvelles
#      ordonnées et en récupérant les valeurs déjà calculées !
#
  
def trapezeFun(f,a,b,n,nmax,tol):
  error = tol+1
  X,h = linspace(a,b,n+1,retstep=True)
  I = trapezeEasy(f,a,b,n)
  while (error > tol and n*2 < nmax):
    Iold = I; hold = h
    n = n*2; h = h/2
    X = arange(a+h,b,hold)
    I = (I/hold + sum(f(X)))*h
    error = abs(Iold-I)
    print("    I = %14.7e (error estimation = %14.7e n = %d)" % (I,error,n))
  return I,n,error
  
#
# -4- Ensuite, on peut penser à Romberg, mais cela n'était PAS dans les consignes du devoir !
#     Ce sera votre mission 4 : si vous l'avez déjà fait : tant mieux pour vous pour le devoir 4 !
#     Mais cela ne sera pas une bonne idée pour l'évaluation du devoir 3...
#
#