Python : interpolation trigonométrique d'une série de points,

Un billet blog de Denis Hulo

I. Introduction

Étant donné un ensemble de n points représentant des données périodiques de période T, on recherche une fonction d'interpolation qui passe par tous ces points.

639389

On souhaite en fait implémenter une fonction Python qui réalise une interpolation trigonométrique d'une série de points supposés uniformément répartis dans un intervalle de longueur T.


Note : l'objectif n'est pas de donner des démonstrations rigoureuses des différentes formules employées, mais surtout de montrer comment les implémenter en Python.



II. Interpolation trigonométrique

On va tout d'abord montrer comment obtenir un polynôme trigonométrique permettant d'approcher une fonction de période 2𝜋. Puis, on va généraliser ce résultat pour n'importe quelle période.


II-A. Polynôme trigonométrique

On peut exprimer tout polynôme trigonométrique comme une somme de sinus et de cosinus :

639593

Avec une série de n points, on s'attend à avoir une fonction polynomiale avec n paramètres ou coefficients à identifier.

C'est à dire comportant 1 paramètre initial a0, ensuite, le nombre de coefficients supplémentaires dépend de la parité de n :

• Si n est pair, on aura m coefficients ak et m-1 coefficients bk, avec m = n/2.
  
• Si n est impair, on aura m coefficients ak et bk, avec m = (n-1)/2.
  

II-A-1. Analogie entre polynôme et série trigonométrique

Soit f une fonction périodique, de période T=2𝜋, et intégrable sur ℝ, on appelle série de Fourier de f la série trigonométrique telle que :

639594


On voit tout de suite l'analogie avec un polynôme trigonométrique :

639593

On note cependant que ce type de polynôme possède un nombre fini de termes alors qu'une série de Fourier en possède un nombre infini.


II-A-2. Formules des coefficients du polynôme trigonométrique

Les coefficients d'une série de Fourier pour une fonction de période T=2𝜋 sont donnés par :

639711

Si maintenant on se place sur l'intervalle [0, 2𝜋] (avec 𝛼=0), et si on effectue le changement de variable x = i.2𝜋/n dans les formules précédentes, on peut alors passer d'une variable continue x à une variable discrète i.

Les intégrales définies sur [0, 2𝜋] se transforment en sommes de n termes multipliées par un facteur 2𝜋/n.

La fonction continue f peut alors être approchée par une fonction polynomiale trigonométrique à n paramètres.

Pour k=0, on obtient ainsi :

639540


Et pour 0 < k ≤ (n-1)/2 :

639541


Pour 0 ≤ k ≤ m :

639542


Pour k = n/2, et n donc pair, on peut remarquer que :

639544


Comme les termes en sinus de la fonction polynomiale s'annulent tous, on obtient également :

639667


On en déduit enfin que pour k = n/2 :

639629

Cette formule est également valable pour k=0.


II-B. Interpolation sur une période quelconque

Si f(x) est périodique de période 2𝜋, alors f(ωx) est périodique de période 2𝜋/ω.

Par conséquent, si la fonction polynomiale f(x)=P(x) permet d'approcher une fonction continue de période 2𝜋, alors f(ωx) permet d'approcher une fonction continue de période 2𝜋/ω.

La fonction d'interpolation recherchée f(ωx) de période T=2𝜋/ω que l'on note g(x) est donc définie par :

639595



III. Implémentation en Python

Nous souhaitons maintenant traduire en Python les formules précédentes pour nous permettre de déterminer la fonction d'interpolation passant par une série de points régulièrement répartis.


III-A. Évaluation des coefficients ak

On se base sur les formules des coefficients ak vues précédemment.

Pour 0 < k ≤ (n-1)/2 :

639541


Pour k=0 ou k=n/2 :

639629


On obtient facilement la fonction Python :

def coefficient_cosinus(y, k):
# calcul du coefficient ak
# y : liste des valeurs de y
# k : indice du coefficient ak

# nombre de termes du polynôme trigonométrique ou nombre de valeurs de y
n = len(y)

# initialisation de la variable ak
ak =0

# parcours des indices i : 0 -> n-1
for i in range(len(y)):
xi=i*2*math.pi/n # xi = i2π/n
ak += y*math.cos(k*xi) # ak = ak + yi*cos(k*xi)

if k>0 and k<=(n-1)/2:
ak=ak*2/n
else:
ak=ak/n

# retourne la valeur du coefficient ak
return round(ak,12)



III-B. Évaluation des Coefficients bk

De la même façon, en se basant sur la formule des coefficients bk :

639542


On obtient le code Python :

def coefficient_sinus(y, k):
# calcul du coefficient bk
# y : liste des valeurs de y
# k : indice du coefficient bk

# nombre de termes du polynôme trigonométrique ou nombre de valeurs de y
n = len(y)

# initialisation de la variable bk
bk =0

# parcours des indices i : 0 -> n-1
for i in range(n):
xi = i*2*math.pi/n # xi = i2π/n
bk += y*math.sin(k*xi) # bk = bk + yi*sin(k*xi)

bk=bk*2/n

# retourne la valeur du coefficient bk
return round(bk,12)



III-C. Évaluation des coefficients du polynôme

En utilisant les fonctions précédentes, on peut alors écrire le code Python permettant d'obtenir la liste des coefficients de la fonction polynomiale à partir des valeurs de y :

def coefs_polynome_trigo(y):
# évalue les coefficients du polynôme trigo.

n=len(y) # nombre de valeurs de y : nombre de paramètres du polynôme
m = n//2

# intialisation de la liste des paires de coefficients (ak, bk)
coefs=[]

# parcours des indices k : 0 -> m
for k in range(m+1):
# évaluation des coefficients ak et bk
ak = coefficient_cosinus(y,k)
bk = coefficient_sinus(y,k)

# ajout de la paire à la liste
coefs.append((ak,bk))

# retourne la liste des couples de coefficients
return coefs



III-D. Evaluation de la fonction d'interpolation

En se basant finalement sur l'expression générale de la fonction d'interpolation :

639595


On aboutit au code Python :

def eval_fonction_interpolation(coefs,T,x):

ω = 2*math.pi/T # w = 2π/T
valeur_polynome = 0 # initialisation de la variable à retourner

# parcours des indices k : 0 -> m
for k in range(len(coefs)):
# valeur_polynome = valeur_polynome + ak.cos(kωx) + bk.sin(kωx)
valeur_polynome += coefs*math.cos(k*ω*x) + coefs*math.sin(k*ω*x)

# retourne la valeur de la fonction d'interpolation pour x
return valeur_polynome



III-E Module de test

Nous donnons pour finir le module de test complet permettant d'obtenir à partir d'une série de points la fonction d'interpolation :

import math # module pour disposer des fonctions mathématiques
import numpy as np
import matplotlib.pyplot as plt # librairie pour générer le graphique

def f(x):
# fonction supposée inconnue que l'on souhaite approcher en réalisant une interpolation trigonométrique de quelques points
return 5 + 7*math.cos(2*x) + 2*math.sin(2*x) + 3*math.cos(4*x) + 4*math.sin(4*x) + 5*math.cos(6*x) + 12*math.sin(6*x) + 27*math.cos(8*x) + 35*math.sin(8*x)


def coefficient_cosinus(y, k):
# calcul du coefficient ak
# y : liste des valeurs de y
# k : indice du coefficient ak

# nombre de termes du polynôme trigonométrique ou nombre de valeurs de y
n = len(y)

# initialisation de la variable ak
ak =0

# parcours des indices i : 0 -> n-1
for i in range(len(y)):
xi=i*2*math.pi/n # xi = i2π/n
ak += y*math.cos(k*xi) # ak = ak + yi*cos(k*xi)

if k>0 and k<=(n-1)/2:
ak=ak*2/n
else:
ak=ak/n

# retourne la valeur du coefficient ak
return round(ak,12)


def coefficient_sinus(y, k):
# calcul du coefficient bk
# y : liste des valeurs de y
# k : indice du coefficient bk

# nombre de termes du polynôme trigonométrique ou nombre de valeurs de y
n = len(y)

# initialisation de la variable bk
bk =0

# parcours des indices i : 0 -> n-1
for i in range(n):
xi = i*2*math.pi/n # xi = i2π/n
bk += y*math.sin(k*xi) # bk = bk + yi*sin(k*xi)

bk=bk*2/n

# retourne la valeur du coefficient bk
return round(bk,12)


def coefs_polynome_trigo(y):
# évalue les coefficients du polynôme trigo.

n=len(y) # nombre de valeurs de y : nombre de paramètres du polynôme
m = n//2

# intialisation de la liste des paires de coefficients (ak, bk)
coefs=[]

# parcours des indices k : 0 -> m
for k in range(m+1):
# évaluation des coefficients ak et bk
ak = coefficient_cosinus(y,k)
bk = coefficient_sinus(y,k)

# ajout de la paire à la liste
coefs.append((ak,bk))

# retourne la liste des couples de coefficients
return coefs


def eval_fonction_interpolation(coefs,T,x):

ω = 2*math.pi/T # w = 2π/T
valeur_polynome = 0 # initialisation de la variable à retourner

# parcours des indices k : 0 -> m
for k in range(len(coefs)):
# valeur_polynome = valeur_polynome + ak.cos(kωx) + bk.sin(kωx)
valeur_polynome += coefs*math.cos(k*ω*x) + coefs*math.sin(k*ω*x)

# retourne la valeur de la fonction d'interpolation pour x
return valeur_polynome


def expr_fonction_interpolation(coefs,T):

ω = 2*math.pi/T # w = 2π/T
expr_polynome = "f(x) = " + str(coefs) + " + " # initialise de la variable expression avec a0

# parcours des indices k : 0 -> len(coefs)-1
for k in range(1, len(coefs)):
expr_polynome += "{0}.cos({1}x) + {2}.sin({1}x) + ".format(coefs,k*ω,coefs)

# retourne l'expression de la fonction d'interpolation
return expr_polynome[:-3]


# initialisation de la liste y
y=[]

# nombre de termes du polynôme trigonométrique recherché (nombre de points de la série)
n = 9

# période de la fonction d'interpolation recherchée. Par exemple, la fonction f(x) = a0 + a1.cos(2x) + b1.sin(2x) + a2.cos(4x) + b2.sin(4x) + ... a une période de T=π
T = math.pi # T=π

# parcours des indices des valeurs de x et de y
for i in range(n):
xi = (i*T/n) # valeur de xi
yi = f(xi) # valeur de yi = f(xi), f est supposée inconnue

# ajout de la valeur yi à la liste y
y.append(yi)

# évaluation des coefficients du polynôme trigonométrique
coefs_polynome=coefs_polynome_trigo(y)

print("I. Expression de la fonction d'interpolation passant par les différents points :")
print()

# affiche l'exrpession de la fonction d'interpolation passant par les différents points
print(expr_fonction_interpolation(coefs_polynome,T))

print();print()
print("II. Courbe représentative de la fonction d'interpolation passant par les différents points..")

# précision de la trace de la fonction d'interpolation
precision_trace = 10 # 10*n : ici on multiplie le nombre de points initial n par 10

# génère 10*n valeurs de x entre 0 et T
x = list(np.linspace(0.0, T, precision_trace*n))

# initialise la liste y1 avec 10*n cases à None : [None, None, None, ..., None]
y1=*(precision_trace*n)
y2=[]

# parcours des indices de x
for i in range(len(x)):
if i % precision_trace==0: # si i modulo 10 = 0 : 0, 10, 20, ...
y1=y[i//precision_trace] # copie de la valeur mesurée yi
#y1= f(x) # copie de la valeur réelle de la fonction

# valeur donnée par la fonction d'interpo.
yi=eval_fonction_interpolation(coefs_polynome,T,x)
# ajout de yi à la liste y2
y2.append(yi)

# trace les points de la série en rouge
plt.plot(x, y1, marker="o", color="r")
# trace la courbe représentative de la fonction d'interpolation
plt.plot(x, y2, color="b")

# définit le titre du graphique : Interpolation trigonométrique
plt.title("Interpolation trigonométrique")

plt.show() # affiche la figure à l'écran


Le code affiche :

I. Expression de la fonction d'interpolation passant par les différents points :

f(x) = 5.0 + 7.0.cos(2.0x) + 2.0.sin(2.0x) + 3.0.cos(4.0x) + 4.0.sin(4.0x) + 5.0.cos(6.0x) + 12.0.sin(6.0x) + 27.0.cos(8.0x) + 35.0.sin(8.0x)


Puis la courbe représentative de la fonction passant par les différents points :

639467


IV. Conclusion

Après avoir montré l'analogie entre les séries de Fourier et les polynômes trigonométriques, nous avons pu obtenir les formules des coefficients du polynôme, puis celle de la fonction d'interpolation.

Enfin, nous avons su traduire ces formules en Python, et les utiliser pour déterminer à partir d'une série de points, la fonction d'interpolation.


Sources :

https://fr.wikipedia.org/wiki/Polyn%C3%B4me_trigonom%C3%A9trique
https://fr.wikipedia.org/wiki/S%C3%A9rie_de_Fourier