sixPointsChebyshev.py
#
# Six points aux abscisses de Chebyshev
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#
from numpy import *
import matplotlib
from matplotlib import pyplot as plt
#from matplotlib.pyplot import *
matplotlib.rcParams['toolbar'] = 'None'
plt.rcParams['figure.facecolor'] = 'white'
#
# -1- Six points aux abscisses de Chebyshev
#
n = 5
X = cos(pi * (2*arange(0,n+1) + 1)/((n+1) * 2))
U = [0.2,0.4,0.6,0.7,0.5,0.8]
puh = polyfit(X,U,5)
x = linspace(-1,1,200)
uh = polyval(puh,x)
plt.plot(x,uh,'-b')
plt.plot(X,U,'or')
plt.plot([x[0],x[-1]],[0,0],'-k')
plt.plot(X,zeros(size(X)),'or')
plt.axis('equal')
plt.axis('off')
plt.figure()
#
# -2- Ajoutons un septième point
#
pu = polyfit([*X,0],[*U,0.3],6)
u = polyval(pu,x)
plt.plot(x,u,'-g')
plt.plot([0],[0.3],'og')
plt.plot(x,uh,'-b')
plt.plot(X,U,'or')
plt.plot([x[0],x[-1]],[0,0],'-k')
plt.axis('equal')
plt.axis('off')
plt.figure()
#
# -3- Dessinons l'erreur
#
error = u - uh
E = error[0]
plt.plot(x,u,'-g')
plt.plot([0.0],[0.3],'og')
plt.plot(x,uh,'-b')
plt.plot(X,U,'or')
plt.plot([x[0],x[-1]],[0,0],'-k')
plt.plot(x,error,'-r');
plt.axis('equal')
plt.axis('off')
plt.figure()
plt.plot(x,error,'-r')
plt.plot([x[0],x[-1]],[0,0],'-k')
plt.plot(X,zeros(size(X)),'or')
plt.plot([x[0]-0.3,x[0]+0.3],[E,E],'-r')
plt.axis('equal')
plt.axis('off')
plt.figure()
#
# -4- Derivons et divisons l'erreur
#
peh = pu - [0,*puh]
eh = polyval(peh, x)
pdeh = peh * [6,5,4,3,2,1,0]
#pdeh = pdeh[0:-1]
pdeh = pdeh[0:6]
print(pdeh)
deh = polyval(pdeh,x)/6
Xd = roots(pdeh)
Ed = polyval(peh,Xd)
E = Ed[0]
plt.plot([x[0],x[-1]],[0,0],'-k')
plt.plot(x,eh,'-r')
plt.plot(x,deh,'-b')
plt.plot(X,zeros(size(X)),'or',markersize=3.5)
plt.plot(Xd,zeros(size(Xd)),'ob',markersize=3.5)
plt.plot(Xd,Ed,'or',markersize=3.5)
plt.axis('off')
plt.axis('equal')
plt.figure()
#
# -5- Et encore quelques petites manipulations
#
# On va superposer les deux expressions de l'erreur
# et la courbe rouge est recouverte par la courbe bleue :-)
#
plt.plot([x[0],x[-1]],[0,0],'-k')
plt.plot(x,E**2-eh**2,'-r')
plt.plot(x,(1-x**2)*(deh**2),'-b')
plt.plot(Xd,zeros(size(Xd)),'ok')
plt.plot([-1,1],[0,0],'or')
plt.axis('off')
plt.axis('equal')
plt.show()