Las bondades de máxima y un ejemplo para el cálculo de diferencias divididas

Aunque máxima no es mi herramienta favorita, hay que reconocer que en ocasiones su inmediatez lo hace muy útil. En este caso os traigo un rápido desarrollo para calcular una tabla de diferencias divididas, en parte gracias a una de las mejores utilidades que ofrece máxima: makelist.

¿Qué son las diferencias divididas?

Cuando queremos calcular el polinomio que interpola una serie de puntos, una de las herramientas más importantes es lo que se da en llamar una tabla de diferencias divididas. Aquí podéis encontrar más información sobre el método:

Método de las diferencias divididas de Newton

¿Cómo calcularlas en máxima?

El posible algoritmo para calcularla tendría dos opciones, una que sería recursiva y otra opción sería calcularlas de forma iterativa. Yo no aconsejaría usar nunca recursividad en Máxima, ya cuesta a veces depurar un código normal en este lenguaje, no quiero imaginarme depurar uno recursivo. Así que tomo la opción iterativa para calcular las diferencias.

Para empezar definimos el formato de entrada para la función, que será una matriz nx2, en el que los valores (x,y) para cada punto serán [(i,1), (i,2)]

Por ejemplo:

(%i1) m:matrix([3.2,22],[2.7,17.8],[1,14.2],[4.8,38.3],[5.6,51.7]);

Para definir la función que calculará las diferencias divididas tenemos el siguiente código (al final de la entrada encontraréis un enlace para ver el código con comentarios):

difdiv(m):= block(
    [j,l,r],
    r: copymatrix(m),
    for j:1 thru (length(r)-1) do  block(
        l:makelist((r[i+1][j+1]-r[i][j+1])/(r[i+j][1]-r[i][1]),i,1,length(r)-j),
        l:append(l,makelist(0,k,1,j)),
    r:addcol(r,l)
    ),
    r
)$

En cada iteración del bucle se calcularán las diferencias divididas de cada nivel y luego se añadirán a la matriz original. De esta forma la salida del algoritmo será la matriz original y las diferencias divididas de cada nivel.

Ahora veamos los comandos en mayor profundidad:

difdiv(m):= block(

Definición de la función

[j,l,r],

Definición de variables locales. j será el contador principal del método, l será una variable en la que se guardarán temporalmente listas y r será una copia de la matriz de entrada sobre la que trabajaremos.

r: copymatrix(m),

Se realiza una copia de la matriz m, para poder trabajar con ella

for j:1 thru (length(r)-1) do block(

Desde j = 1 hasta la longitud de la matriz – 1 se realizan iteraciones sobre el valor j. En cada iteración se realizará el cálculo de un nivel de diferencias divididas, hasta llegar a la última iteración, que hará el cálculo de la última diferencia. En la primera iteración se calcularán las diferencias del tipo f[xi, x(i+1)], la última iteración calculará la diferencia dividida f[x0,…,xn]

l:makelist((r[i+1][j+1]-r[i][j+1])/(r[i+j][1]-r[i][1]),i,1,length(r)-j),

El comando realiza el cálculo para las diferencias divididas de un nivel. El número de diferencias divididas para cada nivel es igual a la longitud de la matriz menos el nivel en el que estamos. Es decir, en la primera batería de cálculos tendremos n-1 diferencias divididas, en el siguiente n-2, hasta llegar a n-j=1, que será el cálculo de f[x0,…,xn].

El primer parámetro de makelist es el cálculo a realizar, en este caso la diferencia dividida f[xi….x(j+1)]. Luego se indica la variable de iteración, que en este caso es i, y sus límites. La iteración se realizará desde i=1 hasta i<longitud de la matriz -j.

l:append(l,makelist(0,k,1,j)),

Añade a la lista que acabamos de crear con diferencias divididas el número de ceros necesarios para que alcance la misma longitud que una columna de la matriz de entrada.

r:addcol(r,l)

Se añade a la matriz de trabajo la lista obtenida con las diferencias divididas como una nueva columna.

En nuestro ejemplo las iteraciones darán como resultado:

Las cuatro iteraciones que se producirán con la matriz ejemplo

En definitiva, aunque Máxima a veces me ponga de los nervios por la dificultad para encontrar fallos en el código, su facilidad de uso de matrices y lista lo convierte en una herramienta muy útil para salir del paso en cálculos rápidos.

Código fuente para esta entrada:

http://pastebin.com/bac44VuW

y como bonus, una forma de calcular el polinomio de interpolación a partir de las diferencias divididas:

http://pastebin.com/XNEywDm4

Anuncios

Publicado el octubre 22, 2012 en Algorítmos, Maxima. Añade a favoritos el enlace permanente. Deja un comentario.

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

A %d blogueros les gusta esto: