tps/chap2/diff_div.py

   1 import matplotlib.pyplot as plt
   2 import numpy as np
   3 from math import *
   4
   5
   6 def coeffs(X,Y):
   7     n = len(X)-1
   8     diff_div=[[0 for _ in xrange(n+1)] for _ in xrange(n+1)]
   9     for i in xrange(n+1):
  10         diff_div[i][0] = Y[i]
  11     for j in xrange(1,n+1):
  12         for i in xrange(n-j+1):
  13             r = float(diff_div[i+1][j-1]-diff_div[i][j-1])/(X[i+j]-X[i])
  14             diff_div[i][j] = r
  15     print diff_div
  16     return diff_div[0]
  17
  18
  19
  20
  21 def prod(list):
  22     r = 1
  23     for x in list:
  24         r *= x
  25     return r
  26
  27
  28 def construit_et_eval_pol(X,Y,x):
  29     d = coeffs(X,Y)
  30     n = len(X)-1
  31     return sum([d[i]*prod([x - X[j] for j in xrange(i)]) for i in xrange(n+1)])
  32
  33
  34
  35 def construit_pol(X,Y):
  36     d = coeffs(X,Y)
  37     n = len(X)-1
  38     return lambda x : \
  39         sum([d[i]*prod([x - X[j] for j in xrange(i)]) for i in xrange(n+1)])
  40
  41
  42
  43
  44 """
  45
  46 print "tp 2.1....."
  47 X = [10,25,60]
  48 XX= [15,40,100]
  49 Y = [2.3, 8, 24.6]
  50 p = construit_pol(X,Y)
  51 Yp =[p(x) for x in XX]
  52 Ypp = [construit_et_eval_pol(X,Y,x) for x in XX]
  53
  54 print Yp, Ypp
  55
  56
  57 x=np.linspace(10,100,100)
  58 plt.plot(x,p(x))
  59 plt.ylabel('temps d\'execution')
  60 plt.xlabel("taille des donnees")
  61 plt.show()
  62
  63
  64
  65
  66
  67 print p(15)
  68 print p(40)
  69 print p(100)
  70
  71
  72 """
  73
  74
  75
  76 print "tp 2.2....."
  77 X1 = [-1+0.25*i for i in xrange(9)]
  78 Y1 = [e**x for x in X1]
  79 p1 = construit_pol(X1,Y1)
  80
  81
  82 X2 = [cos(pi/2*(float(2*i+1)/9)) for i in xrange(9)]
  83 Y2 = [e**x for x in X2]
  84 p2 = construit_pol(X2,Y2)
  85
  86
  87
  88 e1 = lambda x : abs(p1(x) - e**x)
  89 e2 = lambda x : abs(p2(x) - e**x)
  90
  91
  92 l1 = [e1(-1 + 0.001*i) for i in range(2001)]
  93 l2 = [e2(-1 + 0.001*i) for i in range(2001)]
  94
  95 average = lambda l : sum(l)/len(l)
  96
  97 print "l1", max(l1), average(l1)
  98 print "l2", max(l2), average(l2)
  99
 100
 101 x=np.linspace(-1,1,1000)
 102 plt.plot(x,e1(x))
 103 plt.plot(x,e2(x))
 104
 105 plt.ylabel('exp')
 106 plt.xlabel("x")
 107 plt.show()
 108
 109