Trazadores Cúbicos De Tipo 1 Con Python

por | 23 noviembre, 2008

Hola, en general un trazador o función spline en un intervalo [a,b], en matemáticas, son polinomios que están definidos en subintervalos de [a,b] que cumplen ciertas condiciones de regularidad; normalmente cumplen lo que se llama ser función de clase C^n([a,b]) , es decir, que son continuas tanto la función como sus derivadas hasta orden n en [a,b], lo cual es bastante obvio en polinomios.

Los trazadores se utilizan cuando por ejemplo tenemos un conjunto de puntos y queremos encontrar la función que pasa por todos esos puntos y que mejor se aproxima, dicha función a calcular es lo que llamamos trazador o función spline. Esta técnica de los trazadores es muy útil en informática para dibujar figuras no regulares. Por ejemplo, se puede utilizar para dibujar un perro o un edificio con tan sólo dar unos cuantos puntos representativos de dichas figuras y buscando su trazador.

El uso de trazadores en informática es muy útil ya que nos permite encontrar funciones que pasan por una cantidad finita de puntos mediante funciones a trozos, todas ellas polinómicas de grado pequeño, con lo cual el cómputo de operaciones de los algoritmos se reduce drásticamente, además de que dichas funciones a trozos (los trazadores) son funciones “suaves”.

Otro campo matemático parecido a los splines son los NURBS, los cuales nos permiten representar figuras complicadas y libres mediante interacción con el usuario mediante el movimiento de los nodos. Son muy utilizados en campos como la animación, y también se les conoce como B-splines, por ser un caso particular de los que veremos aquí. Si habéis programado con OpenGL sabréis de lo que estoy hablando.
NURB

Pero dejemos de lado estas profundas disquisiciones y profundicemos un poco en los trazadores cúbicos, los cuales son un buen comienzo para adentrarse en este apasionante mundo de la aproximación de figuras libres utilizando las matemáticas. Empecemos definiendo que es una función spline en general.

Sea \Delta=\{ a=x_0 <  x_1   < .... <  x_n= b \} una partición del intervalo [a,b]. S_{\Delta} : [a,b]\longrightarrow\mathbb{R} es una función spline de orden k\in\mathbb{N} asociada a \Delta si S_{\Delta}\in C^{k-1} ([a,b]) y si S_{\Delta} coincide en cada intervalo [ x_i , x_{i+1} ] , i=0,1,...,n-1, con un polinomio de grado \le  k . Notemos que para k=1 se obtiene una función lineal a trozos, que todos llamamos poligonal. En éste artículo me centraré sólo en el caso k=3, conocidos como trazadores cúbicos o splines cúbicas. Así pues nuestras splines serán de clase C^2 y de grado menor o igual que 3. Y en concreto a las de tipo I que definiré en breve. Dado y = ( y_0 , y_1 , .... y_n ) \in\mathbb{R}^{n+1} denotaremos por S_\Delta (y, .) a una función spline cúbica de interpolación que verifica S_\Delta (y,x_i ) = y_i para i = 0,1,...,n-1. Hay que recalcar que estas condiciones no determinan de forma única la spline cúbica. Si os estáis preguntando qué pinta el y es muy sencillo, en nuestros problemas partiremos de una serie de puntos concretos, por ejemplo, consideremos que queremos interpolar los puntos (0,1), (1,2), (4,0); pues bien, cada abscisa de los puntos representan la partición y las ordenadas de cada punto representan la y, que lógicamente estos son los extremos de cada subintervalo de cada función a trozo de que se compone nuestra función spline. Para determinar únivocamente las funciones spline se imponen, en general, 3 tipos de condiciones. Como en éste artículo sólo me voy a centrar en las de tipo I, su condición de unicidad va a ser: S^{\prime\prime}_\Delta (y,a) = S^{\prime\prime}_\Delta (y,b) = 0. Se suelen denominar condiciones naturales.
Ya hablaré de los otros dos tipos en otra ocasión, a modo de información os diré que las de tipo III se utilizan para interpolar funciones periódicas de periodo (b-a).
Para calcular explícitamente los trazadores definiremos primero los momentos de S_{\Delta}(y,.) mediante un vector M cuyas componentes vienen expresadas como M_j = S^{\prime\prime}_{\Delta}(y,x_{j}) para j=0,1,…,n-1.
Denotemos h_j = x_{j+1}-x_j j=0,1,…,n-1 Para calcular los momentos utilizaremos loa siguiente expresión:
\mu_j M_{j-1} + 2M_j + \lambda_j M_{j+1} = d_j  1\le j \le n-1
Donde los vectores \mu,\lambda,d vienen definidos por las expresiones:
\lambda_j = \frac{h_{j+1}}{h_j + h_{j+1}}
\mu_j = 1-\lambda_j
d_j = \frac{6}{h_j + h_{j+1}} \cdot ( \frac{y_{j+1}-y_j}{h_{j+1}}- \frac{y_j-y_{j-1} }{h_j} )
para 1 \le j \le n-1
Para los extremos (j=0, j=n) imponemos las condiciones de tipo I S^{\prime\prime}_\Delta (y,a) =  S^{\prime\prime}_\Delta(y,b) = 0 y obtenemos :
M_0 = M_n = 0
\lambda_0 = d_0 = \mu_n = d_n = 0
Si reunimos todas las fórmulas obtenidas tenemos que llegamos a un sistema lineal AM = d donde A es una matriz tridiagonal en la que su diagonal principal esta formada todos su elementos por el 2, su diagonal superior lo forman el vector \lambda y su diagonal inferior esta formada por el vector \mu .
Con todos estos datos ya podemos pasar a realizar un programa en Python sobre un caso concreto para ver cómo funciona. Las demostraciones de cómo se obtienen las fórmulas anteriores las podéis obtener en cualquier libro de texto universitario que se dedique a la rama matemática de Aproximación Numérica bajo el apartado de interpolación polinómica bajo splines o trazadores.
EJEMPLO
Vamos a determinar la función spline cúbica de tipo I que interpola los puntos:
(0,1), (1,2),(3,0)
Para ello he adaptado el método de trazadores cúbicos de tipo I a un sencillo programa escrito bajo lenguaje Python utilizando las librerías gráficas Matplotlib para poder dibujar la función spline. Su código es el siguiente:

  1. from Numeric import *
  2. from math import *
  3. from numpy.linalg import *
  4. from pylab import *
  5.  
  6. vectorx = array((0.,1.,3.), Float)
  7. vectory = array((1.,2.,0.), Float)
  8. h = array((0.,0.), Float)
  9. lamda = array((0.0, 0.0),Float)
  10. mu = array((0.0, 0.0),Float)
  11. d = array((0.0, 0.0,0.0),Float)
  12. A = array(([0.0, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0]), Float)
  13.  
  14. i = 0
  15. n = 0
  16.  
  17. while i&lt;len(h) :
  18.     h[i] = vectorx[i+1]-vectorx[i]
  19.     i=i+1
  20. i=1
  21. while i&lt;len(h) :
  22.     lamda[i]=h[i]/(h[i-1]+h[i])
  23.     i=i+1
  24.  
  25. i = 1
  26. while i&lt;len(h) :
  27.     mu[i-1]=1.-lamda[i]
  28.     i=i+1
  29.  
  30. i=1
  31. while i&lt;len(h) :
  32.     d[i] = (6./(h[i-1]+h[i]))*((vectory[i+1]-vectory[i])/h[i]-(vectory[i]-vectory[i-1])/h[i-1])
  33.     i=i+1
  34.  
  35. i=0
  36. n=1
  37. while i&lt;len(mu) :
  38.     A[i][n]=lamda[n-1]
  39.     i=i+1
  40.     n=n+1
  41.  
  42. i=1
  43. n=0
  44. while n&lt;len(mu) :
  45.     A[i][n]=mu[i-1]
  46.     i=i+1
  47.     n=n+1
  48.  
  49. i=0
  50. while i&lt;=len(mu) :
  51.     A[i][i]=2.
  52.     i=i+1
  53.  
  54.  
  55. M = solve(A,d)
  56.  
  57.  
  58. def funcionS1() :
  59.     x = arange(0,3,0.1)
  60.     m = len(x)
  61.     s1 = array(([0.,0.]),Float)
  62.     s2 = array(([0.0,0.]),Float)
  63.     S = zeros(m, Float64)
  64.     i = 0
  65.     l = 0
  66.     while i&lt;=1 :
  67.         while l&lt;m :
  68.             s1[i]=(vectory[i+1]-vectory[i])/(h[i+1])-(1./6.)*(2*M[i]+M[i+1])*h[i+1]
  69.             s2[i]=(1./6.)*((M[i+1]-M[i])/h[i+1])
  70.             S[l] = vectory[i]+s1[i]*(x[l]-vectorx[i])+0.5*M[i]*pow(x[l]-vectorx[i],2.)+s2[i]*pow(x[l]-vectorx[i],3.)
  71.             l = l+1
  72.         i = i+1
  73.         plot(x, S, linewidth=1.0)
  74.         xlabel('Abcisa')
  75.         ylabel('Ordenada')
  76.         title('Metodo I De Los Splines ')
  77.         grid(True)
  78.         axhline(linewidth=1, color='r')
  79.         axvline(linewidth=1, color='r')
  80.         show()
  81.  
  82. funcionS1()

El código del programa no tiene mayor misterio, a destacar la instrucción para resolver el sistema lineal AM=d, que es:

  1. M = solve(A,d)

He utilizado la función solve de las librerías Numpy, en concreto, las específicas al álgebra lineal, y que para importarlas a nuestros programas se debe hacer de una forma especial al inicio de nuestro código mediante la instrucción:

  1. from numpy.linalg import *

He decidido utilizar la función solve porque nuestro sistema lineal es muy pequeño, pero en casos grandes es aconsejable utilizar métodos de cálculo de sistemas lineales con matrices de coeficientes tridiagonales, para así disminuir drásticamente el número de operaciones que realice nuestro código.
De igual forma cuando calculamos los puntos interpoladores de nuestra función spline, si la cantidad de puntos a calcular es grande, es aconsejable utilizar el Algoritmo de Horner para su cálculo, lo cual es posible gracias a que nuestras funciones spline son funciones polinómicas definidas a trozos.
El resultado del programa lo podéis ver en la siguiente imagen:
Trazador Cúbico Tipo I

Espero pronto poder escribir sobre los trazadores cúbicos de tipo II y III

Saludos :-h

8 pensamientos en “Trazadores Cúbicos De Tipo 1 Con Python

  1. Jose

    Me gustaria saber como hiciste para poder visualizar el resultado de el calculo que hiciste con el programa.
    Estoy haciendo con ese tipo de trazadores una aplicacion de computacion grafica y tu algoritmo es una manera rapida de intentar visualizar el avance del trabajo

  2. Cristobal

    Hola Jose, la verdad es que viene muy bien, lo hago cargando en el código las librerías Matplotlib al principio con la línea

    # from pylab import *

    Y luego dentro de la función llamada funcionS1 declaro el vector x con arange, que va de 0 a 3 con una diferencia de 0.1
    Para las ordenadas son las componentes del vector S, así dentro de esa misma función con plot dibujo la gráfica, junto con las demás líneas que le siguen para poner nombres a los ejes, colores etc, sin olvidar poner la instrucción show() para poder verla.
    En este blog bajo la etiqueta python encontrarás un algoritmo de la bisección en donde explico esto con más detalle.
    Veo que usas Mac, yo no tengo , pero no creo tengas problemas de instalar la Matplotlib en tu Mac, te la puedes bajar desde su web oficial, fácil de encontrar con Google.
    Si usas Ubuntu la tienes en los repositorios oficiales.

    Me agrada saber de alguien más que yo que se interesa por estas cosas

  3. Sanatas

    Te parecerá bonito, me has hecho desempolvar los apuntes de álgebra y de cálculo para leer esta entrada y menos mal que la programación la llevo al día (aunque no en python) :-*

    P.D. ¿Porque será que estas entradas son las que más me gustan de todo el blog?

  4. Cristobal

    Jejejeje, Sanatas ahora yo también tengo curiosidad por saber el misterio de que te gusten este tipo de entradas
    Siento que tengáis que desempolvar las mates, pero es que si no programo así como que no tiene vidilla la cosa
    Además, así mato dos pájaros de un tiro

  5. Jose

    Primero que todo gracias por responder, ya termine el semestre en la Universidad pero me quedo gustando el hecho de hacer computacion grafica con trazadores cubicos, podriamos ponernos en contacto para llegar a hacer algo.
    Saludos

  6. Elis

    Hola, quería saber cómo se puede hacer en Matlab el gráfico para un trazador cúbico libre y por otro lado sujeto para un conjunto de datos con las condiciones rqueridas de las derivadas primeras.
    Gracias

  7. carlos

    amigo me gustaria que me ayudaras, trato de realizar el mismo codigo para 18 puntos; ke debo modificar para ke sirva

  8. Cristobal

    Carlos me resulta difícil saber qué has hecho
    Bueno, a ver si no meto la pata, a modificar:
    vectorx, vectory->Debes poner los puntos, en el vectorx las abcisas y en el otro las ordenadas.
    La matriz A debe ser 18×18 nula
    lamda, mu y d deberás cambiarles su tamaño a 17
    En la funciónS1 donde pone x=arange debes cambiar los dos primeros números por tus abcisas máxima y mínima.
    Creo que con eso debe funcionar, y si falta algo no creo que sea ya mucho

Los comentarios están cerrados.