Laboratorio de Cálculo Numérico MA-33A
Sesión: Interpolación
Introducción En esta sesión de laboratorio se interpolarán funciones usando el polinomio de interpolación y la función spline.
En la actividad 1, se implementarán diversas fórmulas que permiten calcular el polinomio de interpolación de una función dada.
En la actividad 2, se usará la mejor de las implementaciones estudiadad en la actividad 1, para graficar polinomios de interpolación de diversas funciones y con diferentes mallas de modo de poder apreciar en el computador los errores de interpolación.
En la actividad 3 se implementará la función spline y se graficarán algunas de ellas para compararlas con el polinomio de interpolación.
Finalmente, en la actividad 4 se usará un programa interactivo para estudiar las estabilidades de estas técnicas de interpolación, cuando los puntos de interpolación se mueven un poco.
1.
Diversas Div ersas Impl Implement ementacion aciones es del del polinom polinomio io de de Interpol Interpolación ación
Recordemos que dada una función f : : I ⊆ R → R y una mall mallaa de n puntos { x1 x2 P de grado ≤ (n − 1) que cumple la condición de interpolación ,
P( xi ) = f ( xi )
,
∀i ∈ {1 2 ,
,... ,
n}
existe te un un único polin polinomio omio xn }, exis
,... ,
.
Atención: Note que los puntos de la malla se han numerado de 1 a n (ya que en SCILAB los índices no parten en 0). Por ese motivo el polinomio será de grado ≤ (n − 1).
Para calcular este único polinomio de interpolación, existen diversas fórmulas, a saber: la de Vandermonde , la de Lagrange y la de Newton. A continuación veremos como se implementan cada una de ellas y cuales son sus ventajas y desventajas. Las pruebas la haremos con la función sen ( x) en el intervalo [0 4π] usando una malla equiespaciada de n puntos, con diversos valores de n. Escogemos esta función, ya que tiene todas sus derivadas acotad aco tadas as por 1, y por esto, esto, se pue puede de dem demost ostra rarr (v (ver er curso) curso) que cua cuando ndo n → ∞ , el polinomio de interpolación converge a la función. ,
1
1.1.
Fórmula de Vandermonde
Esta fórmula muy elemental se obtiene escribiendo el polinomio P en la base canónica como P( x) = a1 + a2 x + a3 x2 + ··· + an xn−1
.
Para encontrar las constantes ai se imponen las condiciones de interpolación de los datos, es decir a1 + a2 x1 + a3 x21 + ··· + an xn1−1 .. .
= .. .
f ( x1 ) .. .
a1 + a2 xn + a3 x2n + ··· + an xnn−1
=
f ( xn )
,
,
Las ecuaciones anteriores forman un sistema de ecuaciones lineales donde las incógnitas son los coeficientes del polinomio. En forma matricial este sistema se escribe
1 .. . 1
x21 .. . x2n
x1 .. . xn
··· ..
.
···
xn1−1 .. . n−1 xn
= a1 .. . an
f ( x1 ) .. . f ( xn )
Programación:
Abra el editor de SCILAB y edite allí un archivo de nombre vandermonde.sce. En ese archivo agregue las instrucciones que siguen.
Comience por definir la malla de n puntos (partirémos con n = 2) y los valores de la función en dicha malla. Use las lineas siguientes: n=2; a=0; b=4* %pi; XMalla=linspace(a,b,n)’;//vector columna YMalla=sin(XMalla);
El valor de n lo podrá cambiar más tarde.
Defina ahora la matriz de Vandermonde, columna a columna, usando el siguiente for : MV=[]; for i=0:(n-1) ;
MV=[MV XMalla.^i] ;
end
Ahora agregue las instrucciones SCILAB que resuelven el sistema MV*CoeffA=YMalla, para así calcular los coeficientes del polinomio de interpolación. (Recuerde que Ax = b se resuelve en SCILAB con x=A\b;). CoeffA=.........;//Complételo usted
Ejecute su programa pinchando el botón Load into Scilab del editor Si todo funciona bien debiera obtener (salvo redondeo) MV =
1 1
0 12 57 ,
,
CoeffA =
0 −3 9e − 17 ,
Esto lo puede verificar escribiendo directamente en la ventana SCILAB ( no en el editor ) MV enter CoeffA enter
2
Pruebas numéricas El polinomio de interpolación ya está calculado. Ahora lo graficaremos junto con la función seno para ver si se parecen o no. En toda la prueba numérica trabajaremos en el archivo vandermonde.sce variando el valor de la variable n . Además agregaremos las instrucciones de dibujo siguientes:
Para hacer el gráfico definiremos una malla auxiliar t como el siguiente vector columna: t=linspace(a,b,1001)’;
Para graficar seno definimos sus valores exactos en la malla auxiliar t como yexacto=sin(t);
Ahora, para evaluar el polinomio de interpolación en la malla auxiliar t , recordamos que P (t ) = ∑ni=1 ai t i−1 . Esta sumatoria se puede evaluar mediante un for del modo siguiente: PolVan=0; for i=1:n; PolVan=PolVan+CoeffA(i) * t^(i-1); end
Finalmente debemos hacer el gráfico. La idea es ver las dos curvas simultáneamente . Para esto, usaremos una instrucción más elaborada para graficar, ella es la instrucción plot2d(Mx,My). Esta instrucción recibe dos matrices Mx y My del mismo tamaño e interpreta cada columna como una curva. Estas curvas se dibujan simultaneas (O sea 2 columnas⇒2 curvas, 3 columnas ⇒3 curvas, etc.). En el caso que todas las columnas de Mx sean iguales se puede reemplazar la matriz Mx por un solo vector columna Vx . Esta última versión es la que usaremos aquí. O sea las lineas para dibujar serán: xbasc();xselect();plot2d(t, [yexacto PolVan]); plot2d(XMalla,YMalla,style=-3,strf=’000’); Note que [yexacto PolVan] es la matriz que tiene dos columnas con los dos gráficos que se desean. El segundo plot2d dibuja los puntos de interpolación con cruces (gracias al parámetro opcional style=-3).
Ejecute su programa para valores crecientes de n. Siga la secuencia n=2,4,6,8, etc. hasta que observe que el polinomio de interpolación no respeta la teoría (algo así como que el polinomio deje de pasar por los puntos). Guarde el primero de los gráficos incoherentes con el nombre vandermonde_ ”valor de n” .GIF (donde “valor de n” será el primer valor de n de la secuencia, para el cual el polinomio es másl calculado).
Agregue todos los comentarios relativos a esta experiencia numérica en un Documento Word de nombre Informe.doc. Agregue en este informe el gráfico guardado como archivo .GIF del punto anterior.
1.2.
Fórmula de Lagrange
La fórmula de Lagrange entrega el polinomio de interpolación en forma explícita, mediante la ecuación n
P(t )
=
n
(t − x j )
∑ f ( xi ) ∏ ( xi − x j )
i=1
j =1 j =i
Está fórmula se programa directamente en forma indicial usando un doble for .
Copie el archivo vandermonde.sce a lagrange.sce (basta con hacer “Save as”). En este nuevo archivo, borre todas las lineas que calculan el polinomio de interpolación usando el método de Vandermonde. Conserve solo las lineas que definen la malla, los datos de interpolación, la malla auxiliar y que realizan los gráficos. 3
Programe la fórmula de Lagrange (donde antes calculaba PolVan=...) mediante un doble for . Claramente la fórmula anterior se puede programar así: PolLag=0; for i=1:n pitat=1; II=1:n; II(i)=[];//índices diferentes de i for j=II; pitat=pitat .* (t-XMalla(j)) / (XMalla(...)-XMalla(...)); end //complete lo
que falta PolLag=PolLag + YMalla(i)*pitat;//los f ( xi ) están guardados en YMalla(i) end
Observaciones: La instrucción II(i)=[]; elimina del arreglo II el valor igual a i .
Así se puede hacer la pitatoria para
j =1 n =i j ,...,
mediante el for j=II;. Para comple-
tar el denominador que falta mire la fórmula de Lagrange.
Modifique la instrucción plot2d donde graficaba el polinomio PolVan para ahora graficar el polinomio PolLag. Ejecute el programa y pruebe con valores de n grandes siguiendo la secuencia 10,20,30,...,etc..
Guarde el gráfico de esta interpolación en que comiencen a manifestarse los errores de redondeo en la evaluación. Comente su experiencia en el archivo informe.doc. Comente las diferencias entre esta fórmula y la de Vandermonde, desde el punto de vista teórico y práctico. Agregue el gráfico guardado.
1.3.
Fórmula de Newton
La fórmula de Newton es una elaboración de la fórmula de Lagrange que la transforma en una fórmula incremental. El resultado final es P(t ) = f [ x1 ] + f [ x1 x2 ]( x − x1 ) + ··· + f [ x1 ,
donde los coeficientes f [ x1 ], f [ x1 x2 ] siva según la tabla ,
x1 x2 x3 x4
,... ,
f [ x1
xn ]( x − x1 ) ··· ( x − xn−1 )
,.. . ,
xn ] se llaman diferencias divididas y se calculan en forma recur-
,.. . ,
f [ x1 ]
ց f [ x2 ] → ց f [ x3 ] → ց f [ x4 ] →
f [ x1 x2 ] ,
ց f [ x2 x3 ] → ց f [ x3 x4 ] → ,
f [ x1 x2 x3 ]
,
ց f [ x2 x3 x4 ] →
,
,
,
,
f [ x1 x2 x3 x4 ] ,
,
,
Note que de esta tabla, solo interesan las diferencias divididas de la diagonal. Para programar esta tabla, las diferencias divididas de la diagonal se calcularán en un vector columna DifDiv. Inicialmente este vector contendrá los valores de f ( xi ) y en forma recursiva se irán modificando sus valores, de modo de avanzar hacia la derecha, modificando solo los valores bajo la diagonal. La idea es que al final, este vector contenga las diferencias divididas de la diagonal de la tabla. El esquema de cálculo de estas diferencias divididas se puede visualizar del modo siguiente: 4
i = 2
∆1 ∆2
∆2
∆3
∆3
∆3
∆4
∆4
∆4
∆4
∆5
∆5
∆5
∆5
∆5
∆6
∆6
∆6
∆6
∆6
∆7
∆7
∆7
∆7
∆7
∆8
∆8
∆8
∆8
∆8
∆9
∆9
∆9
∆9
∆9
i = 5
Algoritmo de Reactualización Desde i = 2 hasta n : DifDiv j ←
DifDiv j − DifDiv j−1
∀ j ∈ {i
,... ,
x j − x j−(i−1) n}
Copie el archivo lagrange.sce a lagrange-newton.sce. Aquí no borre la fórmula de Lagrange para poder compararla con la fórmula de Newton.
Programe el cálculo de las diferencias divididas del modo siguiente. DifDiv=YMalla; for i=2:n; j=i:n; DifDiv(j)= (...) ./ (....); end
Nuevamente complete usted lo que falta mirando las fórmulas. n
Programe el cálculo de P (t ) usando la fórmula de Newton P (t ) =
i−1
∑ DDi ∏(t − x j ) i=1
j =1
mediante un doble for . PolNew=0;//este será el polinomio evaluado en t for i=1:n; pitat=1; for j=1:(i-1); pitat=pitat .* ..... ; end; PolNew=PolNew+DifDiv(i)*........; end
Complete usted lo que falta.
2.
Ahora modifique el plot2d para graficar simultáneamente las tres curvas: f ( x), Polinomio de la fórmula Lagrange y de la fórmula Newton. Ejecute el programa para valores grandes de n (donde comienzan a verse los errores) y agregue sus comentarios al s.
Interpolaciones varias de otras funciones
Para aprovechar el programa hecho, visualizaremos polinomios de interpolación obtenidos para otras funciones. Esto lo haremos en un nuevo archivo llamado interpolaciones.sce (guarde su programa anterior con este nuevo nombre). De todas las formas para calcular el polinomio de interpolación nos quedaremos con la fórmula de Newton. Comenzaremos por limpiar el programa de la actividad anterior.
2.1.
Transformación del código a una forma más operativa
Para que el programa de interpolación sea más operativo, definiremos dos funciones SCILAB que nos permitan fácilmente encontrar el polinomio de interpolación de una curva. La primera función calculará las diferencias divididas y la segunda evaluará el polinomio en una malla auxiliar t i . 5
Escriba un nuevo archivo llamado funciones.sci. En este archivo comenzaremos por definir la función difdiv(X,Y), que recibe la malla X y los datos Y y calcula la tabla de diferencias divididas según la estrategia ya estudiada en la sección anterior. Para definir esta función escriba en el archivo funciones.sci lo siguiente: function D=difdiv(X,Y) N=length(X); X=matrix(X,N,1);Y=matrix(Y,N,1);
Esta instrucción transforma el arreglo X en un vector columna con las mismas componentes (esta es una precaución necesaria para que la función sea más robusta). D=Y; for i=2:N j=i:N;
D(j)= (D(j)-D(j-1)) ./ (X(j)-X(j+1-i) );
end
Note que la variable auxiliar D es la que será devuelta por la función con las diferencias divididas.
Ahora definiremos la función EvalNewton que recibe la malla auxiliar t, la malla de interpolación X y las diferencias divididas ya calculadas D. Esta función entregará la evaluación del polinomio de interpolación en cada uno de los puntos definidos por la malla t . Para definir esta función, basta con agregar las siguientes lineas a continuación de las que definían a la función difdiv anterior (o sea la instrucción function termina con la definición de la función anterior): function Pt=EvalNewton(t,X,D) N=length(X);
Pt=0;
for i=1:N p=1;
for j=1:(i-1); p=p .* (t-X(j)); end;
Pt=Pt+D(i)*p; end
Si aun no lo ha hecho, guarde su archivo Lagrange-Newton.sce con el nombre interpolaciones.sce. Allí agregue al principio la instrucción getf(“funciones.sci”);
con esto se cargarán las definiciones de las funciones anteriores.
En la parte donde usted calculaba el polinomio de interpolación, borre todas las lineas innecesarias y cámbielas por solo 2 lineas que usan las dos funciones difdiv y EvalNewton. Es decir donde corresponda escriba DD=difdiv(XMalla,YMalla);
y para evaluar el polinomio en la malla auxiliar t escriba Pol=EvalNewton(t,XMalla,DD);
Si ejecuta su programa en este momento, debiera funcionar igual que al final de la actividad 1. Verifíquelo. Recuerde que ahora hay que graficar solo el polinomio llamado Pol junto a la función seno.
Para interpolar otras funciones, agregue al archivo funciones.sci las siguientes definiciones: function y=f1(x) y=sin(x);//esta es la misma función de antes function y=f2(x) y=sin(4*x);//esta otra función function y=f3(x)
6
y=1 ./(1+x.*x) function y=f4(x) y=sin(x)+sin(4*x)
Para interpolar alguna de estas funciones, defina la función ff=f1; (o f2 , o f3, etc). Cambie en el programa, en todas partes donde decía sin por ff (son 2 cambios en total).
Interpole algunas de estas funciones para diferentes valores de n. Agregue los mejores gráficos, con las explicaciones correspondientes en su informe.doc. Note que seguramente deberá cambiar los rangos de interpolación definidos por a y b. Por ejemplo, la función f3 conviene interpolarla entre a=-5 y b=5 con 10 ≤ n ≤ 20.
3.
Interpolación Spline Cúbica
Dada una partición { x1 x2 xn } (aquí los puntos deben estar ordenados) y un conjunto de datos { y1 y2 yn } 2 la función spline de interpolación de esos datos es una función s ( x) de clase C cuya restricción a cada subintervalo [ xi xi+1 ] es un polinomio de grado 3. ,
,.. . ,
,
,.. . ,
,
,
En el curso MA33A se demuestra que para calcular la Spline en cada intervalo, basta con determinar las derivadas y′i = s′ ( xi ) en los puntos de la malla resolviendo el siguiente sistema
2 h11 1 h1
1 h1
2
1 h2
+ h11
..
.
..
.
..
.
..
.
..
.
..
.
1 hn−1
1 hn−1
2h1
n−1
y′1 y′2 .. . y′n
=
3
m2 h2
m1 h1
+ mh11 .. .
mn−1 + mhnn−−22 hn−1 mn−1 hn−1
Una vez resuelto este sistema, los polinomios que definen la spline quedan definidos en cada intervalo [ xi xi+1 ] por la siguiente fórmula de Hermite ,
′
P(t ) = yi + yi (t − xi ) +
donde m i =
yi − yi−1 hi
y h i = xi − xi−1
m i − y′i hi
2
(t − xi ) +
y′i + y′i+1 − 2mi h2i
(t − xi )2 (t − xi+1 )
(1)
.
Agregaremos al archivo funciones.sci dos funciones análogas a las usadas en interpolación polinomial. La primera, de nombre Yp=Spline_Yp(X,Y) resolverá el sistema que permite calcular los valores de las pendientes y′i (para esto recibirá los puntos de la malla X y los valores de la función Y. La segunda de nombre Spline_eval(X,Y,Yp,t) evaluará la spline en la malla t de puntos auxiliares.
Agregue a su archivo funciones.sci la definición de las funciones Spline_Yp y Spline_eval. La programación es la que se entrega a continuación: function Yp=Spline_Yp(X,Y) n=length(X); X=matrix(X,n,1);//se redimensiona X como vector columna. Y=matrix(Y,n,1);//idem para Y [XX,i]=sort(-X);X=X(i);Y=Y(i);
Esta linea reordena los datos de menor a mayor. (OBS: la función sort ordena de mayor a menor, por eso se usó el signo menos. Ella envía un conjunto de índices que sirve para ordenar los otros vectores consecuentemente) 7
I=1:(n-1); h=X(I+1)-X(I); rh= 1 ./ h; A=diag( [rh;0]+[0;rh] )*2 + diag(rh,1)+diag(rh,-1);
Estas lineas definen la matriz del sistema que permite calcular las derivadas de la función Spline. m=(Y(I+1)-Y(I)) ./ h.^2; B=([m;0]+[0;m])*3;
Estas lineas definen el lado derecho del sistema que permite calcular las derivadas de la función Spline. Yp=A\B;//solución del sistema function S=Spline_eval(X,Y,Yp,t) n=length(X); X=matrix(X,n,1); Y=matrix(Y,n,1); [XX,i]=sort(-X);X=X(i);Y=Y(i); //reordenamiento for i=1:(n-1); a=X(i);b=X(i+1);h=b-a;m=(Y(i+1)-Y(i) ) / h; J=find( t>= a & t<=b );
Esta instrucción es para saber cuales son los puntos t j que están en el intervalo actual S(J)=Y(i)+Yp(i)*(t(J)-a)+(m-Yp(i))/h*(t(J)-a).^2 ... +(Yp(i)+Yp(i+1)-2*m)/h/h*(t(J)-a).^2 .* (t(J)-b);
Estas lineas son la evaluación de la Spline según la fórmula (1) end
Copie su programa interpolaciones.sce en interpolaciones2.sce. Modifíquelo para llamar a las funciones Spline_Yp y Spline_eval apropiadamente. Grafique simultáneamente el polinomio de interpolación y la función spline. Indicaciones: Yp=Spline_Yp(XMalla,...); S=Spline_eval(XMalla,YMalla,Yp,t);
4.
Ejecute este nuevo programa con las diversas funciones usadas en la actividad anterior. Guarde algunos gráficos en Spline_ ”valor de n”.gif e incorpore todas las observaciones que pueda en su informe.doc. En particular vea hasta que valor de n se puede trazar la spline (cuando n sea muy grande, comente la zona donde se calcula el polinomio de Lagrange). Estudie las oscilaciones que presentan los dos métodos cuando se usan para la función f3 .
Estabilidad de los dos métodos de interpolación
En la página del Laboratorio hay 2 archivos que interactuán con sus programas, generando una versión de interpolación interactiva. Con ellos podrá modificar los parámetros de interpolación, usando el mouse y el teclado. Así no solo podrá ver la función de interpolación de diversas funciones sino que además podrá mover los puntos de interpolación y percibir cual es la sensibilidad de los resultados obtenidos cuando los puntos de interpolación se mueven ligeramente. Si al mover los puntos, la curva de interpolación cambia poco se dice que es estable. Por el contrario, si cambia mucho se dice que es inestable.
Guarde los archivos funciones_pagina.sci e interpolaciones_pagina.sce a su carpeta personal. 8
Ejecute el archivo interpolaciones_pagina.sce, use las instrucciones siguientes par modificar los datos de interpolación. Vea que pasa cuando n es pequeño para el polinomio de interpolación y la función spline. Haga lo mismo cuando n es grande. Agregue todos los comentarios pertinentes en su informe.doc
Las teclas preprogramadas en los archivos nuevos son: [+]
Aumenta el valor de n . Los puntos de interpolación serán equiespaciados.
[-]
Disminuye el valor de n . Los puntos de interpolación serán equiespaciados.
[1]...[9] Cambian la función a interpolar (Hay 9 funciones predefinidas) [flechas] Mueven un poco al punto marcado en verde (solo funciona en linux). [a][s][d][w] Funcionan como las flechas para ser usado en windows (como lo hacen algunos juegos). [mouse]
Selecciona un punto verde. Si se mueve el mouse, el punto verde se moverá hasta soltar el mouse.
[r]
Cambia un poco las coordenadas de todos los datos de interpolación en forma aleatoria (equivale a un ruido numérico).
[c]
Dibuja curvas en lugar de funciones. Sobre todo las funciones 8 y 9 han sido programadas como x = f (t ) y y = g(t ).
[f]
Vuelve a dibujar funciones en lugar de curvas.
[l]
Interpola usando el polinomio de Lagrange.
[p]
Interpola usando la función S(p)line natural (vista en clases).
[n]
Interpola usando la función Spline not-a-knot (es una versión avanzada de la spline que no termina con segunda derivada nula en los extremos).
9
