src/MEDOP/tut/medcoupling/pyfunctions/lagrange.py

   1 #!/usr/bin/env python
   2 import scipy
   3 def lagrange(points):
   4     '''
   5     This returns a polynom that fits the points specified in the input
   6     dictionary (lagrange interpolation). In this dictionary, the keys
   7     are x value and the values are y corresponding values
   8     (i.e. y=polynom(x)). The polynom is a scipy polynom and then is a
   9     callable (can be used as a function).
  10     '''
  11     tmp = scipy.poly1d([0])
  12     result=scipy.poly1d([0])
  13
  14     for i in points.keys():
  15         numerator=scipy.poly1d([1])
  16         denom = 1.0
  17         for j in points.keys():
  18             if (i != j):
  19                 tmp = scipy.poly1d([1,-j])
  20                 numerator = numerator * tmp
  21                 denom = denom * (i - j)
  22         tmp = (numerator/denom) * points.get(i)
  23         result = result + tmp
  24
  25     return result
  26
  27 def points_usingfunction(arrX,function):
  28     points={}
  29     for x in arrX:
  30         points[x] = function(x)
  31     return points
  32
  33 def points_usingarray(arrX,arrY):
  34     points={}
  35     for i in range(len(arrX)):
  36         x=arrX[i]
  37         y=arrY[i]
  38         points[x] = y
  39     return points
  40
  41 def sortdict(points):
  42     # Sort this dictionary by keys and returns 2 lists, the list of X
  43     # and the list of Y, the whole ordered by X
  44     keys = points.keys()
  45     keys.sort()
  46     return keys, [points[key] for key in keys]
  47
  48 import pylab
  49 import numpy
  50 def plot(function, start=0., stop=1., step=0.01):
  51     """
  52     The parameter function must be a callable.
  53     """
  54     arrX=numpy.arange(start, stop, step, dtype='float64')
  55     # function is a callable
  56     arrY=map(function,arrX)
  57     pylab.plot(arrX, arrY)
  58     pylab.show()
  59
  60
  61 # ---
  62 # The points does not need to be ordered by x values
  63 def TEST_lagrange_01():
  64     input = {0.:5, 0.2:10, 0.9:10, 0.6:21, 1:8}
  65     polynom = lagrange(input)
  66     print polynom
  67     plot(function=polynom, start=0., stop=1., step=0.001)
  68
  69 def TEST_lagrange_02():
  70     input = {0.:0., 0.5:1., 1.:0.}
  71     polynom = lagrange(input)
  72     print polynom
  73     plot(function=polynom, start=0., stop=1., step=0.001)
  74
  75 # ---
  76 # One can create the input dictionary  from arrays
  77 def TEST_lagrange_usingarrays_01():
  78     arrX = [0., 0.2, 0.9, 0.6, 1]
  79     arrY = [5, 10, 10, 21, 8]
  80     input = points_usingarray(arrX,arrY)
  81     polynom = lagrange(input)
  82     print polynom
  83     plot(function=polynom, start=0., stop=1., step=0.001)
  84
  85 # Another example using numpy
  86 def TEST_lagrange_usingarrays_02():
  87     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
  88     arrY=numpy.zeros(len(arrX), dtype='float64')
  89     arrY[3]=2
  90     input = points_usingarray(arrX,arrY)
  91     polynom = lagrange(input)
  92     print polynom
  93     plot(function=polynom, start=0., stop=1., step=0.001)
  94
  95 # ---
  96 # One can create the input dictionary  from a function applied to an
  97 # array of X values
  98
  99 # simple method for mathematical functions
 100 def TEST_lagrange_usingfunction_01():
 101     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
 102     arrY=numpy.cos(10*arrX)
 103     input = points_usingarray(arrX,arrY)
 104     polynom = lagrange(input)
 105     print polynom
 106     plot(function=polynom, start=0., stop=1., step=0.001)
 107
 108 # General method
 109 xlimit=0.8
 110 def chapeau(x):
 111     if x<xlimit:
 112         y=x
 113     else:
 114         y=2*xlimit-x
 115     return y
 116
 117 def TEST_lagrange_usingfunction_01():
 118     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
 119     input = points_usingfunction(arrX,chapeau)
 120     polynom = lagrange(input)
 121     print polynom
 122     plot(function=polynom, start=0., stop=1., step=0.001)
 123
 124
 125 if __name__ == "__main__":
 126     #TEST_lagrange_01()
 127     TEST_lagrange_02()
 128     #TEST_lagrange_usingarrays_01()
 129     #TEST_lagrange_usingarrays_02()
 130     #TEST_lagrange_usingfunction_01()
 131     #TEST_lagrange_usingfunction_01()