## Ajuste de datos a funciones armónicas

Vamos a emplear el procedimiento de mínimos cuadrados para ajustar un conjunto de datos (ti,xi) a la función

x= A0+A1cos(ωt)+ B1sin(ωt)

El criterio de ajuste se toma como aquél en el que la desviación cuadrática media sea mínima, es decir, debe de ser mínima la suma

$\begin{array}{l}S=\sum _{i=1}^{N}{\left\{{x}_{i}-\left({A}_{0}+{A}_{1}\mathrm{cos}\left(\omega {t}_{i}\right)+{B}_{1}\mathrm{sin}\left(\omega {t}_{i}\right)\right)\right\}}^{2}\\ \frac{\partial S}{\partial {A}_{0}}=0\text{ }N·{A}_{0}+\left(\sum _{i=1}^{N}\mathrm{cos}\left(\omega {t}_{i}\right)\right){A}_{1}+\left(\sum _{i=1}^{N}\mathrm{sin}\left(\omega {t}_{i}\right)\right){B}_{1}=\sum _{i=1}^{N}{x}_{i}\\ \frac{\partial S}{\partial {A}_{1}}=0\text{ }\left(\sum _{i=1}^{N}\mathrm{cos}\left(\omega {t}_{i}\right)\right){A}_{0}+\left(\sum _{i=1}^{N}{\mathrm{cos}}^{2}\left(\omega {t}_{i}\right)\right){A}_{1}+\left(\sum _{i=1}^{N}\mathrm{sin}\left(\omega {t}_{i}\right)\mathrm{cos}\left(\omega {t}_{i}\right)\right){B}_{1}=\sum _{i=1}^{N}{x}_{i}\mathrm{cos}\left(\omega {t}_{i}\right)\\ \frac{\partial S}{\partial {B}_{1}}=0\text{ }\left(\sum _{i=1}^{N}\mathrm{sin}\left(\omega {t}_{i}\right)\right){A}_{0}+\left(\sum _{i=1}^{N}\mathrm{cos}\left(\omega {t}_{i}\right)\mathrm{sin}\left(\omega {t}_{i}\right)\right){A}_{1}+\left(\sum _{1}^{N}{\mathrm{sin}}^{2}\left(\omega {t}_{i}\right)\right){B}_{1}=\sum _{i=1}^{N}{x}_{i}\mathrm{sin}\left(\omega {t}_{i}\right)\\ \left(\begin{array}{ccc}N& \sum _{i=1}^{N}\mathrm{cos}\left(\omega {t}_{i}\right)& \sum _{i=1}^{N}\mathrm{sin}\left(\omega {t}_{i}\right)\\ \left(\sum _{i=1}^{N}\mathrm{cos}\left(\omega {t}_{i}\right)\right)& \left(\sum _{i=1}^{N}{\mathrm{cos}}^{2}\left(\omega {t}_{i}\right)\right)& \sum _{i=1}^{N}\mathrm{sin}\left(\omega {t}_{i}\right)\mathrm{cos}\left(\omega {t}_{i}\right)\\ \sum _{i=1}^{N}\mathrm{sin}\left(\omega {t}_{i}\right)& \sum _{i=1}^{N}\mathrm{cos}\left(\omega {t}_{i}\right)\mathrm{sin}\left(\omega {t}_{i}\right)& \sum _{i=1}^{N}{\mathrm{sin}}^{2}\left(\omega {t}_{i}\right)\end{array}\right)\left(\begin{array}{c}{A}_{0}\\ {A}_{1}\\ {B}_{1}\end{array}\right)=\left(\begin{array}{c}\sum _{i=1}^{N}{x}_{i}\\ \sum _{i=1}^{N}{x}_{i}\mathrm{cos}\left(\omega {t}_{i}\right)\\ \sum _{i=1}^{N}{x}_{i}\mathrm{sin}\left(\omega {t}_{i}\right)\end{array}\right)\end{array}$

Resolvemos el sistema de tres ecuaciones con tres incógnitas mediante el operador \ división a la izquierda para obtener el valor de los coeficientes A0, A1, B1.

Ejemplo:

La temperatura media en cada uno de los meses del año en un cierto lugar ha sido la siguiente

 Mes Temperatura Mes Temperatura Enero Febrero Marzo Abril Mayo Junio Julio Agosto 18.9 21.1 23.3 27.8 32.2 37.2 38.0 36.1 Septiembre Octubre Noviembre Diciembre 34.4 29.4 23.3 18.9

El periodo P=12 meses, la frecuencia angular ω=2π/360 rad/dia.

Para crear una tabla (ti,xi), asignamos el valor de la temperatura media del mes al día correspondiente a la mitad del mes, suponiendo que todos los meses tengan los mismos días.

t=(1:12)-0.5;
x=[18.9 21.1 23.3 27.8 32.2 37.2 38.0 36.1 34.4 29.4 23.3 18.9];
w=2*pi/12;

M=[length(t), sum(cos(w*t)), sum(sin(w*t))
sum(cos(w*t)), sum(cos(w*t).^2), sum(sin(w*t).*cos(w*t))
sum(sin(w*t)),sum(sin(w*t).*cos(w*t)), sum(sin(w*t).^2)]

C=[sum(x);sum(x.*cos(w*t));sum(x.*sin(w*t))];
A=M\C

hold on
plot(t,x,'ro','markersize',4,'markerfacecolor','r')
y=@(t) A(1)+A(2)*cos(w*t)+A(3)*sin(w*t);
fplot(y,[t(1),t(end)]);
xlabel('mes')
ylabel('temperatura °C')
title('Ajuste a funciones armónicas')
grid on
hold off

Corremos el script en la ventana de comandos

    M =
12.0000   -0.0000    0.0000
-0.0000    6.0000    0.0000
0.0000    0.0000    6.0000
A =
28.3878
-9.2559
-2.8002

La función que mejor ajusta a los datos es

$x=28.3878-9.2559\mathrm{cos}\left(\frac{\pi }{6}t\right)-2.8002\mathrm{sin}\left(\frac{\pi }{6}t\right)$

Observamos que los elementos de la matriz de los coeficientes del sistema de tres ecuaciones con tres incógnitas M son nulos excepto los de la diagonal principal. Esto se debe a las propiedades de las funciones armónicas (véase la figura más abajo). Cuando realizamos N=12 observaciones equiespaciadas tn=n-0.5, n=1,2...N. separadas Δt=P/N, siendo P=12 el periodo, se verifica que

$\begin{array}{l}\sum _{i=1}^{N}\mathrm{sin}\left(\omega {t}_{i}\right)=0\text{ }\sum _{i=1}^{N}\mathrm{cos}\left(\omega {t}_{i}\right)=0\\ \frac{\sum _{i=1}^{N}{\mathrm{sin}}^{2}\left(\omega {t}_{i}\right)}{N}=\frac{1}{2}\text{ }\frac{\sum _{i=1}^{N}{\mathrm{cos}}^{2}\left(\omega {t}_{i}\right)}{N}=\frac{1}{2}\\ \sum _{i=1}^{N}\mathrm{sin}\left(\omega {t}_{i}\right)\mathrm{cos}\left(\omega {t}_{i}\right)=0\text{ }\end{array}$

Para demostrarlo, calculamos la suma de la progresión geométrica

$\begin{array}{l}\sum _{n=1}^{12}\mathrm{exp}\left(i\frac{2\pi }{12}\left(n-\frac{1}{2}\right)\right)=\mathrm{exp}\left(-i\frac{2\pi }{6}\right)\sum _{n=1}^{12}\mathrm{exp}\left(i\frac{2\pi }{12}n\right)=\\ \mathrm{exp}\left(-i\frac{2\pi }{6}\right)\frac{1-{\left(\mathrm{exp}\left(i\frac{2\pi }{12}\right)\right)}^{12}}{1-\mathrm{exp}\left(i\frac{2\pi }{12}\right)}=\mathrm{exp}\left(-i\frac{2\pi }{6}\right)\frac{1-\mathrm{exp}\left(2\pi i\right)}{1-\mathrm{exp}\left(i\frac{2\pi }{12}\right)}=0\end{array}$

>> t=(1:12)-0.5;
>> sum(sin(2*pi*t/12))
ans =  -3.8858e-16
>> sum(cos(2*pi*t/12))
ans =  -6.6613e-16
>> sum(sin(2*pi*t/12).*cos(2*pi*t/12))
ans =  -3.8858e-16

Este resultado se puede generalizar para N medidas igualmente espaciadas en el periodo P

Despejamos los coeficientes A0, A1, B1en el sistema de tres ecuaciones con tres incógnitas

$\begin{array}{l}\left(\begin{array}{ccc}N& 0& 0\\ 0& N/2& 0\\ 0& 0& N/2\end{array}\right)\left(\begin{array}{c}{A}_{0}\\ {A}_{1}\\ {B}_{1}\end{array}\right)=\left(\begin{array}{c}\sum _{i=1}^{N}{x}_{i}\\ \sum _{i=1}^{N}{x}_{i}\mathrm{cos}\left(\omega {t}_{i}\right)\\ \sum _{i=1}^{N}{x}_{i}\mathrm{sin}\left(\omega {t}_{i}\right)\end{array}\right)\\ {A}_{0}=\frac{\sum _{i=1}^{N}{x}_{i}}{N}\text{ }\\ {A}_{1}=\frac{2}{N}\sum _{i=1}^{N}{x}_{i}\mathrm{cos}\left(\omega {t}_{i}\right)\text{ }{B}_{1}=\frac{2}{N}\sum _{i=1}^{N}{x}_{i}\mathrm{sin}\left(\omega {t}_{i}\right)\end{array}$

El análisis anterior se puede extender a un modelo general

x= A0+A1cos(ωt)+ B1sin(ωt)+ A2cos(2ωt)+ B2sin(2ωt)+...+ Amcos(mωt)+ Bmsin(mωt)

Para datos igualmente espaciados, los coeficientes valen

$\begin{array}{l}{A}_{0}=\frac{\sum _{i=1}^{N}{x}_{i}}{N}\text{ }\\ {A}_{k}=\frac{2}{N}\sum _{i=1}^{N}{x}_{i}\mathrm{cos}\left(k\omega {t}_{i}\right)\text{ }{B}_{k}=\frac{2}{N}\sum _{i=1}^{N}{x}_{i}\mathrm{sin}\left(k\omega {t}_{i}\right)\text{ }k=1,2,...m\end{array}$

### Posición del asteriode Pallas

Gauss estaba interesado en calcular de un modo preciso la órbita del asteroide Pallas, a partir de las observaciones de su posición.

 Ascensión Declinación 0 30 60 90 120 150 180 210 240 270 300 330 408 89 -66 10 338 807 1238 1511 1583 1462 1183 804

Disponía de 12 pares de datos (xi,yi) que deseaba interpolar con seis funciones armónicas.

y= A0+A1cos(ωx)+ B1sin(ωx)+ A2cos(2ωx)+ B2sin(2ωx)+...+ A6cos(6ωx)+ B6sin(6ωx)

x=0:30:330;
y=[408 89 -66 10 338 807 1238 1511 1583 1462 1183 804];
hold on
plot(x,y,'ro','markersize',4,'markerfacecolor','r')

n=length(y);
w=2*pi/360;
a0=sum(y)/n;
k=1:6;
a=cos(w*k'*x)*y'*2/n;
b=sin(w*k'*x)*y'*2/n;

x=0:360;
y=a0+a'*cos(w*k'*x)+b'*sin(w*k'*x);
plot(x,y)

xlim([0 360])
ylabel('Declinación (minutos)')
title('Posición del asteroide Pallas')
grid on
hold off

### Ejercicios

Ajustar un conjunto de datos (ti,xi) a la función

x= A0+A1cos(ωt)+ B1sin(ωt)

 Mes Iluminación Mes Iluminación Enero Febrero Marzo Abril Mayo Junio Julio Agosto 144 188 245 311 351 359 308 287 Septiembre Octubre Noviembre Diciembre 260 211 159 131

ω=2π/12

 t x 0 0.15 0.3 0.45 0.6 0.75 0.9 1.05 1.2 1.35 2.2 1.595 1.031 0.722 0.786 1.2 1.805 2.369 2.678 2.614

ω=2π/1.5