Cómo ajustar un conjunto de datos a una función cualquiera

Tenemos un conjunto de puntos que queremos ajustar a una cierta función. Por ejemplo, tenemos los puntos:

x = [1 2 3 4 5]';
y = [1.24 0.59 0.31 0.19 0.09]';

Y queremos ajustarlos a la función y = k·exp(-a·x). Para ello, tenemos que estimar los parámetros "k" y "a".

---

Primero, creamos la función a ajustar en un archivo .m, por ejemplo "f01.m":

function z = f01(u,x)
z = u(1)*exp(-u(2)*x);

donde u(1) = k, y u(2) = a .

A continuación, creamos la función objetivo que vamos a optimizar, y que sume los errores al cuadro entre los valores medidos y los calculados. La siguiente función la creamos en el archivo "fobj.m":

function z = fobj(u,x,ym)

ycal = f01(u,x);
z = sum((ym-ycal).^2);

Por último, calculamos los parámetros "k" y "a" con la función "fminsearch":

[u,f] = fminsearch(@(u) fobj(u,x,y),[1 1])

donde [1 1] es el valor inicial de los parámetros, con los que el optimizador empezará a buscar la solución, la primera salida de la función "u" es un vector con los parámetros estimados, y la segunda salida de la función "f" es el valor de la función objetivo en el punto mínimo.

Matlab nos devuelve el siguiente resultado:

u =

    2.4430    0.6864

f =

    0.0022

---

En vez de usar "fminsearch" y tener que crear una función objetivo, podríamos haber usado directamente la funcion "lsqcurvefit". Pero la verdad, es que el método anterior me suele funcionar mejor.

[u,f] = lsqcurvefit(@(u,x) f01(u,x),[1,1],x,y)

donde [1 1] son los valores iniciales de u(1) y u(2), con los que la función empieza a tantear. Por cierto, cuando escribimos @(u,x), hay que poner primero el vector de parámetros a estimar, y segundo la variable independiente.

---

Por último, representamos tanto los puntos reales como la curva obtenida:

xx = [0:0.1:6];
yy = f01(u,xx);

figure(1)
hold on
box on
plot(xx,yy,'-','Color',[1 0.6 0.78],'LineWidth',2)
plot(x,y,'o','MarkerSize',5,'Color',[0.48 0.06 0.89],'MarkerFaceColor',[0.48 0.06 0.89])
legend({'real','estimado'},1,'EdgeColor',[1 1 1]);
axis([0 6 0 3])
set(gca,'XTick',[0:2:6])
set(gca,'YTick',[0:1:3])
set(gca,'Fontsize',10)

Nota: He usado 'EdgeColor',[1 1 1] para que la leyenda no esté encerrada en una caja. Y 'MarkerFaceColor',[0.48 0.06 0.89] para que los círculos estén rellenos.

Cómo trasladar y girar una curva

Mover una curva es otra de las cosas más sencillas de hacer. Por ejemplo, si tenemos la siguiente curva que parte del punto (0,0) y queremos moverla al punto (6,4), podemos hacer la siguiente transformación:

% curva original
x = [0:0.1:4*pi];
y = sin(x);

% translación al punto (6,4)
x1 = x + 6;
y1 = y + 4;

figure(1)
hold on
box on
plot(x,y,'-b')
plot(x1,y1,'-m')
axis equal
legend('original','movida',2)

---

Si lo que queremos es girar la curva un ángulo tita, la operación también es muy sencilla. Aplicamos la siguiente transformación (que se calcula muy facilmente por geometría). Por ejemplo, para un ángulo θ = τ/6 radianes:

% rotacion de tita radianes
tita = 2*pi/6;
x1 = (x-y*tan(tita))*cos(tita);
y1 = y/cos(tita) + (x-y*tan(tita))*sin(tita);

figure(2)
hold on
box on
plot(x,y,'-b')
plot(x1,y1,'-m')
axis equal
legend('original','girada')

Hay que tener cuidado, porque las ecuaciones anteriores no están definidas si θ vale τ/4 (π/2) o múltiplos. En esos casos, hay que girar la curva directamente de la forma:

% giro de 1/4 tau (1/2 pi) radianes
x1 = -y;
y1 = x;

% giro de 2/4 tau (pi) radianes
x2 = -x;
y2 = -y;

% giro de 3/4 tau (3/2 pi) radianes
x3 = y;
y3 = -x;

figure(3)
hold on
box on
plot(x,y,'-b')
plot(x1,y1,'-m')
plot(x2,y2,'-r')
plot(x3,y3,'-k')
legend('original','girada 1/4 \tau','girada 2/4 \tau', 'girada 3/4 \tau')

Cómo dibujar círculos y polígonos regulares

Para dibujar un círculo con Matlab, una posibilidad muy fácil, es hacerlo en coordenadas polares. Es decir:  

R = 1;
tita = (0:0.01:2.01*pi);
x = R*cos(tita);
y = R*sin(tita);

figure(1)
plot(x,y,'-b')
axis equal
axis([-1.5 1.5 -1.5 1.5])
set(gca,'XTick',[-1.5:0.5:1.5])
set(gca,'XTickLabel',{'-1.5','-1.0','-0.5','0.0','0.5','1.0','1.5'})
set(gca,'YTick',[-1.5:0.5:1.5])
set(gca,'YTickLabel',{'-1.5','-1.0','-0.5','0.0','0.5','1.0','1.5'})
set(gca,'Fontsize',10)

Por cierto, el valor de "tita" lo llevo hasta 2.01*pi, en vez de hasta solo 2.00*pi, para que el círculo se cierre (y no quede un pequeñito espacio en blanco). Además, he puesto "axis equal" para evitar que Matlab distorsione los ejes.

En realidad, no hemos dibujo un círculo, sino un polígono de 630 lados. Lo que a simple vista, es suficientemente parecido a un círculo.

Bueno, pues entonces, dibujar un polígono regular de n lados, se hace igual. Por ejemplo, para un pentágono:

n = 5;
R = 1;
tita = [0:(2*pi/n):2*pi];
x = R*cos(tita);
y = R*sin(tita);

figure(1)
plot(x,y,'-b')
axis equal
axis([-1.5 1.5 -1.5 1.5])
set(gca,'XTick',[-1.5:0.5:1.5])
set(gca,'XTickLabel',{'-1.5','-1.0','-0.5','0.0','0.5','1.0','1.5'})
set(gca,'YTick',[-1.5:0.5:1.5])
set(gca,'YTickLabel',{'-1.5','-1.0','-0.5','0.0','0.5','1.0','1.5'})
set(gca,'Fontsize',10)

Para evitar que el pentágono esté girado, podemos hacer que el valor de "tita" empiece en π/2:

n = 5;
R = 1;
tita = [0:(2*pi/n):2*pi]+pi/2;
x = R*cos(tita);
y = R*sin(tita);

Ahora ya es muy fácil hacer muchas figuras. Por ejemplo, una estrella de 5 puntas:

n = 5;
R = 1;
tita = [0:(4*pi/n):4*pi]+pi/2;
x = R*cos(tita);
y = R*sin(tita);

Cómo trabajar con coordenadas polares

Muchas funciones son más fáciles de ser representadas en coordenadas polares. Por ejemplo, la cardiode:

ρ = R·(1+cos(θ+τ/4))

donde θ va de 0 a τ. (Ya sabes, τ = 2·π.)

---

Para dibujar una función en polares, podemos usar directamente la función "polar", como muestro a continuación:

R = 1;
tita = 0:0.01:2.0*pi;
rho = R*(1+cos(tita+pi/2));

figure(1)
polar(tita,rho,'m-')

---

Sin embargo, para editar el gráfico, es más fácil hacerlo en coordenadas cartesianas. Para ello, primero transformamos las variables "tita" y "rho", en cartesianas ("x" e "y"):

x = rho.*cos(tita);
y = rho.*sin(tita);

Y a continuación, representamos "x" e "y" normalmente:

figure(1)
hold on
box on
plot(x,y,'-','LineWidth',2,'Color',[0.5 0.05 0.9])
axis equal
axis([-1.5 1.5 -2.5 0.5])
set(gca,'XTick',[-1.5:0.5:1.5])
set(gca,'XTickLabel',{'-1.5','-1.0','-0.5','0.0','0.5','1.0','1.5'})
set(gca,'YTick',[-2.5:0.5:0.5])
set(gca,'YTickLabel',{'-2.5','-2.0','-1.5','-1.0','-0.5','0.0','0.5'})
set(gca,'Fontsize',10)

Nota: La línea "axis equal" es importante para que Matlab no deforme el tamaño de los ejes, como suele hacer por defecto

Por último, si queremos añadir los ejes, podemos hacerlo directamente así:

plot([-2 2],[0 0],'--','Color',[0.5 0.5 0.5])
plot([0 0],[-3 2],'--','Color',[0.5 0.5 0.5])

Cómo ajustar un conjunto de puntos a un polinomio

Pocas cosas son más fáciles que ajustar un conjunto de puntos a un polinomio. Por ejemplo, tenemos la siguiente serie de puntos, que queremos ajustar a un polinomio de orden dos:

x = [1 2 3 4 5];
y = [3.6 7.1 9.1 11.3 12.5];

Para ello, usamos la función "polyfit":

n = 2;
alfa = polyfit(x,y,n)

alfa es un vector de tamaño 3x1, formado por los coeficientes del polinomio:

y = alfa[1]·x² + alfa[2]·x + alfa[3]

En este caso, alfa vale:

alfa =

   -0.3143    4.0857   -0.0800

---

Para representar gráficamente el polinomio correspondiente, podemos usar la función "polyval", como muestro a continuación:

x1 = [0:0.1:6];
y1 = polyval(alfa,x1);

figure(1)
hold on
box on
plot(x,y,'or')
plot(x1,y1,'-b')
axis([0 6 0 15])
legend('real','estimado',2)