Prac2 Lagrange - Apuntes octave polinomio de Langrage PDF

Preview

Summary

Práctica 2: Interpolación polinómica2 Comandos nuevos2.1 MatricesLos elementos de una matriz se introducen, por filas, entre corchetes. Las filas se separan unas de otras por un punto y coma (;). Los elementos de cada fila se introducen separados por comas o por espacios en blanco. Por ejemplo, para...

Description

Práctica 2: Interpolación polinómica 2.1

Comandos nuevos

2.1.1

Matrices

Los elementos de una matriz se introducen, por filas, entre corchetes. Las filas se separan unas de otras por un punto y coma (;). Los elementos de cada fila se introducen separados por comas o por espacios en blanco. Por ejemplo, para obtener la matriz   1 2 3 4 5 6 se debe teclear: M = [1, 2, 3; 4, 5, 6]

M(i, j)

sirve para obtener el elemento de una matriz que ocupa en la fila i la posición de la columna j. Es decir, a la orden M(2, 1)

el programa responde: ans

= 4

Los comandos sobre manejo de matrices que utilizaremos en esta práctica son: • zeros(m, n) es la matriz de dimensión m×n cuyas componentes son todas nulas. • M(i, :) extrae la fila i-ésima de la matriz M. • M(i, j:k) extrae los elementos de la fila i-ésima que ocupan las posiciones de las columnas desde la j-ésima a la k-ésima. Por ejemplo, se extraen los elementos M52 , M53 y M54 de una matriz M con la orden >> M(5, 2:4) Además, se pueden yuxtaponer matrices con las instrucciones siguientes: • [A

B] es la matriz m×(n+p)-dimensional si A es (m×n)-dimensional y B es (m×p)-dimensional.

• [A; B] es la matriz (m+p)×n-dimensional si A es (m× n)-dimensional y B es (p×n)-dimensional. Nota 2.1 Lo expuesto aquí para matrices sirve también para vectores, pues un vector es un caso particular de matriz, ya sea de una sola fila o de una sola columna.

1

Práctica 2

2

2.1.2

Comando diff actuando sobre vectores

En la práctica siguiente veremos cómo emplear el comando diff para obtener derivadas de expresiones simbólicas que definan funciones de una variable. Sin embargo, aquí lo emplearemos sobre vectores de componentes numéricas. En ese caso, diff sirve para calcular la diferencia entre cada dos componentes consecutivas del vector que tenga como argumento, generando como respuesta un nuevo vector con una componente menos que el de partida. Por ejemplo, la orden diff([1, 3, -1, 0, 2])

proporciona como respuesta ans = 2

2.1.3

-4

1

2

Algo de programación

Aunque la programación en Octave no es el objetivo principal de las prácticas de esta asignatura, sí es cierto que en esta práctica emplearemos un fichero conteniendo un algoritmo —que implementa un método numérico— que se escribe de forma más elegante y sencilla si se emplea cierta estructura de programación (un bucle). Por ello, introducimos, de forma breve, las estructuras de repetición en su forma más simple. Estructuras de repetición (bucles) for k = 1:n sentencias endfor

Nota 2.2 • Si el bucle se teclea de forma interactiva en la ventana de comandos (es decir, sin utilizar el editor) las órdenes se ejecutan sólo tras teclear el último endfor que completa el bucle. • El salto del incremento del contador puede ser distinto de 1 sustituyendo k = 1:n por k = 1:m:n, donde m denota el salto. • El contador no tiene porque empezar en 1, podría ser k = p:n (o bien k = p:m:n, como caso más general). • El incremento puede ser negativo, para eso en k = p:m:n n debe ser menor que p. • Dentro del bucle no debe cambiar el valor de k (no son aconsejables instrucciones del tipo k = k+1). • No se deben utilizar ni i ni j para denotar el paso del bucle —pues en Octave representan a la unidad imaginaria de los números complejos—. • Los bucles anidados se ejecutan desde los más internos hacia los más externos.

2.2 2.2.1

Interpolación polinómica Polinomio de Lagrange

Aunque la interpolación de Lagrange es la más sencilla desde el punto de vista teórico, sin embargo es la menos utilizada en la práctica, por eso Octave no tiene ningún comando que la implemente. Por esta razón, en esta sección, se incluye el código de un archivo de función que calcula la expresión del polinomio de Lagrange conocidos los nodos de interpolación. Intercalados con el código se incluyen comentarios que explican las órdenes introducidas y que reproducen paso a paso la generación de la tabla de diferencias divididas (como la haríamos a mano). Usaremos ese archivo para resolver algunos de los ejercicios de esta práctica.

2.2. INTERPOLACIÓN POLINÓMICA function pol_lagr(X, Y) % funcion que obtiene la expresion del polinomio de Lagrange dados los vectores de las componentes de % los nodos de interpolacion % % pol_lagr(X, Y) tiene dos argumentos: % X es el vector de las abscisas de los nodos % Y es el vector de las ordenadas de los nodos % % Octubre 2021 % n = length(X); % numero de nodos M = zeros(n, n+1); % matriz nula de tantas filas como nodos % el numero de columnas es el de la tabla de diferencias divididas incluyendo la columna de las abscisas de los nodos M(:, 1) = X; % primera columna de M. Esta formada por las abscisas de los nodos M(:, 2) = Y; % segunda columna de M. Esta formada por las ordenadas de los nodos % % Diferencias divididas de primer orden % dd = diff(Y)./diff(X); % equivale a la orden dd = diff(M(:, 2))./diff(M(:, 1)); M(:, 3) = [dd'; 0]; % tercera columna de M. Esta formada por las diferencias divididas de primer orden . Se completa en la ultima fila con un cero % % Restantes columnas de la tabla de diferencias divididas % for k = 4:n + 1 dd = diff(dd)./(X(k-1:n) - X(1:n-k+2)); %en el paso k se forman las diferencias divididas de orden k−2 M(:, k) = [dd'; zeros(k-2, 1)]; % columna k−esima de la matriz M endfor % disp(''); disp('Tabla completa de diferencias divididas: '); disp(''); disp(M); disp(''); % % Coeficientes de la expresion de Newton del polinomio de Lagrange. % coef = M(1, 2:n+1); % Primera fila de la matriz M (a partir de la segunda columna). La forman la primera componente de Y y las diferencias divididas que apareceran en la expresion del polinomio % disp(''); disp('Los coeficientes de la formula de Newton son (primera fila de la tabla desde Y(1)): ') disp(''); disp(coef); disp('') % % Finalmente usamos la libreria symbolic para obtener la expresion del polinomio de Lagrange pkg load symbolic; warning off; % syms x; % monomios = x - X(1:n-1); % monomios que aparecen en la expresion del polinomio de Lagrange aux_0 = zeros(1, n - 1); aux = sym('aux_0'); % aux es ahora un vector nulo simbolico for k =1: n-1 aux(k) = prod(monomios(1:k)); % prod calcula el producto de las componentes del % vector monomios desde la primera a la k−esima.

3

Práctica 2

4 endfor pol = expand(coef(1) + sum(coef(2:n).*aux(1:n-1))); % las lineas de la 51 a la 58 equivalen las lineas de la 60 a la 63 (comentadas) %pol = x.^0.*coef(n); %for k = n:−1:2 %pol = expand(pol.*(x − X(k − 1)) + coef(k − 1)); %endfor disp('el polinomio de Lagrange es:') disp('') disp(pol); % Formula de Newton del polinomio de Lagrange % % Representacion grafica del polinomio resaltando los nodos de interpolacion % ezplot(pol, sym([X(1) - 0.2, X(n) + 0.2])) % representamos el polinomio obtenido hold on; plot(X, Y, '*') % resaltamos con asteriscos los nodos de interpolacion hold off endfunction

2.2.2

Otros tipos de interpolación

Aunque el polinomio que interpola los datos de una tabla es único, su expresión (como hemos visto en clases de teoría) puede venir dada de formas diferentes. La forma de Lagrange es de interés fundamentalmente teórico, pues se usa para deducir fórmulas, por ejemplo, de aproximación de integrales —como verás en el último tema del programa de esta asignatura— y en el ámbito práctico sólo se usa para interpolar tablas con un número pequeño de nodos. La fórmula de Newton es la más apropiada para la computación y también se usa para construir fórmulas para resolver ecuaciones diferenciales (aunque esto excede el propósito de esta asignatura). El hecho de que el polinomio de Lagrange se use poco a efectos prácticos se debe, entre otras causas, a que cuando se usan muchos nodos de interpolación, el polinomio es de orden alto y, en los puntos próximos a los extremos del intervalo de interpolación, se producen grandes oscilaciones (efecto Runge). Como la oscilación es una debilidad inherente a la interpolación polinomial, en el caso de un número grande de nodos, existen métodos mejores. Por eso, para reducir ese “ruido”, en esta sección, introducimos comandos que construyen polinomios en cada subintervalo generado por cada dos nodos consecutivos, [xi , xi+1 ] . El resultado es una función que, interpolando todos los nodos de una tabla, será polinómica a trozos. Esto puede hacerse de distintas formas, dos de las cuales veremos a continuación (usaremos estos métodos de interpolación, aunque no obtendremos de forma explícita las funciones interpoladoras). Interpolación polinómica lineal a trozos En este tipo de interpolación, en cada subintervalo [xi , xi+1 ] se considera una función lineal. El resultado es, pues, una función lineal a trozos (poligonal). Para este caso Octave sí tiene implementada una función llamada interp1 cuya sintaxis es: interp1(X, Y, V, 'linear')

donde X e Y denotan, respectivamente, los vectores de las abscisas y ordenadas de los nodos de interpolación. El tercer argumento, V, denota el vector de abscisas de los puntos en los que queremos evaluar la función interpoladora. (El último argumento es opcional, pues éste es el método de interpolación por defecto de Octave). Interpolación polinómica a trozos, de orden 3, con segunda derivada continua: trazador (spline) cúbico En este caso, en cada subintervalo [xi , xi+1 ] se considera un polinomio de orden 3 tal que, en cada nodo, la unión entre polinomios es “suave” pues la función resultante tiene derivada segunda continua. Los fundamentos teóricos de este método de interpolación exceden los objetivos de esta asignatura; no obstante, lo emplearemos para visualizar, en un caso práctico, cómo la aproximación obtenida con el trazador mejora la proporcionada por el polinomio de Lagrange. Con las mismas notaciones de antes, la sintaxis del comando que tiene implementado Octave en este caso es: interp1(X, Y, V, 'spline')

2.3. EJERCICIOS

2.2.3

5

Aplicaciones de la interpolación

En el mundo real la interpolación tiene múltiples aplicaciones. Se emplea por ejemplo en el ajuste de audífonos, en cartografía para la creación de mapas de isolíneas; esto es, mapas que proporcionan mediciones obtenidas de estaciones meteorológicas que se usan para estimar patrones del clima. En medicina, por ejemplo, se puede predecir la respuesta de un paciente a una dosis de un medicamento (sin habérsela administrado) si se conoce su respuesta a diferentes dosis, previamente administradas y, por supuesto, en el procesamiento de imágenes por ordenador. Las aplicaciones de diseño gráfico realizado con ordenadores requieren que se generen curvas de forma rápida y que los cambios en una porción de la curva no afecten al resto de la misma, o afecten de forma mínima, tanto por razones computaciones como estéticas. Esto descarta el uso de polinomios de Lagrange, debido a la oscilación que se genera cuando el número de nodos es alto. Las funciones interpoladoras, que suelen emplear estas aplicaciones informáticas, son polinomios cúbicos a trozos. Así, si se desea modificar una parte de la curva, se puede dejar el resto sin apenas variación, basta asegurarnos de que la unión entre los polinomios adyacentes es “suave”. De esta forma, la curva se puede modificar trozo a trozo. Puesto que los cálculos son mínimos la curva se modifica rápidamente y el resultado se puede ver de forma casi instantánea. Además, como los datos para la modificación se obtienen de las coordenadas de los nodos, el usuario no necesita tener ningún conocimiento matemático sobre cómo se realizan dichos cálculos. Las curvas en el espacio se generan de forma parecida, aunque aquí el mayor problema es la pérdida de la tercera dimensión al proyectarla sobre un plano (la pantalla del ordenador). Existen para ello diversas técnicas de proyección. Pero, dado que en el programa de la asignatura no se tratan funciones de dos variables, no abordaremos este tema.

2.3

Ejercicios

1. Calcula las diferencias divididas de primer orden de los nodos de la tabla: xi yi

1 -1

3 6

4 -2

7.5 0

8 15

2. Aplica el código de la sección 2.2.1 a los nodos de la tabla del ejercicio anterior y comprueba que el polinomio de Lagrange obtenido es correcto. Para esto, añade al código del archivo pol_lagr.m la línea, o líneas necesarias, que realicen la comprobación. 3. Se han tomado los tiempos de ejecución de un programa informático sobre n datos de entrada, para diversos valores de n, obteniéndose la tabla: n t (segs)

10 10

20 50

30 150

40 300

Calcula el tiempo de ejecución que emplea el programa cuando usa n = 25 datos de entrada. Hazlo usando la función pol_lagr y también a mano.

Práctica 2

6 4. Para dibujar el contorno superior del pato de la figura se utilizarán los datos de la tabla: xi yi

0.9 1.3

1.3 1.5

xi yi

12 0.6

12.6 0.5

1.9 1.85

2.1 2.1

13 0.4

13.3 0.25

2.6 2.6

3 2.7

3.9 2.4

4.4 2.15

4.7 2.05

5 2.1

6 2.25

7 2.3

8 2.25

9.2 1.95

10.5 1.4

11.3 0.9

11.6 0.7

que se han obtenido usando el sistema de coordenadas superpuesto a la figura del pato. Nótese que se usan más puntos del contorno en las zonas en las que la curva cambia rápidamente.

Figure 2.1: Gráfico obtenido de Métodos Numéricos, J.D.F and R.B, Thomson, pág 103 (a) Esboza el contorno superior del pato usando los datos de la tabla. Para ello, representa los puntos usando un asterisco para cada uno de ellos. (b) Calcula y representa gráficamente el polinomio de Lagrange que interpola los datos de la tabla. La representación hazla en una ventana gráfica diferente a la del apartado anterior. Observarás que al ser el polinomio de grado 20 se produce el efecto Runge (grandes oscilaciones que provocan un perfil muy distinto al que se esperaría del lomo de un pato). Se mejorará la aproximación en los dos apartados siguientes del ejercicio. (c) Representa gráficamente la función de interpolación lineal a trozos, para los datos de la tabla. Hazlo en una ventana gráfica diferente a las dos anteriores. (d) Calcula, usando interp1, el valor de la función que has dibujado, en los puntos x = 1.2, x = 2.9 y x = 11.25. (e) Representa gráficamente el trazador (spline) cúbico, para los datos de la tabla. Hazlo en una ventana gráfica diferente a las anteriores. Pista: Usa plot, el operador dos puntos (:) e interp1 con la opción spline. Incluye un título en cada gráfica que la identifique....