Las funciones que que aparecen en este capítulo no son necesariamente parte de un toolkit de Matlab. No es más que una clasificación artificial de las funciones con una finalidad concreta que no justifica la creación de un capítulo a parte. Esta sección es mucho menos amplia de lo que debería, se ha sacrificado a propósito a favor de temas más básicos como el cálculo el álgebra o el dibujo de gráficas. Mientras todo lo escrito hasta aquí es casi definitivo este capítulo no se considerará nunca como terminado.
6.1 Estadística descriptiva y análisis de datos
Matlab es un lenguaje muy utilizado en análisis de datos, tanto experimentales como producidos por otros programas. En el caso de las series de datos experimentales habrá que importar un archivo almacenado en un formato compatible y asignar todos los valores a una variable; más adelante aprenderemos a hacerlo.
Supongamos el caso ideal de que ya tenemos todos los datos cargados en una variable. Matlab tiene una cantidad enorme de rutinas que facilitan el análisis de datos de modo interactivo. las funciones son tan de alto nivel que es posible sacar cualquier resultado secundario sin necesidad de pensar un script; directamente con la consola. las funciones más útiles son:
- mean
- Calcula la media de una muestra de datos: x=1/n∑n=1nxi.
- std
- Calcula la desviación típica de una muestra de datos: s=1/n-1∑i=1n(xi-x)2.
- median
- Calcula la mediana de la muestra.
- min
- Valor mínimo de la muestra.
- max
- Valor máximo de la muestra.
- sort
- Ordena los elementos de menor a mayor.
- center
- Resta el valor medio de la serie de datos.
6.1.1 Ajuste de curvas por mínimos cuadrados.
Supongamos que tenemos una serie de datos obtenidos mediante un experimento y queremos hacer cálculos con ellos. El principal impedimento es que ya no disponemos de una serie continua de datos sino de medidas discretas. Las dificultades que se presentan son numerosas; no podemos poner dichos datos en un sistema de ecuaciones no lineales porque no es diferenciable, tampoco en una ecuación diferencial. La solución es convertir estos datos en una muestra contínua que nuestras herramientas numéricas puedan manejar.
La práctica más común es la de utilizar un ajuste polinómico. El polinomio resultante se obtiene resolviendo el problema de mínimos cuadrados. El grado del polinomio es uno de los datos que debemos escoger. Suele ser una regla válida que polinomios de mayor órden consiguen menor error pero no es siempre cierta. Es bastante usual utilizar erroneamente el ajuste por mínimos cuadrados; sirven para crear un modelo polinómico de una función, no para conseguir una curva a través de unos puntos; para esto está la interpolación. Por ejemplo, supongamos que disponemos de la siguiente serie de datos:
- polyfit
- Devuelve los coeficientes del polinomio del grado que queramos y soluciona el problema de ajuste por mínimos cuadrados.
Introducimos ahora las dos series de datos e intentamos ajustarlo por mínimos cuadrados para crear un modelo lineal de los datos adquiridos. Lo conseguimos con la siguiente líneas de código: fit=polyfit(linspace(0,1,11),y,1)
fit=
1.5928
1.9828
Ahora tenemos dos números en un vector llamado fit, este vector son los dos coeficientes del polinomio que ajusta los datos; una recta. Ahora representaremos gráficamente los datos del siguiente modo plot(linspace(0,1,11),y,'b*')
hold on
plot(linspace(0,1,11),polyval(fit,linspace(0,1,11)))
El resultado es la figura 6.1

Figure 6.1: Ajuste por mínimos cuadrados de una serie de puntos
Resolver un problema de mínimos cuadrados es en el fondo resolver un problema mal condicionado. Si se plantean las ecuaciones del problema llegamos a un sistema lineal mal condicionado. El criterio de resolución suele ser el de minimizar el error cuadrático de la solución dada (por eso el nombre de mínimos cuadrados). Este es exactamente el mismo problema planteado en la pseudoinversa, operación que también se realiza por mínimos cuadrados.
Nada nos obliga a ajustar los datos mediante un polinomio; la condición de minimización del error cuadrático puede utilizarse con cualquier función aunque siempre será más sencillo hacerlo con polinomios. Una práctica muy habitual es utilizar un ajuste exponencial, táctica para la que Matlab no ha escrito ninguna función. El motivo es bien sencillo, si tenemos una nube de puntos que creemos sigue un modelo exponencial lo único que debemos hacer es calcular el logaritmo de la serie de datos para lo que recuperaremos automáticamente el ajuste polinómico.
6.1.1.1 ¿Qué calcula el ajuste polinómico por mínimos cuadrados?
Cuando hemos hablado de mínimos cuadrados en realidad nos referimos a los Mínimos Cuadrados Lineales Generales. El problema planteado es el de ajustar una serie de datos cualquiera a una función de la forma:
y(x)=a1+a2x+a3x2+…+aNxN-1
Este no es más que el tipo más sencillo de ajuste. Podríamos utilizar también series de funciones armónicas o funciones arbitrarias de la forma
en las que X representa una base de funciones. El objetivo es minimizar el funcional error representado por la expresión
| χ2= |
|
鼂R> 꼂R> 꼂R> 뼯FONT> |
|
|
? ? ? ?NT> |
|
donde σ es la desviación estándar de la muestra i-ésima. Esta formulación permite generar una matriz para una aplicación lineal:
Matriz que posee muchas más filas que columnas. La buena noticia es que acabamos de generar un problema tan simple como una aplicación lineal (ax=b), la mala es que la matriz no es cuadrada y no servirán los métodos de resolución convencionales. Una de las posibles estrategias de resolución es utilizar una SVD para resolver el sistema.
6.2 Interpolación y aproximación de funciones
La aproximación local o global de funciones es una parte esencial del cálculo y del análisis. Esta sección supone que ya conocemos los desarrollos de Taylor o las transformadas de Fourier. Algunos de los métodos planteados no son más que la extensión discreta de métodos comunes.
Aunque no sea del todo riguroso dividiremos este tema en tres partes según las características del desarrollo de nuestros datos. En cálculo numérico siempre nos veremos obligados a describir una función de un modo discreto, ya sea por el valor que toman en algunos puntos o por los coeficientes de un desarrollo. Dichas descripciones son más o menos adecuadas según la finalidad que busquemos.
Hablaremos de la interpolación polinómica a trozos cuando sean dados los valores de la función en unos puntos o nodos fijos. Utilizaremos funciones sencillas para intentar hallar una función continua que pase por todos los puntos. Hablaremos de interpolación polinómica cuando queramos modelar una función conocida o desconocida mediante funciones más simples buscando el mínimo error posible. Se diferencia de la anterior en que nuestro dato es una función y no una serie de puntos, esto nos permitirá escoger los nodos para buscar el mínimo error de interpolación. Por último trataremos a parte los desarrollos en serie de funciones sea cual sea la base (funciones armónicas, polinomios)
En el primer caso buscaremos convertir unos datos de naturaleza discreta en una función continua, en el segundo y el tercer caso intentaremos describir una función de un modo discreto, ya sea por sus valores o por unos coeficientes, por ejemplo para hallarla cuando es el resultado de una ecuación en derivadas parciales.
Los métodos numéricos no son exclusivos de cada uno de los objetivos, los splines, por ejemplo, sirven tanto para hallar una curva contínua como para describir una función mediante unos coeficientes. El uso de estas distinciones es que Matlab se basa más en la utilidad del método numérico que en su naturaleza. Si queremos splines para interpolar utilizaremos una función distinta de la dedicada a hallar sus coeficientes.
6.2.1 Interpolación polinómica a trozos
La interpolación polinómica a trozos será para nosotros la técnica de hallar una curva que pase por unos puntos dados. Se dice que es a trozos porque se utilizan polinomios de bajo orden definidos a intervalos.
- interp1
- Usa los puntos para interpolar en una dimension. Soprota interpolación discreta, lineal, cúbica, hermitiana y por splines cúbicos.
Un ejemplo gráfico es el siguiente script: xf=0:0.05:10;yf=sin(2*pi*xf/5);
xp=0:10 ;yp=sin(2*pi*xp/5);
lin=interp1(xp,yp,xf);
spl=interp1(xp,yp,xf,'spline');
cub=interp1(xp,yp,xf,'cubic');
near=interp1(xp,yp,xf,'nearest');
title('Comparacion de las opciones de interp1')
plot(xf,yf,xf,lin,xf,spl,xf,cub,xf,near)
legend('real','lineal','splines','cubica','discreta')
Cuyo resultado es la figura 6.2.

Figure 6.2: Comparación de los métodos de interpolación
- interp2
- Interpolación polinómica bidimensional.
De entre todos los métodos de interpolación a trozos el más efectivo suele ser siempre el uso de splines. Un spline es una curva cúbica; de este modo reduciremos en gran medida el error en el caso que los puntos sean de una función suave, cosa que no sabemos. Cuando interpolamos por trozos es importante conocer algo de información sobre la función que representan, no sirve ninguna formula universal
Sería muy lógico en este punto plantearse la siguiente pregunta. ¿Por qué escoger una función definida a intervalos y no un único polinomio de grado igual al número de puntos menos 1? La respuesta aparece cuando intentamos realizar el cálculo. Los polinomios de alto orden tienen tendencia a oscilar en sus extremos, es el llamado fenómeno de Runge. Como patrón diremos que todo lo relacionado con polinomios de alto orden suele acarrear problemas. Trataremos con más profundidad la interpolación polinómica en un dominio finito en la sección 6.2.3.
6.2.1.1 Splines
La interpolación polinómica a trozos más común en Matlab es la interpolación mediante splines. Un spline es una curva definida mediante polinomios en intervalos finitos, la interpolación será entonces del orden que deseemos según el polinomio aunque lo más frecuente será utilizar splines cúbicos. Como en muchos otros casos los splines no sirven únicamente para convertir en continua una serie de datos discreta, función que cumple interp1. Pueden servir también para analizar más profundamente los datos estimando sus derivadas e incluso para resolver ecuaciones en derivadas parciales.
Para quien desee algo de rigor matemático definimos un spline como el conjunto de polinomios pj(x) que ajustan una determinada función f(x) en [a,b], donde cada polinomio tiene validez en un intervalo ξj ∈ [a,b]. Además de imponer que la curva pase por todos los puntos podemos imponer también que dichos polinomios se unan suavemente al final de sus intervalos de definición. Podemos representar el conjunto mediante sus coeficientes y los nodos ( pp form)
o como B-spline, generalización de las curvas de Bézier.
Matlab, en su spline toolkit, proporciona funiones para calcular y evaluar con toda libertad los coeficientes de una interpolación a trozos mediante splines, piecewise spline interpolation. Algunas de las funciones disponibles son las siguientes:
- spline
- Llamada con dos parámetros retorna la matriz de coeficientes de la interpolación. Llamada con tres parámetros retorna el valor de la interpolación en los valores dados.
- ppval
- Evalúa la interpolación mediante splines en los puntos dados dada su representación en coeficientes.
- csapi
- Define
- csape
- Define
Un posible uso de los splines es la interpolación de curvas en más de una dimensión. Es un método muy habitual en diseño y es la base de las herramientas de CAD con las NURBS (Non-Uniform Rational B-Splines). Mediante puntos de control y unas pocas transformaciones geométricas los splines definen casi todas las piezas producidas industrialmente.
6.2.1.2 Regeneración de funciones mediante datos discretos
Como ya hemos mencionado con anterioridad, uno de los problemas de tratar con series discretas de datos es precisamente su no continuidad. Hay varias maneras de pasar los datos a una función y hemos analizado que el ajuste polinómico no siempre es una buena opción. Una solución más o menos universal es la de crear una función mediante una interpolación y un function handle. Para ello sólo tenemos que definir las series de datos y pasarlos como argumento a la vez que definimos una variable independiente adicional. Por ejemplo, suponemos que tenemos la siguiente serie de datos: x=[1 2 3 4 5 6 7 8]
y=[2 4 3 5 4 6 5 7]
Y queremos generar una función que sea capaz de proporcionar los datos de forma continua. La solución es la creación del siguiente function handle: >> newfunc=@(z) interp1d([1 2 3 4 5 6 7 8],...
[2 4 3 5 4 6 5 7],z,'spline');
A todos los efectos esto es una nueva función que pasa por todas las parejas de puntos (x,y): >> newfunc(3.443)
ans = 3.9162
La técnica de reducir variables mediante function handles no es válida en Octave, donde los FH son estrictamente funciones independientes.
6.2.2 Transformadas rápidas de Fourier
Un desarrollo posible de una función periodica es el desarrollo de Fourier. En él las funciones ortogonales que nos sirven de base son funciones armónicas que llamaremos Φk. El desarrollo de la función u es entonces:
donde u son los coeficientes de Fourier que suelen ser números complejos. Si la función admite el desarrollo anterior también admitirá una serie truncada de fourier. Esta serie tiene la forma:
que constituye la serie truncada de Fourier de orden N. Los coeficientes se calculan exactamente del mismo modo que la anterior con la diferencia que truncamos el desarrollo a partir de un cierto número de onda. En el cálculo numérico nunca trabajaremos con desarrollos con una cantidad infinita de términos; esto es dominio del cálculo simbólico.
Si aproximamos la función inicial por un polinomio interpolante definido a partir de una serie de puntos llegamos a la serie discreta de Fourier. Se dice discreta porque tiene como datos los puntos en los que se evalúe la función. Exactamente del mismo modo podríamos estar trabajando con una serie de puntos experimentales. Esto generará un desarrollo discreto porque sólo con unos cuantos puntos no tenemos suficiente información como para llegar a una serie infinita.
Como estamos obteniendo toda la información posible normalmente se habla de transformada discreta de Fourier. El algoritmo para calcularla es la tansformada rápida de Fourier o Fast Fourier Transform (FFT).
- fft
- Aplica la transformada rápida de Fourier a una serie de puntos utilizando como funciones base Φx(x)=eikx.
- ifft
- Calcula la antitransformada rápida de Fourier.
Las fft's se usan muchísimo en los métodos espectrales, resolución de ecuaciones en derivadas parciales lineales y no lineales, análisis de datos y filtrado de señales. Muchos de los problemas de la física con condiciones de contorno periodicas son mucho más sencillos en su formulación espectral.1
Disponemos también de otras funciones que encapulan tareas comunes en las que están involucradas transformadas de fourier
- fftshift
- Mueve el primer número de onda (frecuencia cero) al centro del espectro.
- fftfilt
- Aplica directamente un filtro espectral.
- fft2
- Transformada rápida de Fourier bidimensional. También disponemos de una antitransformada para esta función
- fftn
- Transformada rápida de Fourier n-dimensional. Existe ifftn.
Las librerías de transformadas rápidas de Fourier suelen disponer de transformadas del seno y del coseno para cubrir los casos de condiciones de contorno no periodicas. Los drivers para estas bibliotecas no son tan completos en Matlab y tendremos que convertir la transformada exponencial en transformada de seno o de coseno manualmente.