doc/tut/medcoupling/pyfunctions/lagrange.py

   1 #!/usr/bin/env python3
   2 # Copyright (C) 2012-2020  CEA/DEN, EDF R&D
   3 #
   4 # This library is free software; you can redistribute it and/or
   5 # modify it under the terms of the GNU Lesser General Public
   6 # License as published by the Free Software Foundation; either
   7 # version 2.1 of the License, or (at your option) any later version.
   8 #
   9 # This library is distributed in the hope that it will be useful,
  10 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 # Lesser General Public License for more details.
  13 #
  14 # You should have received a copy of the GNU Lesser General Public
  15 # License along with this library; if not, write to the Free Software
  16 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 #
  18 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  19 #
  20
  21 import scipy
  22 def lagrange(points):
  23     '''
  24     This returns a polynom that fits the points specified in the input
  25     dictionary (lagrange interpolation). In this dictionary, the keys
  26     are x value and the values are y corresponding values
  27     (i.e. y=polynom(x)). The polynom is a scipy polynom and then is a
  28     callable (can be used as a function).
  29     '''
  30     tmp = scipy.poly1d([0])
  31     result=scipy.poly1d([0])
  32
  33     for i in points:
  34         numerator=scipy.poly1d([1])
  35         denom = 1.0
  36         for j in points:
  37             if (i != j):
  38                 tmp = scipy.poly1d([1,-j])
  39                 numerator = numerator * tmp
  40                 denom = denom * (i - j)
  41         tmp = (numerator/denom) * points.get(i)
  42         result = result + tmp
  43
  44     return result
  45
  46 def points_usingfunction(arrX,function):
  47     points={}
  48     for x in arrX:
  49         points[x] = function(x)
  50     return points
  51
  52 def points_usingarray(arrX,arrY):
  53     points={}
  54     for i in range(len(arrX)):
  55         x=arrX[i]
  56         y=arrY[i]
  57         points[x] = y
  58     return points
  59
  60 def sortdict(points):
  61     # Sort this dictionary by keys and returns 2 lists, the list of X
  62     # and the list of Y, the whole ordered by X
  63     keys = sorted(points.keys())
  64     return keys, [points[key] for key in keys]
  65
  66 import pylab
  67 import numpy
  68 def plot(function, start=0., stop=1., step=0.01):
  69     """
  70     The parameter function must be a callable.
  71     """
  72     arrX=numpy.arange(start, stop, step, dtype='float64')
  73     # function is a callable
  74     arrY=list(map(function,arrX))
  75     pylab.plot(arrX, arrY)
  76     pylab.show()
  77
  78
  79 # ---
  80 # The points does not need to be ordered by x values
  81 def TEST_lagrange_01():
  82     input = {0.:5, 0.2:10, 0.9:10, 0.6:21, 1:8}
  83     polynom = lagrange(input)
  84     print(polynom)
  85     plot(function=polynom, start=0., stop=1., step=0.001)
  86
  87 def TEST_lagrange_02():
  88     input = {0.:0., 0.5:1., 1.:0.}
  89     polynom = lagrange(input)
  90     print(polynom)
  91     plot(function=polynom, start=0., stop=1., step=0.001)
  92
  93 # ---
  94 # One can create the input dictionary  from arrays
  95 def TEST_lagrange_usingarrays_01():
  96     arrX = [0., 0.2, 0.9, 0.6, 1]
  97     arrY = [5, 10, 10, 21, 8]
  98     input = points_usingarray(arrX,arrY)
  99     polynom = lagrange(input)
 100     print(polynom)
 101     plot(function=polynom, start=0., stop=1., step=0.001)
 102
 103 # Another example using numpy
 104 def TEST_lagrange_usingarrays_02():
 105     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
 106     arrY=numpy.zeros(len(arrX), dtype='float64')
 107     arrY[3]=2
 108     input = points_usingarray(arrX,arrY)
 109     polynom = lagrange(input)
 110     print(polynom)
 111     plot(function=polynom, start=0., stop=1., step=0.001)
 112
 113 # ---
 114 # One can create the input dictionary  from a function applied to an
 115 # array of X values
 116
 117 # simple method for mathematical functions
 118 def TEST_lagrange_usingfunction_01():
 119     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
 120     arrY=numpy.cos(10*arrX)
 121     input = points_usingarray(arrX,arrY)
 122     polynom = lagrange(input)
 123     print(polynom)
 124     plot(function=polynom, start=0., stop=1., step=0.001)
 125
 126 # General method
 127 xlimit=0.8
 128 def chapeau(x):
 129     if x<xlimit:
 130         y=x
 131     else:
 132         y=2*xlimit-x
 133     return y
 134
 135 def TEST_lagrange_usingfunction_01():
 136     arrX=numpy.arange(start=0., stop=1., step=0.1, dtype='float64')
 137     input = points_usingfunction(arrX,chapeau)
 138     polynom = lagrange(input)
 139     print(polynom)
 140     plot(function=polynom, start=0., stop=1., step=0.001)
 141
 142
 143 if __name__ == "__main__":
 144     #TEST_lagrange_01()
 145     TEST_lagrange_02()
 146     #TEST_lagrange_usingarrays_01()
 147     #TEST_lagrange_usingarrays_02()
 148     #TEST_lagrange_usingfunction_01()
 149     #TEST_lagrange_usingfunction_01()