§ 3	Ajustements au sens des moindres carrés

§ 3.1	Ajustement par une combinaison linéaire de fonctions à une variable

Exemple:	Ajustement d'une droite

Enoncé du problème d'ajustement au sens des moindres carrés

Une expérience de physique consiste à mesurer le volume y d'un gaz, à pression constante, en fonction de la température x:

<<Tableaux`;
afficheTableau[{"x [°C]","y [(cm^3)]"},None,{{"(Subscript[x, 1])","(Subscript[x, 2])","..","(Subscript[x, n])"},{"(Subscript[y, 1])","(Subscript[y, 2])","..","(Subscript[y, n])"}}]

n désigne le nombre de mesures effectuées.  Figure pour n=5:

x={-3,-1,2,5,8};y={3.2`,4.7`,5,6.3`,6.7`}-3;
data=Transpose[{x,y}];
Show[Graphics[{PointSize[0.015`],Point/@data},Axes->True,AxesOrigin->{0,0},PlotRange->{0,4},ImageSize->{400,300},Ticks->{{{x[[1]],"\!\(\*SubscriptBox[\(x\), \(1\)]\)"},{x[[2]],"\!\(\*SubscriptBox[\(x\), \(2\)]\)"},{x[[3]],"\!\(\*SubscriptBox[\(x\), \(3\)]\)"},{x[[4]],"\!\(\*SubscriptBox[\(x\), \(4\)]\)"},{x[[5]],"\!\(\*SubscriptBox[\(x\), \(5\)]\)"}},{{y[[1]],"\!\(\*SubscriptBox[\(y\), \(1\)]\)"},{y[[2]],"\!\(\*SubscriptBox[\(y\), \(2\)]\)"},{y[[3]],"\!\(\*SubscriptBox[\(y\), \(3\)]\)"},{y[[4]],"\!\(\*SubscriptBox[\(y\), \(4\)]\)"},{y[[5]],"\!\(\*SubscriptBox[\(y\), \(5\)]\)"}}}]]

Pour des raisons physiques, on postule que y est une fonction affine de x:

y=a x+b

Etant donné l'imprécision des mesures, la loi théorique n'est vérifiée qu'approximativement:

Les écarts suivants sont appelés "résidus"

(Subscript[e, 1])=(Subscript[y, 1])-(a (Subscript[x, 1])+b), (Subscript[e, 2])=(Subscript[y, 2])-(a (Subscript[x, 2])+b), ...,(Subscript[e, n])=(Subscript[y, n])-(a (Subscript[x, n])+b)

"Ajuster la droite au sens des moindres carrés" signifie qu'il faut déterminer a, b de manière que la somme des carrés des résidus soit minimale

(Subscript[e, 1]^2)+(Subscript[e, 2]^2)+...+(Subscript[e, n]^2)=min

Réduction du problème d'ajustement à un problème de projection

(Subscript[e, 1]^2)+(Subscript[e, 2]^2)+...+(Subscript[e, n]^2)=((((Subscript[y, 1])-a (Subscript[x, 1])-b))^2)+((((Subscript[y, 2])-a (Subscript[x, 2])-b))^2)+...+((((Subscript[y, n])-a (Subscript[x, n])-b))^2)=(((∥(({{(Subscript[y, 1])-a (Subscript[x, 1])-b}, {(Subscript[y, 2])-a (Subscript[x, 2])-b}, {...}, {(Subscript[y, n])-a (Subscript[x, n])-b}}))∥))^2)=(((∥(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {...}, {(Subscript[y, n])}}))-a(({{(Subscript[x, 1])}, { (Subscript[x, 2])}, {...}, {(Subscript[x, n])}}))-b(({{1}, {1}, {...}, {1}}))∥))^2)

Autrement dit, dans (ℝ^n), on donne le vecteur

(OverVector[y])=(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {...}, {(Subscript[y, n])}}))

et le sous-espace G, engendré par les vecteurs

(OverVector[(Subscript[g, 1])])=(({{(Subscript[x, 1])}, { (Subscript[x, 2])}, {...}, {(Subscript[x, n])}})),      (OverVector[(Subscript[g, 2])])=(({{1}, {1}, {...}, {1}})),dont les éléments sont     (OverVector[g])=a  (OverVector[(Subscript[g, 1])])+b (OverVector[(Subscript[g, 2])])

Déterminez les nombres a, b tels que

(((∥(OverVector[y])-(OverVector[g])∥))^2) =  min

Nous avons étudié, dans le § 2, que la solution est la projection orthogonale (OverVector[g]) de (OverVector[y]) sur G déterminée par les conditions suivantes

(OverVector[g])=a (OverVector[(Subscript[g, 1])])+b (OverVector[(Subscript[g, 2])])
 (OverVector[e])=(OverVector[y])-(OverVector[g])=(OverVector[y])-a (OverVector[(Subscript[g, 1])])-b (OverVector[(Subscript[g, 2])])
 (OverVector[e]) .(OverVector[(Subscript[g, 1])]) =0       et      (OverVector[e]) .(OverVector[(Subscript[g, 2])]) =0

On obtient ainsi les équations normales

((OverVector[y])-a (OverVector[(Subscript[g, 1])])-b (OverVector[(Subscript[g, 2])])).(OverVector[(Subscript[g, 1])]) =0       et    ( (OverVector[y])-a (OverVector[(Subscript[g, 1])])-b (OverVector[(Subscript[g, 2])])).(OverVector[(Subscript[g, 2])]) =0

(OverVector[y]).(OverVector[(Subscript[g, 1])]) -a (OverVector[(Subscript[g, 1])]).(OverVector[(Subscript[g, 1])]) -b (OverVector[(Subscript[g, 2])]).(OverVector[(Subscript[g, 1])]) =0       et    (OverVector[y]).(OverVector[(Subscript[g, 2])]) -a (OverVector[(Subscript[g, 1])]).(OverVector[(Subscript[g, 2])]) -b (OverVector[(Subscript[g, 2])]) .(OverVector[(Subscript[g, 2])]) =0

(OverVector[(Subscript[g, 1])]).(OverVector[(Subscript[g, 1])]) a+(OverVector[(Subscript[g, 2])]).(OverVector[(Subscript[g, 1])]) b=(OverVector[y]).(OverVector[(Subscript[g, 1])])        et  (OverVector[(Subscript[g, 1])]).(OverVector[(Subscript[g, 2])]) a+(OverVector[(Subscript[g, 2])]).(OverVector[(Subscript[g, 2])]) b=  (OverVector[y]).(OverVector[(Subscript[g, 2])])

(({{(OverVector[(Subscript[g, 1])]).(OverVector[(Subscript[g, 1])]) , (OverVector[(Subscript[g, 2])]).(OverVector[(Subscript[g, 1])]) }, { (OverVector[(Subscript[g, 1])]).(OverVector[(Subscript[g, 2])]), (OverVector[(Subscript[g, 2])]).(OverVector[(Subscript[g, 2])])}}))(({{a}, {b}}))=(({{(OverVector[y]).(OverVector[(Subscript[g, 1])])  }, {(OverVector[y]).(OverVector[(Subscript[g, 2])]) }}))

Les coefficients du système sont

(OverVector[(Subscript[g, 1])]).(OverVector[(Subscript[g, 1])])=(Subscript[x, 1]^2)+(Subscript[x, 2]^2)+...+(Subscript[x, n]^2)=(∑_(i=1)^n)(Subscript[x, i]^2)

(OverVector[(Subscript[g, 2])]).(OverVector[(Subscript[g, 1])])=(Subscript[x, 1])+(Subscript[x, 2])+...+(Subscript[x, n])=(∑_(i=1)^n)(Subscript[x, i])=(OverVector[(Subscript[g, 1])]).(OverVector[(Subscript[g, 2])])

(OverVector[(Subscript[g, 2])]).(OverVector[(Subscript[g, 2])])=(1^2)+(1^2)+...+(1^2)=n

(OverVector[y]).(OverVector[(Subscript[g, 1])])=(Subscript[y, 1])(Subscript[x, 1])+(Subscript[y, 2])(Subscript[x, 2])+...+(Subscript[y, n])(Subscript[x, n])=(∑_(i=1)^n)(Subscript[x, i])(Subscript[y, i])

(OverVector[y]).(OverVector[(Subscript[g, 2])])=(Subscript[y, 1])+(Subscript[y, 2])+...+(Subscript[y, n])=(∑_(i=1)^n)(Subscript[y, i])

Les équations normales prennent la forme

(({{(∑_(i=1)^n)(Subscript[x, i]^2), (∑_(i=1)^n)(Subscript[x, i]) }, {(∑_(i=1)^n)(Subscript[x, i]), n}}))(({{a}, {b}}))=(({{(∑_(i=1)^n)(Subscript[x, i])(Subscript[y, i]) }, {(∑_(i=1)^n)(Subscript[y, i]) }}))

Exemple numérique et résolution des équations normales

x={1,3,4,6,8,9,11,14};
y={1,2,4,4,5,7,8,9};

<<Tableaux`;
afficheTableau[{"x","y"},None,{x,y}]

n=Length[x]

(Overscript[(Subscript[g, 1]), ⟶])=(Overscript[x, →])

g2=Table[1,{n}]

m={{x.x,x.g2},{g2.x,g2.g2}};MatrixForm[m]

v={y.x,y.g2}

{a,b}=LinearSolve[m,v]

Les valeurs de a et b peuvent aussi être estimées approximativement à partir d'un graphique:

Finalement, nous avons résolu le problème en 5 étapes:
  1.	formulation du problème d'ajustement;
  2.	réduction du problème d'ajustement à un problème de moindres carrés;
  3.	réduction du problème de moindres carrés à un problème de projection orthogonale;
  4.	réduction du problème de projection orthogonale à un système d'équations normales;
  5.	résolution du système d'équations.

Autre méthode de calcul : la fonction Fit[..] de Mathematica

?Fit

donnees=Transpose[{x,y}]

Comme la lettre x est déjà utilisée pour désigner la liste des abscisses, on doit choisir une autre lettre pour désigner la variable. Si on choisit comme variable le symbole t, les fonctions de base sont

(Subscript[f, 1])(t)=t,           (Subscript[f, 2])(t)=1

et la fonction que l'on cherche est une combinaison linéaire des deux fonctions de base

a t+b=a (Subscript[f, 1])(t)+b (Subscript[f, 2])(t)

f[t_]=Fit[donnees,{t,1},t]

En prenant un autre exemple,  Fit[donnees,{(t^2),(1/t)},t]   signifie que l'on combine les fonctions t ↦ (t^2)  et  t ↦ (1/t) pour obtenir la fonction d'ajustement  t↦ a·(t^2)+ b· (1/t)  dans laquelle a et b sont les deux inconnues du problème d'ajustement.

Remarque: pour définir la fonction, utiliser une assignation immédiate "=" plutôt qu'une assignation différée ":=".

Exercice 3.1 - 1            [Sans ordinateur] Propriétés des résidus

On ajuste la droite y=a x+b aux n points de mesure ((Subscript[x, 1]), (Subscript[y, 1])), ..., ((Subscript[x, n]), (Subscript[y, n])) au sens des moindres carrés. On note (Subscript[e, i])=(Subscript[y, i])-a (Subscript[x, i])-b les résidus pour i ∈ {1, 2, ..., n}. Démontrez que

a) 	la somme des résidus est nulle:
	(∑_(i=1)^n)(Subscript[e, i])=0

b) 	la somme des résidus, pondérée par les abscisses, est nulle:
	(∑_(i=1)^n)(Subscript[x, i])(Subscript[e, i])=0

Indication: Il s'agit de conséquences immédiates des équations normales.

Exercice 3.1 - 2   [Sans ordinateur] Ajustements à un paramètre

Reprenons les données de l'exemple numérique:

a) 	Ajustez la droite y=a x à ces données, au sens des moindres carrés.
	Directive: répétez toute la démarche en 5 étapes.

b) 	Ajustez la droite y=b à ces données, au sens des moindres carrés.
	Interprétez le résultat obtenu.

Exercice 3.1- 3      Ajustement d'une parabole

Voici une statistique des taux de natalité aux USA. Le taux de natalité est le nombre d'enfants nés vivants en une année pour mille habitants

x=Range[1915,1955,5];
y={25.0,23.7,21.3,18.9,16.9,17.9,19.5,23.6,24.6};

ListPlot[Transpose[{x,y}],Joined->True,AxesOrigin->{1910,15}]

On y remarque l'effet de la crise de 1929. Fixons l'origine du temps au premier juillet 1935

ListPlot[Transpose[{x,y}],Joined->True,AxesOrigin->{0,15}]

On demande d'ajuster la parabole y=a (x^2)+b x+c à ces données.
Directive: répétez toute la démarche en 5 étapes.
[Avec Mathematica] Faites une représentation graphique.

Exercice 3.1 - 4      Droite de régression de x en fonction de y

Parmi les étudiants de sexe masculin d'une université, on en a tiré 12 au hasard. On a mesuré leurs masses x en kilogrammes et leurs tailles respectives y en centimètres:

x={70,68,63,72,60,66,70,74,65,62,67,65};
y={155,152,150,180,135,156,168,178,160,132,145,139};

Ajustons d'abord la droite y=a x+b au sens des moindres carrés; cette droite est aussi appelée droite de régression de y en fonction de x; on obtient

Fit[Transpose[{x,y}],{t,1},t]

Ajustons ensuite la droite x=m y+p au sens des moindres carrés; cette droite est aussi appelée droite de régression de x en fonction de y; on obtient

Fit[Transpose[{y,x}],{t,1},t]

Questions:

a)	Pourquoi les deux fonctions affines ne sont-elles pas réciproques ?

b)	Ajustez la droite y=a x+b en minimisant la somme des carrés des résidus horizontaux.

§ 3.2	Moindres carrés pondérés             [Pour OS PAM seulement]

On peut donner plus d'importance à certaines mesures en les répétant. Reprenons l'exemple numérique du § 3.1 en répétant deux fois la première mesure et trois fois la dernière:

x={1,1,3,4,6,8,9,11,14,14,14};
y={1,1,2,4,4,5,7,8,9,9,9};

Alors, dans l'expression des moindres carrés apparaissent des coefficients appelés "poids":

2(((1-a 1-b))^2)+(((2-a 3-b))^2)+(((4-a 4-b))^2)+(((4-a 6-b))^2)+(((5-a 8-b))^2)+(((7-a 9-b))^2)+(((8-a 11-b))^2)+3(((9-a 14-b))^2)

Avec de tels poids, la réponse obtenue est différente:

donnees=Transpose[{x,y}];
Fit[donnees,{t,1},t]

Plus généralement, on peut définir un problème de "moindres carrés pondérés"

(Subscript[c, 1])(Subscript[e, 1]^2)+(Subscript[c, 2])(Subscript[e, 2]^2)+...+(Subscript[c, n])(Subscript[e, n]^2)=(Subscript[c, 1])((((Subscript[y, 1])-a (Subscript[x, 1])-b))^2)+(Subscript[c, 2])((((Subscript[y, 2])-a (Subscript[x, 2])-b))^2)+...+(Subscript[c, n])((((Subscript[y, n])-a (Subscript[x, n])-b))^2)

où les poids (Subscript[c, 1]), (Subscript[c, 2]), ..., (Subscript[c, n]) sont des nombres réels positifs (pas nécessairement entiers).
La méthode de résolution consiste à redéfinir un nouveau produit scalaire pondéré

⟨(Overscript[u, ->])|(Overscript[v, ->])⟩=(Subscript[c, 1])(Subscript[u, 1])(Subscript[v, 1])+(Subscript[c, 2])(Subscript[u, 2])(Subscript[v, 2])+...+(Subscript[c, n])(Subscript[u, n])(Subscript[v, n])

Avec ce nouveau produit scalaire pondéré, la norme des vecteurs ne sera plus la norme euclidienne usuelle mais sera aussi pondérée par les coefficients souhaités:

(((∥(Overscript[e, ->])∥))^2)=⟨(Overscript[e, ->])|(Overscript[e, ->])⟩=(Subscript[c, 1])(Subscript[e, 1]^2)+(Subscript[c, 2])(Subscript[e, 2]^2)+...+(Subscript[c, n])(Subscript[e, n]^2)

Pour que la norme de ce vecteur soit minimale, il suffit de calculer la projection orthogonale en remplaçant le produit scalaire usuel par le nouveau produit scalaire pondéré.

Dans notre exemple numérique,

x={1,3,4,6,8,9,11,14};
y={1,2,4,4,5,7,8,9};
c={2,1,1,1,1,1,1,3};
Clear[u,v,ps];
ps[u_List,v_List]:=(c*u).v

n=Length[x]

g2=Table[1,{n}]

m={{ps[x,x],ps[x,g2]},{ps[g2,x],ps[g2,g2]}};MatrixForm[m]

v={ps[y,x],ps[y,g2]}

{a,b}=LinearSolve[m,v];N[{a,b}]

Exercice 3.2-1          [Pour OS PAM seulement]

Une statistique a montré que les couples (x,y) ne prenaient que les valeurs données dans le tableau qui suit, avec les fréquences correspondantes c

x={0,1,2,3,4,5,6,7};
y={100,103,107,109,111,114,117,121};
c={0.05,0.1,0.15,0.2,0.2,0.15,0.1,0.05};

Ajustez la droite y=a x+b au sens des moindres carrés pondérés en prenant pour poids les fréquences (Subscript[c, i]).
Directive: répétez toute la démarche en 5 étapes.

§ 3.3	Ajustement par une combinaison linéaire de fonctions à plusieurs variables

Exemple:	Ajustement d'un plan

Le graphique d'une fonction de deux variables, par exemple

z=f(x,y)=-2 x-3 y+25

se présente comme une surface dans l'espace.

Plot3D[-2 x-3 y+25,{x,0,5},{y,0,5},ViewPoint->{3,1,1},AxesLabel->{"x","y","z"}]

Dans le cas particulier où la fonction est affine

z=f(x,y)=a x+b y+c

la surface est un plan dont l'équation est

a x+b y-z+c=0

Etant donné un "nuage de points" et un plan, on appelle "résidus" les écarts entre les points et le plan mesurés parallèlement à l'axe des z; on cherche le plan tel que la somme des carrés des résidus soit minimale.

Choisissons comme variables les symboles s et t. On peut considérer la fonction cherchée comme étant une combinaison linéaire des trois fonctions de base suivantes

(Subscript[f, 1])(s,t)=s,               (Subscript[f, 2])(s,t)=t,                 (Subscript[f, 3])(s,t)=1

En effet,

a (Subscript[f, 1])(s,t)+b (Subscript[f, 2])(s,t)+c (Subscript[f, 3])(s,t)=a s+b t+c

Exercice 3.3-1

On a choisi 12 garçons au hasard et on a mesuré
	x = l'âge, en années;
	y = le tour de taille, en cm;
	z = le tour de poitrine, en cm.

x={8,10,6,11,8,7,10,9,10,6,12,9};
y={57,59,49,62,51,50,55,48,52,42,61,57};
z={64,71,53,67,55,58,77,57,56,51,76,68};

Ajustez la fonction z=a x+b y+c aux données, au sens des moindres carrés.
  a)	 En répétant toute la démarche en 5 étapes.
  b)	Avec la méthode Fit[..] de Mathematica.

§ 3.4 Ajustements trigonométriques

L'ajustement de données périodiques fait généralement appel aux fonctions trigonométriques.

Exemple

Ajustez les paramètres  (Subscript[c, 0]), ..., (Subscript[c, 5]) de la fonction

f(x)=(Subscript[c , 0])+(Subscript[c, 1]) cos(x)+(Subscript[c, 2])sin(x)+(Subscript[c, 3]) cos(2x)+(Subscript[c, 4])sin(2x)+(Subscript[c, 5]) cos(3x)

aux données

{(0,0),(((2 π)/9),0.2),(2((2 π)/9),0.5),(3((2 π)/9),0.3),(4((2 π)/9),-0.2),(5((2 π)/9),-0.4),(6((2 π)/9),-0.1),(7((2 π)/9),0.1),(8((2 π)/9),-0.1)}

Dans un premier temps, résolvons ce problème comme un ajustement au sens des moindres carrés selon la méthode générale du § 3.1:

p={{0,0},{((2 π)/9),0.2},{2((2 π)/9),0.5},{3((2 π)/9),0.3},{4((2 π)/9),-0.2},{5((2 π)/9),-0.4},{6((2 π)/9),-0.1},{7((2 π)/9),0.1},{8((2 π)/9),-0.1}};

x=Transpose[p][[1]]

n=Length[x]

Les vecteurs de base de l'espace sur lequel on effectue la projection sont

g[0]=Table[1,{n}]

g[1]=Cos[x]

g[2]=Sin[x]

g[3]=Cos[2x]

g[4]=Sin[2x]

g[5]=Cos[3x]

k=6

Le système des équations normales est

m=Table[FullSimplify[g[i].g[j]],{j,0,k-1},{i,0,k-1}];MatrixForm[m]

y=Transpose[p][[2]]

v=Table[g[j].y,{j,0,k-1}]

c=LinearSolve[m,v]

La fonction d'ajustement est la meilleure approximation des données par un polynôme trigonométrique

f[t_]=c.{1,Cos[t],Sin[t],Cos[2t],Sin[2t], Cos[3t]}

Plot[f[t],{t,0,2 π},Epilog->{PointSize[0.02],Point/@p},PlotRange->1.2{Min[y],Max[y]}]

Première idée de la méthode "transformée de Fourier rapide"

Dans le cas où il y a k fonctions de base

1,cos(x),sin(x),cos(2x),sin(2x), cos(3x), sin(3x),cos(4x), ...

et n abscisses uniformément réparties sur une période

0,((2π)/n),2((2π)/n),3((2π)/n),...,(n-2)((2π)/n),(n-1)((2π)/n)

avec k≤n, alors la méthode d'ajustement donne une matrice diagonale. La résolution du système nécessite alors moins de calculs. En poursuivant l'effort d'optimisation, il est possible d'organiser les calculs de telle manière que, même appliquée à des données volumineuses, la méthode reste efficace. On parle alors de "transformée de Fourier rapide".

Généralisation

Fonctions de base de Fourier

Généralisons à des fonctions périodiques de période T. En posant ω=((2 π)/T), choisissons les k fonctions suivantes

(Subscript[f, 0])(t)=1,
(Subscript[f, 1])(t)= cos(ω t), (Subscript[f, 2])(t)= sin(ω t),
(Subscript[f, 3])(t)=cos(2 ω t), (Subscript[f, 4])(t)=sin(2 ω t), 
(Subscript[f, 5])(t)=cos(3 ω t), (Subscript[f, 6])(t)=sin(3 ω t),
(Subscript[f, 7])(t)=cos(4 ω t),...,
(Subscript[f, k-1])(t)=⋯

c'est-à-dire

(Subscript[f, 0])(t)=1,   (Subscript[f, j])(t)={({{cos((((j+1)/2))ω t), si j est impair}, {sin((j/2) ω t), si j est pair}}),   j=1,2,...,k-1

Echantillonnage ou discrétisation (ou "digitalisation")

Nous supposons que les abscisses d'échantillonnage sont uniformément réparties sur une période. Plus précisément, dans l'intervalle de référence  [0,  T], on choisit

(Subscript[t, j])=j(T/n)  ,      j=0, 1,... , n-1

On suppose que k ≤ n.

En ces abscisses, on donne les valeurs correspondantes

(Subscript[y, 0]),(Subscript[y, 1]),(Subscript[y, 2]),...,(Subscript[y, n-1])

Vecteurs de base

Les vecteurs de base de l'espace vectoriel sur lequel on effectue la projection orthogonale sont  { (OverVector[(Subscript[g, 0])]),(OverVector[ (Subscript[g, 1])]), ..., (OverVector[(Subscript[g, k-1])])}

(OverVector[(Subscript[g, 0])])=(({{1}, {1}, {1}, {1}, {...}, {1}})),(OverVector[(Subscript[g, 1])])=(({{cos(0)}, {cos(((2 π)/n))}, {cos(2((2 π)/n))}, {cos(3((2 π)/n))}, {...}, {cos((((n-1)2 π)/n))}})),..,(OverVector[(Subscript[g, 3])])=(({{cos(0)}, {cos(((4 π)/n))}, {cos(2((4 π)/n))}, {cos(3((4 π)/n))}, {...}, {cos((((n-1)4 π)/n))}})),..,(OverVector[(Subscript[g, k-1])])=(({{(Subscript[f, k-1])((Subscript[t, 0]))}, {(Subscript[f, k-1])((Subscript[t, 1]))}, {(Subscript[f, k-1])((Subscript[t, 2]))}, {(Subscript[f, k-1])((Subscript[t, 3]))}, {...}, {(Subscript[f, k-1])((Subscript[t, n-1]))}}))

Théorème 3.4

Avec les choix effectués ci-dessus, les vecteurs { (OverVector[(Subscript[g, 0])]), (OverVector[(Subscript[g, 1])]), ...,(OverVector[ (Subscript[g, k-1])])} sont orthogonaux, c'est-à-dire

(OverVector[(Subscript[g, i])]).(OverVector[(Subscript[g, j])])=0        pour       i≠j

En d'autres termes, les équations normales sont un système diagonal. (Voir exercices).

Corollaire

Pour ajuster la fonction

f(t)=(Subscript[c, 0])(Subscript[f, 0])(t)+(Subscript[c, 1])(Subscript[f, 1])(t)+...+(Subscript[c, k-1])(Subscript[f, k-1]) (t)

aux données, le calcul des coefficients peut s'effectuer simplement au moyen de la formule

(Subscript[c, j])=(((OverVector[y]).(OverVector[(Subscript[g, j])]))/((OverVector[(Subscript[g, j])]).(OverVector[(Subscript[g, j])]))),    j=0, 1, ..., k-1

Ces coefficients sont appelés "coefficents de Fourier".
La fonction d'ajustement est périodique de période T .
Les résidus suivants sont minimaux au sens des moindres carrés:

(Subscript[e, 0])=(Subscript[y, 0])-f((Subscript[t, 0])), ...,(Subscript[e, k-1])=(Subscript[y, k-1])-f((Subscript[t, k-1]))

Application

Dans le cas où le nombre de points d'échantillonnage est élevé, la résolution des équations normales peut prendre un temps excessif. On peut éviter de former la matrice m, le vecteur v et d'exécuter LinearSolve[m, v]. En lieu et place, on utilise la forme explicite de la réponse:

c=Table[((y.g[j])/(g[j].g[j])),{j,0,k-1}]

Il est encore possible d'améliorer la vitesse de calcul (par exemple, on peut donner les valeurs de g[i].g[j]) mais cela sortirait du cadre de ce cours. La méthode de calcul optimalisée s'appelle "transformée de Fourier rapide" (en anglais "Fast Fourier Transform" abrégé FFT).

Cas particulier de l'interpolation

Dans le cas où k=n, les résidus sont nuls et la fonction d'ajustement passe par les points donnés:

(Subscript[y, 0])=f((Subscript[t, 0])), ...,(Subscript[y, n-1])=f((Subscript[t, n-1]))

Au lieu d'ajustement, on parle alors d'interpolation. Un chapitre entier est consacré à ce thème.

Exercice 3.4 - 1

Voici le tableau des températures diurnes moyennes collectées en un lieu imaginaire

A ces données, ajustez les paramètres de la fonction

f(t)=(Subscript[c , 0])+(Subscript[c, 1]) cos(ω t)+(Subscript[c, 2])sin(ω t)

Utilisez la méthode de calcul "rapide".

Dans un même graphique, superposez les données et la fonction ajustée.

Indications:
Pour le temps, on prendra le mois comme unité. Nous aurons donc

T=12;                 ω=((2 π)/12)

Pour respecter les abscisses imposées par la méthode, nous prendrons t=0 à la mi-janvier;...; t=11 à la mi-décembre.

Exercice 3.4 - 2

Ecrivez explicitement les assertions du théorème 3.4 sous la forme d'énoncés trigonométriques. Par exemple,

(OverVector[(Subscript[g, 0])])·(OverVector[(Subscript[g, 1])])=(∑_(i=0)^(n-1))cos(i((2 π)/n))=0

Exercice 3.4 - 3             [Facultatif]

Démontrez les relations suivantes

(OverVector[(Subscript[g, 0])])·(OverVector[(Subscript[g, 1])])=(∑_(i=0)^(n-1))cos(i·1·((2 π)/n))=0

(OverVector[(Subscript[g, 0])])·(OverVector[(Subscript[g, 2])])=(∑_(i=0)^(n-1))sin(i·1·((2 π)/n))=0

(OverVector[(Subscript[g, 0])])·(OverVector[(Subscript[g, 3])])=(∑_(i=0)^(n-1))cos(i·2·((2 π)/n))=0

(OverVector[(Subscript[g, 0])])·(OverVector[(Subscript[g, 4])])=(∑_(i=0)^(n-1))sin(i·2·((2 π)/n))=0

(OverVector[(Subscript[g, 0])])·(OverVector[(Subscript[g, 5])])=(∑_(i=0)^(n-1))cos(i·3·((2 π)/n))=0,           etc.

Indications
Le point O(0,0) désignant l'origine, définissons les points (Subscript[E, i]^j) tels que

(Overscript[(Subscript[OE, i]^j), ⟶])=(({{cos(i·j·((2π)/n))}, {sin(i·j·((2π)/n))}}))

En prenant les exemples j=1, j=2,j=3, représentez graphiquement et vérifiez géométriquement que

(∑_(i=0)^(n-1))  (Overscript[(Subscript[OE, i]^j), ⟶])  =(Underoverscript[ ∑, (j=0), (n-1)]) (({{cos(i·j·((2π)/n))}, {sin(i·j·((2π)/n))}}))=(OverVector[0])

Déduisez que, pour chaque j,

(∑_(i=0)^(n-1))cos(i·j·((2π)/n))=0          et           (∑_(i=0)^(n-1))sin(i·j·((2π)/n))=0

Exercice 3.4 - 4             [Facultatif]

Démontrez que les vecteurs du théorème 3.4

(OverVector[(Subscript[g, 0])]), (OverVector[(Subscript[g, 1])]),..., (OverVector[(Subscript[g, k-1])])

sont orthogonaux.

Indications
Utilisez les résultats de l'exercice 3.4-3 et les relations trigonométriques

cos(x)cos(y)=(1/2)(cos(x+y)+cos(x-y))

cos(x)sin(y)=(1/2)(sin(x+y)-sin(x-y))

sin(x)sin(y)=(1/2)(-cos(x+y)+cos(x-y))

Exercice 3.4 - 5            [Travail dirigé]

Dans la région de Zurich, la température moyenne mensuelle en degrés Celsius est la suivante:

Question a)
	Ajustez la famille de fonctions

f(t)=a+b cos( ω t)+c sin(ω t)

à ces températures.

Indications pour la question a)
	1)	Pour ajuster les 12 points donnés par une simple sinusoïde, nous ne prenons que les 3 premières
		fonctions de base données dans la théorie : 1, cos(((2 π)/365)t), sin(((2 π)/365)t).
	2)	Associer les températures au 15 de chaque mois ne convient qu'approximativement car
		les mois ont des durées inégales; par ailleurs, l'ajustement se calcule plus aisément pour la suite des 
		abscisses 0, (T/12), ((2 T)/(12,)), ..., ((11 T)/12); divisons l'année de 365 jours en 12 intervalles égaux

x=Range[0,360,(365/12)]

ce qui nous donne les points d'interpolation suivants

3)	pour interpréter les résultats du calcul, on se réfère au tableau suivant;
		il s'agit d'associer les dates voulues à la numérotation des jours par la variable  t (ou  x)

xr=Round[x]

tb=Transpose[{xr,{"16 janvier","15 février","18 mars","17 avril","18 mai","17 juin","17 juillet","17 août","16 septembre","17 octobre","16 novembre","17 décembre"}}];tableauGraph[tb];

Question b)
	Ajustez la fonction

f(t)=a+d cos( ω t-φ)

à ces températures.

Indications pour la question b)
Le point de départ est le résultat de la partie a)
f(t)=a+b cos(ω t)+c sin(ω t)  que nous allons mettre sous la forme désirée. Pour ce faire, calculons les coordonnées polaires du vecteur suivant

(({{b}, {c}}))=d(({{cos(φ)}, {sin(φ)}}))

où  d   est la norme du vecteur (d est dénommé amplitude)

d=(Sqrt[(b^2)+(c^2)])

φ désigne l'angle entre les vecteurs (({{1}, {0}}))     et      (({{b}, {c}})):

cos(φ)=(b/d),              sin(φ)=(c/d)

L'angle  φ,  dénommé déphasage, peut s'écrire

φ=sign(c)·arccos((b/d))

En utilisant la relation trigonométrique (voir Formulaires et tables)

cos(α-β)=cos(α)cos(β)+sin(α)sin(β)

la fonction ajustée peut être transformée sous la forme souhaitée

f(t)=a+d( cos( ω t)(b/d)+ sin(ω t)(c/d))=a+d( cos( ω t)cos(φ)+ sin(ω t)sin(φ))=a+d cos( ω t-φ)

Question c)
	Donnez une interprétation physique et graphique des grandeurs suivantes
	1)	a= température annuelle moyenne (en degrés Celsius);
	2)	d= amplitude thermique annuelle (en degrés Celsius);
	3)	(φ/ω)= décalage temporel (en jours);
	4)	φ= déphasage (en radians).