Este capítulo trata de todas las funciones dentro de la biblioteca relacionadas con el cálculo infinitesimal y el análisis matemático. Matlab no es un lenguaje de cálculo simbólico, no calcula ni derivadas ni primitivas. Todas las rutinas de este capítulo están enfocadas al cálculo numérico.
5.1 Funciones elementales.
Ya hemos comentado que Matlab posee una importante biblioteca de funciones, es el momento de conocerla con un poco más de detalle
5.2 Polinomios
En Matlab los polinomios se almacenan como vectores de coeficientes. Cada componente del vector es un coeficiente ordenado del término de mayor a menor orden. La manera más habitual de iniciar un polinomio es utilizar la siguiente función:
- poly
- Inicia un polinomio. Si el argumento es una matriz cuadrada devuelve el polinomio característico, si es un vector devuelve un polinomio cuyas raíces son los elementos del vector.
También nos serán muy útiles las funciónes:
- roots
- Devuelve las raíces del polinomio cuyos coeficientes son los elementos del vector argumento.
- polyval
- Calcula el valor del polinomio en el punto dado
- polyout
- (Octave) Escribe un polinomio en un formato más leíble:
Por ejemplo, si queremos escribir el polinomio s5-8s4+3s2+16: >> polinomio = [1,-8,0,+3,0,16]
>> polyval(polinomio,44)
ans= 134937280
Los polinomios se pueden manipular analíticamente con las siguientes funciones:
- polyder
- Deriva un polinomio.
- polyint
- Integra analíticamente un polinomio.
- conv
- Multiplica dos polinomios.
- deconv
- Dados dos polinomios a e y halla b y r que cumplen:
y=(a*b)+r
- residue
- Descompone el cociente de dos polinomios en fracciones con un único polo. Por ejemplo, dados los polinomios p=s2+s+1 y q=s3-5s2+8s-4 la descomposición sería:
5.2.1 Evaluación de funciones.
Los polinomios son una parte importantísima del cálculo numérico porque son su nexo fundamental con el cálculo analítico. Uno de los elementos fundamentales en el cálculo numérico es la aproximación de una función analítica por otra con una expresión numérica lo más equivalente posible. Todas las funciones analíticas no existen realmente dentro de un programa de cálculo numérico por limitaciones del propio ordenador; la aritmética exacta es algo prohibitivo. La evaluación de funciones se realiza por series de potencias almacenadas por sus coeficientes; este es el funcionamiento de la evaluación de cualquier función analítica en un programa de cálculo numérico. Las series de potencias acotan el error con el número de términos y el orden de la serie; a mayor precisión más términos a evaluar. Por ejemplo, si Matlab quiere obtener el valor de la función seno en un punto evaluará la serie:
O si quiere evaluar la función de Bessel de primera especie:
| Jn(x)= |
漂R> 缂R> 缂R> 輯FONT> |
|
?> ?> ?> ?NT> |
|
|
|
|
Si este es el funcionamiento interno de Matlab será muy beneficioso que tengamos esto en cuenta cada vez que tengamos que evaluar funciones especialmente pesadas. No necesitamos aritmética exacta, evaluar una precisión mayor a la de la máquina es malgastar recursos.
Supongamos que tenemos una función con una gran cantidad de funciones trigonométricas. Es probable que alguna vez en nuestra vida tengamos funciones de dos y tres líneas de código. Evaluarla entera es una pérdida de recursos. Imaginemos además que estemos integrando una ecuación diferencial cuyo término es esta misma función. Evaluando todos los términos no estamos ganando precisión; estamos sumando y restando polinomios de un orden mayor al que queremos con valores que se cancelan continuamente. Utilizando la interpolación polinómica podemos acelerar cualquier código hasta en un órden de magnitud.
5.3 Derivadas
Todas las funciones de este apartado se refieren a la derivada numérica. No son derivadas de una función analítica sino el resultado de aplicar las ecuaciones de diferencias de orden enésimo a una serie discreta de puntos.
- diff
- Aplica la fórmula de diferencias enésimas a una serie de puntos
- gradient
- Calcula el gradiente numérico de una matriz. Esta función es muy interesante por su relación íntima con la función quiver. Si queremos expresar un determinado campo vectorial a partir de su función potencial, siendo el potencial la función φ(x,y) entonces el campo es en cada punto (∂xφ,∂yφ), es decir, ∇φ. Para pintar este campo vectorial sólo tenemos que hacer quiver(gradient(phi)) .
5.4 Integrales
La integración numérica es un vasto campo del cálculo numérico y así se refleja en Matlab. Podemos hacer tanto integrales en una dimensión como en dos. Las rutinas más relevantes son:
- quad
- Integración numérica de una función analítica introducida a través de un archivo de función o por una sentencia inline, es en este tipo de funciones donde el inlining tiene más utilidad. Por ejemplo, para integrar ∫010exp(x4) dx haremos quad(inline('exp(-x4)'),0,10). La función nos permite introducir la tolerancia requerida para el resultado, cuando pasemos este argumento es importante que no sea mayor que la precisión numérica con la que estemos trabajando o haremos que el programa produzca un error. También podemos pasarle otros argumentos para aumentar la precisión como marcar ciertos puntos donde la función puede presentar singularidades.
- quadl
- Su uso es idéntico que la función anterior pero contiene un algoritmo de cuadratura de mayor precisión. Obviamente a mayor precisión mayor tiempo de cálculo.
- dblquad
- quad2dc y quad2dg en Octave. Calcula la integral doble de una función de dos variables. Obviamente la frontera de integración serán cuatro constantes.
- triplequad
- No existe aún en Octave. Calcula la integral de volúmen de una función de tres variables.
- trapz
- Aplica la regla del trapecio para calcular la integral numérica de una serie discreta de puntos.
5.5 Ecuaciones diferenciales ordinarias.
Se podría decir que hay un método de integración para cada problema. Existen multitud de esquemas de integración y escoger uno de ellos no es fácil. Uno de los objetivos de las funciones de Matlab es intentar ser lo más universales posible. La función perfecta sería la que sin configuración alguna integrara cualquier EDO con una precisión aceptable; desgraciadamente no existe y probablemente no exista nunca.
Una consideración esencial antes intentar resolver una EDO es saber si el paso temporal viene condicionado por la precisión numérica o por criterios de estabilidad. El segundo caso es bastante más preocupante que el primero. Cuando la función introduce gradientes elevados que acortan significativamente el paso de la iteración se dice que el problema es stiff. En estos casos será apropiado utilizar un esquema de integración implícito que a costa de mayor coste computacional por paso permite llegar al resultado con muchas menos iteraciones. Saber si un problema es stiff o no es tan laborioso o más que resolverlo. Para saber si utilizar un esquema explícito o implícito la estrategia será la siguiente:
- Represenataremos la función a integrar para comprobar si tiene discontinuidades o gradientes muy elevados
- Si no tiene ninguna singularidad utilizaremos un esquema explícito (no stiff) calculando el tiempo de proceso.
- Si el tiempo nos parece demasiado elevado cambiaremos a un esquema implícito
- Si el tiempo sigue siendo demasiado elevado y la programación es correcta puede ser que la función sea demasiado costosa de evaluar. Para solucionarlo escribiremos la función en otro lenguaje de programación como C, C++ o Fortran y la convertiremos en un archivo que Matlab entienda.
La precisión numérica es también muy importante pero el mismo esquema de integración es capaz de controlarla. Daremos como válido que a mayor precisión numérica mayor dificultad para obtener una solución estable.
Las diferencias entre Matlab y Octave en este caso son más que significativas. Octave dispone de las funciones de integración de Matlab pero existen sólo por motivos de compatibilidad, son mucho más lentas y no se pueden configurar. Sin embargo Octave integra la suite ODEPACK, uno de los paquetes de integración de EDOs más robusto y efectivo.
5.5.1 Octave
Octave cuenta con las funciones de integración de Matlab pero están programadas como archivos .m y no como parte de una biblioteca binaria. Es muy recomendable que utilcemos la función lsode (Livermore Solver for Ordinary Differential Equations), parte de ODEPACK. Utiliza un método multipaso Adams Para problemas no stiff y una BDF (Backward Differentiation Formula) para los problemas stiff.1
- lsode
- (Octave) Es el driver general de integración de ecuaciones diferenciales ordinarias.
[X, ISTATE, MSG] = lsode (FCN, X_0, T, T_CRIT)
Resuelve el problema de Cauchy de la forma:
dadas las condiciones iniciales.
- FCN
- Matriz de dos elementos tipo cadena de caracteres. El primero es el nombre de la función que queremos integrar y el segundo la función que calcula el jacobiano.
La función a integrar debe definirse necesariamente mediante una función, ya sea dentro del script o en un archivo a parte. La cabecera tendrá sólo dos variables y su forma será XDOT=FUNC(X,T)
donde tanto X como XDOT pueden ser escalares o vectores. El jacobiano también tendrá la forma JAC=FJAC(X,T)
donde JAC será la matriz que implemente las funciones derivadas parciales según el orden establecido por la variable X. Se usa para acelerar la integración en problemas stiff.
- X_0
- Vector de condiciones iniciales.
- T
- Vector que determina los instantes temporales en los que se dará una solución.
- T_CRIT
- Puntos en los que la rutina no debe integrar por la existencia de singularidades.
- X
- Vector de soluciones en los instantes T
Si la integración no ha sido satisfactoria ISTATE devolverá el valor 2 y MSG nos dará más información.
- lsode_options
- Función de configuración de lsode
lsode_options index{lsode options} (OPT, VAL)
Llamada con un único argumento (OPT) nos dirá el valor de la opción de configuración que hayamos solicitado. Pasando dos argumentos podremos cambiarlos. Las opciones disponibles son los de la tabla siguiente:
Table 5.1: Valores de configuración de lsode_options
Pronto veremos un ejemplo de cómo usar esta función.
5.5.2 Matlab
El planteamiento de las rutinas para ODEs de Matlab es completamente distinto. En vez de tener una función que sirve más o menos para todo nos permite escojer el método de integración. Disponemos de los presentes en la tabla 5.2:
Table 5.2: Funciones de integración de EDOs de Matlab
La llamada básica es la misma para todas ellas y tiene la forma: [T,Y] = solver(odefun,tspan,y0)
[T,Y] = solver(odefun,tspan,y0,options)
[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)
odefun es la llamada convencional a una función, preferiblemente un function handle. A diferencia de lsode la función debe tener la forma:
y su resultado debe ser un vector columna. tspan es un vector de dos componentes con el tiempo inicial y final de integración; finalmente y0 es el vector con la condición inicial. Los parámetros de salida son los usuales, T e Y son respectivamente los vectores tiempo y solución del problema. En breve propondremos un ejemplo
También es necesario que conozcamos la función de configuración odeset. La cantidad de opciones es bastante mayor a las disponibles en lsode de modo que lo más sensato será consultar la ayuda.
5.5.3 Solución de la ecuación de Van der Pol
La ecuación diferencial de Van der Pol es:
| mx′′-α x′+β |
( |
x′ |
) |
3+kx=0 |
No es más que el movimiento de una masa oscilante con amortiguamiento lineal y no lineal a la vez. Adimensionalizando el problema llegamos a la ecuación.
x′′+x+?(x′2-1)x′=0
Esta ecuación tiene la particularidad de que a medida que el parámetro ? aumenta el problema se vuelve stiff. Con un ? de orden unidad el sistema puede considerarse no stiff mientras que si aumenta hasta ser del orden del millar la ecuación introduce gradientes muy acusados. Este comportamiento tan particular la hace perfecta para ensayar los distintos métodos de integración. Lo primero es descomponer la ecuación en un problema que Matlab pueda entender:
con condiciones de contorno y1(0)=2 e y:2(0)=0. Resolveremos el problema para ?=1 y ?=1000 porque Matlab ya cuenta con estas funciones (vdp1 y vdp1000). Para agilizar los cálculos escribiremos las funciones para Octave en C++: #include <octave/oct.h>
DEFUN_DLD (vdp1,args, ,
"Ecuacion de Van der Pol para mu=1 ")
{
ColumnVector xdot (2);
ColumnVector x (args(0).vector_value());
float mu=1;
xdot(0) = x(1);
xdot(1) = mu*x(1)*(1-x(0)*x(0))-x(0);
return octave_value (xdot);
}
Y la función para ?=1000: #include <octave/oct.h>
DEFUN_DLD (vdp1000,args, ,
"Ecuacion de Van der Pol para mu=1000 ")
{
ColumnVector xdot (2);
ColumnVector x (args(0).vector_value());
float mu=1000;
xdot(0) = x(1);
xdot(1) = mu*x(1)*(1-x(0)*x(0))-x(0);
return octave_value (xdot);
}
