lagrange.py
#
# interpolation de Lagrange : naif ou efficace ?
# tic, toc : implémentation Python des fonctions tic, toc de Matlab
#
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#
from numpy import *
import time
import matplotlib
from matplotlib import pyplot as plt
# =========================================================================
def tic():
global startTime
startTime = time.time()
# =========================================================================
def toc(message = ''):
global startTime
stopTime = time.time()
if message:
message = ' (' + message + ')' ;
print("Elapsed time is %.6f seconds %s" % ((stopTime - startTime),message) )
startTime = 0
# =========================================================================
def lagrange_naive(x,X,U):
n = min(len(X),len(U))
m = len(x)
uh = zeros(m)
phi = zeros(n)
for j in range(m):
uh[j] = 0
for i in range(n):
phi[i] = 1.0
for k in range(n):
if i != k:
phi[i] = phi[i]*(x[j]-X[k])/(X[i]-X[k])
uh[j] = uh[j] + U[i] * phi[i]
return uh
# =========================================================================
def lagrange(x,X,U):
n = min(len(X),len(U)) # Evite un message d'erreur :-)
phi = ones((n,len(x))) # Preallocation (important !)
for i in range(n):
for k in list(range(0,i)) + list(range(i+1,n)):
phi[i,:] = phi[i,:]*(x-X[k])/(X[i]-X[k])
return U @ phi # Produit matriciel (plus rapide !)
# ============================= mainProgram ===============================
matplotlib.rcParams['toolbar'] = 'None'
plt.rcParams['figure.facecolor'] = 'silver'
plt.figure('Polynomial interpolation')
#
# -1- x : c'est quoi ?
#
X = [0,1,2]
U = [0,0,1]
x = linspace(0,2,9)
uh = lagrange_naive(x,X,U);
plt.plot(x,zeros(size(x)),'.r',markersize=15)
plt.plot(x,uh,'.-b',markersize=15)
plt.plot(X,U,'.k',markersize=25)
plt.show()
#
# -2- Interpoler une fonction
#
# n = 15
# h = 1/n
# x = linspace(-1,1,2000)
# Xunif = arange(-1+h/2,1+h/2,h)
# Xcheb = cos( pi*(2*arange(0,2*n)+1) / (4*n) )
#
# u = lambda x : 1/(x**2 + 0.09)
# Uunif = u(Xunif)
# Ucheb = u(Xcheb)
# plt.figure('Polynomial interpolation')
# plt.plot(x,lagrange(x,Xunif,Uunif),'-b',label='Abscisses équidistantes')
# plt.plot(x,lagrange(x,Xcheb,Ucheb),'-r',label='Abscisses de Tchebychev')
# plt.plot(Xunif,Uunif,'ob',Xcheb,Ucheb,'or')
# plt.xlim((-1,1))
# plt.ylim((0,max(Uunif)*2))
# plt.title('Abscisses : n = %d (uniforme) %d (Tchebychev)'
# % (len(Xunif),len(Xcheb)))
# plt.legend(loc='upper right')
# plt.show()
#
# -3- Interpoler plusieurs fonctions
#
n = 15
h = 1/n
x = linspace(-1,1,2000)
Xunif = arange(-1+h/2,1+h/2,h)
Xcheb = cos( pi*(2*arange(0,2*n)+1) / (4*n) )
functions = [lambda x : 1/(x**2 + 0.09),
lambda x : abs(x - 0.1)]
for u in functions:
Uunif = u(Xunif)
Ucheb = u(Xcheb)
plt.figure('Polynomial interpolation')
plt.plot(x,lagrange(x,Xunif,Uunif),'-b',label='Abscisses équidistantes')
plt.plot(x,lagrange(x,Xcheb,Ucheb),'-r',label='Abscisses de Tchebychev')
plt.plot(Xunif,Uunif,'ob',Xcheb,Ucheb,'or')
plt.xlim((-1,1))
plt.ylim((0,max(Uunif)*2))
plt.title('Abscisses : n = %d (uniforme) %d (Tchebychev)'
% (len(Xunif),len(Xcheb)))
plt.legend(loc='upper right')
plt.show()
#
# -4- Estimation du temps CPU
#
# Elapsed time is 0.000219 seconds (Lagrange n = 2000)
# Elapsed time is 0.017949 seconds (LagrangeNaive n = 2000)
# Elapsed time is 0.088590 seconds (Lagrange n = 2000000)
# Elapsed time is 16.298496 seconds (LagrangeNaive n = 2000000)
# Elapsed time is 1.415670 seconds (Lagrange n = 20000000)
# Elapsed time is 159.760674 seconds (LagrangeNaive n = 20000000)
#
for n in [2000,20000]:
x = linspace(-1,1,n)
tic()
lagrange(x,[0,1,2],[0,1,2])
toc("Lagrange n = %d" % n)
tic()
lagrange_naive(x,[0,1,2],[0,1,2])
toc("LagrangeNaive n = %d" % n)
# =========================================================================