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)
   8     diff_div=np.zeros((n,n))
   9     for i in range(n):
  10         diff_div[i][0] = Y[i]
  11     for j in range(1,n):
  12         for i in range(n-j):
  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     print "les coeffs sont",d
  38     n = len(X)-1
  39     return lambda x : \
  40         sum([d[i]*prod([x - X[j] for j in xrange(i)]) for i in xrange(n+1)])
  41
  42
  43
  44
  45
  46
  47 print "tp 2.1....."
  48 X = [10,25,60]
  49 XX= [15,40,100]
  50 Y = [2.3, 8, 24.6]
  51 p = construit_pol(X,Y)
  52
  53 Yp =[p(x) for x in XX]
  54 Ypp = [construit_et_eval_pol(X,Y,x) for x in XX]
  55
  56 print "Yp", Yp#, Ypp
  57
  58
  59 x=np.linspace(10,100,100)
  60 plt.plot(x,p(x))
  61 plt.ylabel('temps d\'execution')
  62 plt.xlabel("taille des donnees")
  63 plt.show()
  64
  65
  66
  67
  68
  69 print p(15)
  70 print p(40)
  71 print p(100)
  72
  73
  74 """
  75
  76
  77
  78 print "tp 2.2....."
  79 X1 = [-1+0.25*i for i in xrange(9)]
  80 Y1 = [e**x for x in X1]
  81 p1 = construit_pol(X1,Y1)
  82
  83
  84 X2 = [cos(pi/2*(float(2*i+1)/9)) for i in xrange(9)]
  85 Y2 = [e**x for x in X2]
  86 p2 = construit_pol(X2,Y2)
  87
  88
  89
  90 e1 = lambda x : abs(p1(x) - e**x)
  91 e2 = lambda x : abs(p2(x) - e**x)
  92
  93
  94 l1 = [e1(-1 + 0.001*i) for i in range(2001)]
  95 l2 = [e2(-1 + 0.001*i) for i in range(2001)]
  96
  97 average = lambda l : sum(l)/len(l)
  98
  99 print "l1", max(l1), average(l1)
 100 print "l2", max(l2), average(l2)
 101
 102
 103 x=np.linspace(-1,1,1000)
 104 plt.plot(x,e1(x))
 105 plt.plot(x,e2(x))
 106
 107 plt.ylabel('exp')
 108 plt.xlabel("x")
 109 plt.show()
 110
 111 """