T.D. Informatique / simulation numérique

TD_simulnum.md

Code en Python3 pour les T.D. d'informatique de mars 2014

On trouvera :

• TD_simulnum_gauss.py : méthode du pivot de Gauss
• TD_simulnum_riemann.py : méthodes d'intégration
• TD_simulnum_solve.py : recherche du zéro d'une fonction

TD_simulnum_gauss.py

# -*- coding: utf8 -*-

# Pivot de Gauss
# ==============

# la fonction prend en argument m une matrice sous la forme d'une liste
# de listes (chaque liste constituant une ligne de la matrice) :
#
# Par exemple la matrice identité 3x3 sera : [[1,0,0],[0,1,0],[0,0,1]]
#
def pivot(m):
  nbr_rows = len(m) # nombre de lignes
  for i in range(nbr_rows):
    # on cherche le meilleur pivot, c'est-à-dire le plus grand terme
    # (en valeur absolue) pour la colonne i parmi les lignes non encore
    # utilisées ; compte-tenu de l'algorithme, les colonnes de gauche
    # (sauf pour la première itération évidemment) comporteront
    # nécessairement une valeur nulle.
    pivot = i # par défaut le pivot se trouve sur la ligne i
    for j in range(i, nbr_rows):
      if abs(m[j][i]) > abs(m[pivot][i]):
        pivot = j
    # s'il n'y a aucune colonne non nulle, la matrice est singulière
    if m[pivot][i] == 0:
      raise Exception("Singular matrix")
    # on échange les lignes i et pivot
    tmp = m[i]      # variable temporairement utile pour l'échange
    m[i] = m[pivot]
    m[pivot] = tmp
    # on normalise la ligne i (en mettant à 1 la colonne i)
    p = m[i][i]
    m[i] = [ x/p for x in m[i] ]
    # on annule la colonne i dans toutes les autres lignes
    for j in range(0,nbr_rows):
      if j != i:
        m[j] = [ x - m[j][i]*m[i][k] for k,x in enumerate(m[j]) ]
  return m


# Exemples
# ========
#
# 1) Inversion de la matrice
#     
#      /  2 -1  0 \                            /  2 -1  0 | 1  0  0 \
#      | -1  2 -1 |  en appliquant le pivot à  | -1  2 -1 | 0  1  0 |
#      \  0 -1  2 /                            \  0 -1  2 | 0  0  1 /

print(pivot([[2,-1,0,1,0,0],[-1,2,-1,0,1,0],[0,-1,2,0,0,1]]))


# 2) Résolution d'un système d'équations linéaires :
#            ,
#           /    x -   y + 2.z  =   5
#          +|  3.x + 2.y +   z  =  10
#           \  2.x - 3.y - 2.z  = -10
#            `
#    par application de l'algorithme du pivot à la matrice 
#
#           / 1 -1  2 |   5 \
#           | 3  2  1 |  10 |
#           \ 2 -3 -2 | -10 /

print(pivot([[1,-1,2,5],[3,2,1,10],[2,-3,-2,-10]]))

TD_simulnum_riemann.py

# -*- coding: utf8 -*-

# Intégration d'une fonction entre deux bornes selon plusieurs
# types d'approximation :
#
#   * integr_rect1( f, a, b, epsilon )
#     Intégration par des rectangles dont la hauteur coïncide avec
#     la courbe sur sa gauche
#
#   * integr_rect2( f, a, b, epsilon )
#     Intégration par des rectangles dont la hauteur est la hauteur moyenne
#     des deux points de la courbe à gauche et à droite
#     (cela revient en fait à construire un trapèze)
#
#   * integr_rect3( f, a, b, epsilon )
#     Intégration par des rectangles dont la hauteur est la plus petite
#     des hauteurs des deux points de la courbe à gauche et à droite
#
#   * integr_rect4( f, a, b, epsilon )
#     Intégration par des rectangles dont la hauteur est la plus grande
#     des hauteurs des deux points de la courbe à gauche et à droite
#
#   * integr_trap( f, a, b, epsilon )
#     Intégration par des trapèzes (en fait identique à integr_rect2)

def integr_rect1(f, a, b, epsilon):
  s = 0                  # somme des aires
  n = int((b-a)/epsilon) # nombre d'intervalles de largeur epsilon
  for i in range(n):
    x = a + i*epsilon
    s += f(x) * epsilon  # hauteur x largeur
  return s

def integr_rect2(f, a, b, epsilon):
  s = 0                  # somme des aires
  n = int((b-a)/epsilon) # nombre d'intervalles de largeur epsilon
  for i in range(n):
    x = a + i*epsilon
    s += ( f(x)+f(x+epsilon) ) / 2 * epsilon  # hauteur x largeur
  return s

def integr_rect3(f, a, b, epsilon):
  s = 0                  # somme des aires
  n = int((b-a)/epsilon) # nombre d'intervalles de largeur epsilon
  for i in range(n):
    x = a + i*epsilon
    s += min( f(x),f(x+epsilon) ) * epsilon  # hauteur x largeur
  return s

def integr_rect4(f, a, b, epsilon):
  s = 0                  # somme des aires
  n = int((b-a)/epsilon) # nombre d'intervalles de largeur epsilon
  for i in range(n):
    x = a + i*epsilon
    s += max( f(x),f(x+epsilon) ) * epsilon  # hauteur x largeur
  return s

integr_trap = integr_rect2

# exemples
# ========

from math import sin, atan, pi
f = lambda x: x*atan(x)
print( "integr_rect1    =",integr_rect1(f,0,1,1e-5) )
print( "integr_rect2    =",integr_rect2(f,0,1,1e-5) )
print( "integr_rect3    =",integr_rect3(f,0,1,1e-5) )
print( "integr_rect4    =",integr_rect4(f,0,1,1e-5) )
print( "integr_trap     =",integr_trap(f,0,1,1e-5) )
print( "valeur attendue =",pi/4-.5)

f = lambda x: sin(x)**2
print( "integr_rect1    =",integr_rect1(f,0,pi/2,1e-5) )
print( "integr_rect2    =",integr_rect2(f,0,pi/2,1e-5) )
print( "integr_rect3    =",integr_rect3(f,0,pi/2,1e-5) )
print( "integr_rect4    =",integr_rect4(f,0,pi/2,1e-5) )
print( "integr_trap     =",integr_trap(f,0,pi/2,1e-5) )
print( "valeur attendue =",pi/4)

TD_simulnum_solve.py

# -*- coding: utf8 -*-


# résolution par dichotomie
# =========================
#
#   arguments : fonction, borne gauche, borne droite, epsilon
#
def solve_dichot(f,a,b,epsilon):
  a,b = min(a,b), max(a,b) # on remet a et b dans le bon ordre
  signe_a = f(a) < 0         # True si a est négatif
  signe_b = f(b) < 0         # True si b est négatif
  if not signe_a ^ signe_b:
    raise Exception("les deux bornes n'encadrent pas un zéro de la fonction")
  while b-a > epsilon:
    c = (a+b)/2
    signe_c = f(c) < 0
    if signe_a ^ signe_c:
      b = c # f(a) et f(c) sont de signes opposés
    else:
      a = c # f(a) et f(c) sont de signes identiques
  return (a+b)/2




# Newton-Raphson (avec la dérivée en argument)
# ============================================
#
#   arguments : fonction, dérivée, x initial, epsilon
#
def solve_newton1(f,ff,x,epsilon):
  d = 8*epsilon # initialisation à une valeur "élevée"
  while abs(d) > epsilon:
    d = f(x)/ff(x)
    x = x - d
  return x




# Newton-Raphson (sans la dérivée en argument)
# ============================================
#
#   arguments : fonction, x initial, epsilon
#
def solve_newton2(f,x,epsilon):
  d = 8*epsilon # initialisation à une valeur "élevée"
  while abs(d) > epsilon:
    deriv = ( f(x+epsilon/2.0) - f(x-epsilon/2.0) ) / epsilon
    d = f(x)/deriv
    x = x - d
  return x




# test des fonctions
# ------------------
from math import log
print( solve_dichot(lambda x: log(x) - 1, 2.0, 4.0, 1e-8) )
print( solve_newton1(lambda x: log(x) - 1, lambda x: 1.0/x, 4.0, 1e-8) )
print( solve_newton2(lambda x: log(x) - 1, 4.0, 1e-8) )

from math import cos, sin
print( solve_dichot(lambda x: cos(x), 0.5, 3.0, 1e-8) )
print( solve_newton1(lambda x: cos(x), lambda x: -sin(x), 0.5, 1e-8) )
print( solve_newton2(lambda x:cos(x), 0.5, 1e-8) )

import cProfile
cProfile.run("solve_dichot(lambda x: log(x) - 1, 2.0, 4.0, 1e-8)")
cProfile.run("solve_newton1(lambda x: log(x) - 1, lambda x: 1.0/x, 4.0, 1e-8)")
cProfile.run("solve_newton2(lambda x: log(x) - 1, 4.0, 1e-8)")