""
6.2.3 Aproximación de funciones
Esta técnica, aunque coceptualmente muy parecida a la interpolación polinómica a trozos, suele tener usos completamente distintos. La aproximación de una función o de una serie de puntos por un único polinomio en un dominio dado se acerca más al desarrollo en serie de Fourier que a la interpolación polinómica a trozos. Los polinomios de Lagrange, de Legendre o de Chebyshev fueron creados para desarrollar funciones continuas en dominios finitos mediante una base polinómica. Tienen un uso esencial en la evaluación de funciones complejas en dominios finitos y en los métodos espectrales de resolución de ecuaciones en derivadas parciales.
La interpolación polinómica no suele utilizarse para ajustar una serie de puntos dados por culpa del fenómeno de Runge. Supongamos que tenemos la siguiente serie de puntos y en función otra serie equiespaciada de puntos x. x =
Columns 1 through 8:
-1.00000 -0.87500 -0.75000 -0.62500 -0.50000 -0.37500 -0.25000 -0.12500
Columns 9 through 16:
0.00000 0.12500 0.25000 0.37500 0.50000 0.62500 0.75000 0.87500
Column 17:
1.00000
y =
Columns 1 through 8:
0.058824 0.075472 0.100000 0.137931 0.200000 0.307692 0.500000 0.800000
Columns 9 through 16:
1.000000 0.800000 0.500000 0.307692 0.200000 0.137931 0.100000 0.075472
Column 17:
0.058824
Supongamos ahora que queremos ajustarla por un polinomio de grado N-1 suponiendo N el número de puntos de la serie. El polinomio será entonces:
problema que se cerrará con la condición de que p(xi)=yi. Al final se llega a la misma ecuación de siempre:
Ac=y
donde A es la matriz de Vandermonde generada con los puntos xi, c el vector de incógnitas de los coeficientes y y el vector de puntos de la serie yi. El vector de coeficientes se genera con este código: >> c=vander(x) \y'
c =
6.4739e+03
3.6848e-11
-2.1040e+04
-1.0123e-10
2.7271e+04
1.0430e-10
-1.8231e+04
-5.0934e-11
6.8268e+03
1.2373e-11
-1.4741e+03
-1.4184e-12
1.8795e+02
6.5051e-14
-1.5402e+01
-8.4500e-16
1.0000e+00
Si representamos este polinomio de orden 16 con los puntos que tiene interpolar junto con la función solución y(x)=1/1+16x2(figura 6.3):
Figure 6.3: Demostración del fenómeno de Runge
Como vemos, la función interpolante cumple los requisitos impuestos, pasa por todos los puntos; pero no hace nada más bien. ¿Cuál es entonces el problema? ¿Qué hay que hacer para solucionarlo? La interpolación polinómica es un problema global, es decir, intenta aproximar toda una función con una cantidad limitada de datos (los puntos de los que tenemos valores). A la función interpolante los árboles no le dejan ver el bosque, si quiere ajustar con unos puntos dados (nodos) es incapaz de capturar la función de la que provienen.
La solución es comprender las implicaciones globales del problema. Ya no estamos intentando hacer pasar una curva por una serie de puntos, estamos intentando resolver la aproximación de la misma curva. La solución la encontramos desde el punto de vista global... ¿Por qué los nodos deben ser equiespaciados? ¿Y si acercamos los nodos entre ellos (clustering) en las zonas donde los errores son más elvados? Sean cuales sean los tipos de polinomios que utilicemos así como el orden del que sean hay nodos más adecuados para minimizar el error del problema global, por ejemplo los nodos de Chebyshev de la forma:
xj=cos(jπ/N), j=0,1,...,N
puntos de Chabyshev-Lobatto o puntos extremos de Chebyshev. Esta formula es la proyección de los puntos equiespaciados en una circunferencia de radio unidad en el eje de abcisas. Si ahora en vez utilizar los nodos equiespaciados utilizamos los nodos óptimos de Chebyshev llegamos a que el polinomio de interpolación es sensiblemente mejor (figura 6.4):
Figure 6.4: Uso de los nodos óptimos de Chebyshev para reducir el error de interpolación
Vemos entonces que recurrir a una elección óptima de nodos permite utilizar polinomios como base de desarrollos de funciones con un error más que aceptable con todas las ventajas que implica trabajar con polinomios en nuestros cálculos.
- Importante:
- La interpolación polinómica es un problema global. No depende únicamente del número de puntos sino de su elección.
Una demostración de hasta dónde llega la importancia de elegir bien los nodos es pensar que el error que se comete con un polinomio de grado mayor no necesariamente se reduce con lo que estamos malgastando tiempo y potencia de cálculo sin ganar precisión. Hay que hacer las cosas con cuidado.
6.3 Resolución de ecuaciones no lineales y optimización.
Este apartado merece un pequeño comentario de entrada. Matlab posee una más que amplia selección de funciones orientadas a la optimización. Matlab contiene funciones para hacer casi todo lo imaginable como buscar ceros de funciones n-dimensionales y no lineales, mínimos de funciones condicionados, programación cuadrática... El gran problema es que todas estas herramientas están en un toolkit que no forma parte de la instalación básica del programa. En el caso que estemos delante de una versión comercial o de estudiante del Matlab es posible que no dispongamos de ninguna de las funciones anteriormente mencionadas.
No tiene mucho sentido entrar los detalles del Optimization Toolkit porque la ayuda en formato HTML que se incluye en el paquete es más que suficiente. Se podría decir que es una buena introducción a la optimización en general.
Este tema es difícil de encajar en el planteamiento general del libro. Siempre hemos buscado un planteamiento agnóstico respecto a los dos intérpretes muy difícil de mantener en este caso. Octave sólo contiene las funciones básicas de optimización, la colección es pobre si se compara con el Optimization Toolkit; pero en la mayoría de los casos es más que suficiente.
Veremos entonces sólo una introducción básica a la optimización, minimización de funciones, programación lineal y no lineal y búsqueda de raíces de sistemas de ecuaciones no lineales.
6.3.1 Resolución de ecuaciones no lineales. Root finding.
La resolución de sistemas de ecuaciones no lineales podría considerarse como una disciplina independiente de la optimización. Algunos de los métodos de resolución no tienen absolutamente nada que ver con los algoritmos de minimización de funciones. Sin embargo hemos decidido mantener un órden consistente con el de los paquetes de funciones de Matlab.
Root finding es el término inglés para la práctica de la resolución de una ecuación que no tiene una solución analítica o con una solución exacta demasiado costosa de encontrar. Fue uno de los primeros casos en los que se aplicó el cálculo numérico porque es donde se hace imprescindible.
Todos los métodos se basan en la aproximación mediante un desarrollo en serie de la función en un punto inicial. Lo que los diferencia es el tipo de desarrollo y cómo aproximan la función en el punto. El método más básico de todos es el método de Newton que se basa en la aproximación lineal de la función en cada punto para obtener la siguiente raíz para iterar. El método de la secante no es más que una aproximación grosera de la derivada mediante dos puntos. Luego nos encontramos con los métodos Regula-falsi, Ridders... o más sofisticados como el Van Wijngaarden-Dekker-Brent.
Todos los métodos se usan de la misma manera, se les da una función, un punto inicial y se cruzan los dedos. Lo único que debemos saber del Root finding en Matlab es que todo lo lleva la función...
- fzero
- Busca la solución más cercana al punto inicial dado de cualquier ecuación dada.
Lo único que debemos recordar es que como cualquier función a la que hay que dar otra función como argumento nos obligará a utilizar un function handle. Por ejemplo, queremos encontrar el cero de la siguiente ecuación:
lnx-sinx=0
Esta ecuación no tiene una solución analítica con lo que el uso del cálculo numérico se hace imprescindible. Suele ayudar representar gráficamente las dos funciones (figura 6.5) para saber si se cruzan en algún punto; si buscamos una solución inexistente es probable que el método nos devuelva soluciones complejas o espúreas.
Figure 6.5: Representación de las funciones a resolver.
La función que se resolverá será siempre de la forma f(x)=0, de modo que nuestro function handle será: >> eq=@(x) log(x)-sin(x);
>> sol=fzero(eq,1);
>> sol
sol =
2.2191
que es consistente con la solución intuida basándonos en la gráfica.
Una más que buena práctica es pedir una solución dentro de un intervalo puesto que en muchos casos las ecuaciones no lineales pueden tener varias e incluso infinitas soluciones. Sólo hay que cambiar el punto inicial de iteración por un intervalo de búsqueda: >> eq=@(x) log(x)-sin(x);
>> sol=fzero(eq,[1 4]);
>> sol
sol =
2.2191
6.3.1.1 Más incompatibilidades entre Matlab y Octave.
En la versión 2.1 de Octave el soporte para los function handles es bastante limitado. Muchas de las funciones, sobretodo las que no son más que llamadas a rutinas en fortran o C, deben obtener las funciones por su nombre y no mediante un function handle. La función fzero en Octave es una víctima bastante desafortunada de esta limitación. Si intentamos resolver la ecuación anterior por el primer método propuesto recibiremos un error, a priori incomprensible, mientras que con el segundo script llegaremos a la solución.
El motivo es el exótico comportamiento de fsolve en Octave. Si le pasamos sólo un parámetro, fzero llamará a la función fsolve de la que hablaremos en el siguiente apartado. fsolve no tiene soporte para function handles de modo que tendremos que pasarle la función por su nombre despues de haberla definido, por ejemplo: >> function y=eq(x)
... y=log(x)-sin(x);
... end
>> fzero('eq',1)
ans = 2.2191
Tampoco es una gran dificultad añadida ya que Octave soporta la definición de funciones mediante el intérprete a diferencia de Matlab. En cambio, si definimos un intervalo de búsqueda de soluciones sí le podremos pasar un function handle como argumento porque usa un método de resolución alternativo2: >> eq1=@(x) log(x)-sin(x);
>> fzero(eq1,[1,4])
ans = 2.2191
6.3.2 Búsqueda de soluciones de sistemas de ecuaciones no lineales.
Se define como solución de un sistema de ecuaciones como la x que cumple:
f(x)=0
siendo f la función vectorial representa al sistema.
El objetivo de la búsqueda de soluciones de sistemas de ecuaciones no lineales es siempre el mismo. Llegar a una solución válida del sistema con el mínimo coste computacional posible. La no linealidad de las ecuaciones puede provocar todo tipo de catástrofes como la no convergencia o que el error se dispare hacia el infinito. Siempre se parará automáticamente la iteración y se nos dará un mensaje de error. El coste computacional es una complicación menor. Depende de dos aspectos, la elección del punto inicial y el algoritmo de evaluación del gradiente de la función.
Los métodos más sencillos de resolución de ecuaciones no lineales se basan en la aproximación lineal de la función en un punto. Si hacemos el desarrollo de Taylor de primer orden en el punto inicial tenemos que:
f(x)=f(x0)+(x-x0)J(x0)
fórmula a la que le siguen términos de mayor órden. J(x0) es el gradiente de la función evaluado en el punto x0. Si suponemos que el resultado de la aproximación es la verdadera raíz de nuestra función:
f(x0)+(x-x0)J(xo)=0
obtendremos el siguiente punto de la iteración, en este caso x=x1:
x1=x0-J-1(x0)f(x0)
Iteración que se llevará a cabo hasta que se una norma de la solución sobrepase una cota de error. Este método es conocido como Newton-raphson. Esto nos recuerda que una de las operaciones con mayor coste computacional es precisamente la inversión de una matriz. Además la matriz es el gradiente de la función en un punto, gradiente que puede ser en algunos casos casi imposible de calcular. Nos encontramos ante el problema de evaluar de una manera óptima un gradiente para luego evaluarlo e invertirlo. Todo esto considerando que el punto inicial nos lleve más o menos directamente a la solución que buscamos.
Nada de esto debe hacernos perder la referencia de que lo realmente importante de los métodos de resolución es la facilidad con la que converjan a una solución más que la rapidez con que lo hagan. Al igual que en el caso unidimensional, Matlab apuesta por la simplicidad. La función que nos realizará la tarea es:
- fsolve
- Busca una solución de un sistema de ecuaciones dado. Esta función es diferente según se use en Matlab o en Octave; la versión de Matlab llamará al sistema de ecuaciones mediante su nombre o con un function handle mientras que la versión de Octave sólo podra llamarlo por nombre.
Las diferencias entre ambas sólo se manifiestan en el modo en el que llamaremos a la función, no en el modo en el que la definiremos porque es imposible expresar todo un sistema de ecuaciones sólo con un function handle. Será imprescindible definir una función propiamente dicha lo que en Matlab significará crear un archivo de función.
6.3.2.1 Algunas de las cosas que pueden salir mal
Resolver ecuaciones escalares no lineales o sistemas de ecuaciones no lineales no es numéricamente muy complejo, pero los problemas son muy a menudo mal condicionados. Infinidad de factores pueden llevarnos a soluciones erróneas o a problemas de convergencia. Algunos de ellos son muy simples como el siguiente: >> eq=@(x) log(x)-sin(x);
>> fzero(eq,[0 4])
Warning: Log of zero.
> In @(x) log(x)-sin(x)
In fzero at 218
??? Error using ==> fzero
Function values at interval endpoints must be finite and real.
Obviamente debido a que ln0=-∞.
Otros problemas debidos a la convergencia pueden ser un auténtico dolor de cabeza
6.3.3 Minimización de funciones.(+)
6.3.4 Minimización de funcionales.(+)
- 1
- La gran ventaja de las transformadas rápidas de Fourier es que son una operación especialmente rápida cuando tenemos series de 2n puntos. En estos casos se hacen del órden de NlogN operaciones para calcularla. Esto las hace muy útiles en la resolución de ecuaciones en derivadas parciales lineales (ecuación de ondas) o por métodos espectrales. Un caso especialmente importante es el de la ecuación de Poisson (∇2φ=f(x)).
Supongamos que tenemos las ecuaciones de Navier-Stokes en dos dimensiones. Se demuestra que se pueden reescribir en función de ω y de ψ, la vorticidad y la función de corriente de la siguiente manera:
∇2ψ=-ω
Imaginemos entonces que queremos resolver un problema de turbulencia 2D isótropa con condiciones de contorno periódicas. Esto nos obliga a resolver una ecuación de Poisson por cada paso temporal, operación bastante costosa porque requiere la inversión de una matriz. Podemos ahorrarnos gran cantidad de operaciones si hacemos la transformada de Fourier de la segunda ecuación:
∂2ψ(i,j)exp(ı(kx+ly)) |
|
∂ x2 |
|
+ |
∂2ψ(i,j)exp(ı(kx+ly)) |
|
∂ y2 |
|
=-ωexp(ı(kx+ly)) |
ψ(i,j)(k2+l2)=ω(i,j)
Que es un sistema de ecuaciones de resolución trivial. Acto seguido hacemos la antitransformada de los coeficientes ψ(i,j) y ya podemos pasar al paso de la ecuación parabólica.
- 2
- fsolve utiliza el método Powell mientras que fzero utiliza un método Brent.
""
