Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

All published worksheets from http://sagenb.org

Views: 168731
Image: ubuntu2004

Interpolación de Lagrange



Recordemos el siguiente teorema conocido por el alumno.

Teorema. Sean (a1,b1),…,(an,bn)(a_1,b_1), \ldots, (a_n,b_n), nn puntos distintos del plano donde las primeras coordenadas son todas dintintas. ai≠aja_i\not = a_j si i≠ji\neq j. Entonces existe un único polinomio ff de grado a lo sumo n−1n-1 tal que f(ai)=bif(a_i)=b_i, i=1,…,ni=1,\ldots, n.

Vamos a ver que este teorema no es más que el Teorema Chino de los Restos: la aplicación K[x]⟶Kf↦f(ai)\begin{matrix}\mathbb{K}[x]&\longrightarrow & \mathbb{K}\\ f&\mapsto & f(a_i)\end{matrix} es un homomorfismo suprayectivo de anillos cuyo núcleo es el ideal (x−ai)(x-a_i) y, por tanto, K[x](x−ai)≃K.f+(x−ai)∼f(ai)\begin{matrix}\frac{\mathbb{K}[x]}{(x-a_i)} &\simeq &\mathbb{K}.\\ f+(x-a_i)&\sim&f(a_i)\end{matrix}

Dados nn puntos (a1,b1),…,(an,bn)(a_1,b_1),\ldots, (a_n,b_n), el problema de buscar un polinomio ff tal que f(ai)=bif(a_i)=b_i, i=1,…,ni=1,\ldots, n equivale a saber si el punto (b1,…,bn)(b_1,\ldots, b_n) está en la imagen de la aplicación: K[x]→K[x](x−a1)×⋯×K[x](x−a1)f↦(f+(x−a1),…,f+(x−an))=(f(a1),…,f(an)).\begin{matrix}\displaystyle\mathbb{K}[x] &\rightarrow &\displaystyle\frac{\mathbb{K}[x]}{(x-a_1)}\times\cdots \times \frac{\mathbb{K}[x]}{(x-a_1)}&\\ f & \mapsto & \big(f+(x-a_1),\ldots,f+(x-a_n)\big)&=(f(a_1),\ldots, f(a_n))\end{matrix}.

La primera parte del teorema chino de los restos, afirma que, si los ideales (x−ai)(x-a_i) son comaximales dos a dos entonces siempre hay solución. Notese que (x−ai)(x-a_i) y (x−aj)(x-a_j) son comaximales si y solo si ai≠aja_i\neq a_j, pues, en ese caso, 1=x−ai−(x−aj)aj−ai∈(x−ai,x−aj).1=\frac{x-a_i -(x-a_j)}{a_j-a_i}\in (x-a_i,x-a_j).

Por otro lado, la segunda parte del teorema chino de los restos afirma que el núcleo del homomorfismo anterior es el ideal  a=((x−a1)⋯(x−an))\mathbf{a}=((x-a_1)\cdots (x-a_n)). Así, tenemos el isomorfismo: K[x]((x−a1)⋯(x−an))→K[x](x−a1)×⋯×K[x](x−a1)f+a↦(f+(x−a1),…,f+(x−an))\begin{matrix}\displaystyle\frac{\mathbb{K}[x]}{((x-a_1)\cdots(x-a_n))} &\rightarrow &\displaystyle\frac{\mathbb{K}[x]}{(x-a_1)}\times\cdots \times \frac{\mathbb{K}[x]}{(x-a_1)}\\ f +\mathbf{a}& \mapsto & \big(f+(x-a_1),\ldots,f+(x-a_n)\big)\end{matrix} por lo que hay una única clase FF módulo a\mathbf{a} del problema. Esta clase tiene un único representante dado por un polinomio de grado a lo más n−1n-1, que se obtiene de dividir cualquier representante de FF por (x−a1)⋯(x−an)(x-a_1)\cdots (x-a_n).

Además, en este caso, la solución dada tipicamente en el teorema de interpolación de Lagrange es la misma que define el teorema Chino de los restos. Sea: fi(x)=(x−a1)⋯(x−ai−1)(x−ai+1)⋯(x−an)(ai−a1)⋯(ai−ai−1)(ai−ai+1)⋯(ai−an)f_i(x) = \frac{(x-a_1)\cdots(x-a_{i-1})(x-a_{i+1})\cdots (x-a_n)}{(a_i-a_1)\cdots(a_i-a_{i-1})(a_i-a_{i+1})\cdots (a_i-a_n)} estos polinomios cumplen que: fi(ak)={0sii≠k1sii=kf_i(a_k) = \left\{\begin{matrix}0 &\textrm{si}& i\neq k\\1&\textrm{si}& i=k\end{matrix} \right.

por lo que, la solución buscada es: f=b1f1(x)+…+bnfn(x)f = b_1f_1(x)+\ldots +b_nf_n(x)

K = QQ['x']
Kfield = K.fraction_field()
x = K.gen()

def lagrange_interpolation(puntos, valores):
    n = len(puntos)
    polinomio_nucleo = prod([x-puntos[i] for i in range(n)])
    interpolador = 0
    for i in range(n):
        pol_base = polinomio_nucleo // (x-puntos[i])
        pol_base = pol_base * valores[i] / pol_base(x = puntos[i])
        interpolador = interpolador + pol_base
    return interpolador

Cálculo de determinantes usando el Teorema Chino de los restos

Supongamos que tenemos una matriz M=(a11…a1n⋮⋱⋮an1…ann)M = \begin{pmatrix}a_{11}&\ldots&a_{1n}\\\vdots&\ddots&\vdots\\a_{n1}&\ldots &a_{nn}\end{pmatrix} cuyas entradas son polinomios aij∈Q[x]a_{ij}\in\mathbb{Q}[x].

Si calculamos el determinante de MM usando el método de Gauss tenemos el problema de que el grado de las funciones racionales que aparecen en la eliminación gausiana crece exponenciamente con el tamaño nn de la matriz.

%cython

from sage.all import copy

def haz_ceros(M,k):
    # Hacer ceros en la matriz M debajo del elemento k,k
    # Devuelve la matriz y un elemento 0 o 1 si se ha hecho un intercambio de filas
    n = M.nrows()
    # buscamos un pivote
    pivote = k
    while pivote< n and M[pivote,k] ==0:
        pivote += 1
    if pivote == n:
        #No hay nada que reducir
        return M, 0
    #Modificamos la matriz
    if pivote != k:
        M.swap_rows(pivote, k)
        swap = 1
    else:
        swap = 0
    for i in range(k+1,n):
        f = ~M[k,k]
        M.add_multiple_of_row(i, k, -M[i,k]*f,start_col=k )
    return M, swap

def escalona(M):
    n = M.nrows()
    swap = 0
    for k in range(n-1):
        M, s = haz_ceros(M,k)
        swap = swap + s
    return M, swap

El código anterior escalona una matriz dada por filas y también devuelve el número de intercambios que hemos hecho en las filas de la matriz. Apliquémoslo a unos ejemplos

M = matrix(Kfield, 3, 3, lambda i,j: ZZ[x].random_element(1))
M
\newcommand{\Bold}[1]{\mathbf{#1}}\left(−1010x+3x−1−4x−x+15xx+1−x+1x−74\begin{array}{rrr} -10 & 10 x + 3 & x - 1 \\ -4 x & -x + 1 & 5 x \\ x + 1 & -x + 1 & x - 74 \end{array}\right)
escalona(M)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(−1010x+3x−10−4x2−115x+1−25x2+275x00−192x3+14665x2+157x−74110−4x2−115x+1\begin{array}{rrr} -10 & 10 x + 3 & x - 1 \\ 0 & -4 x^{2} - \frac{11}{5} x + 1 & -\frac{2}{5} x^{2} + \frac{27}{5} x \\ 0 & 0 & \frac{-\frac{19}{2} x^{3} + \frac{1466}{5} x^{2} + 157 x - \frac{741}{10}}{-4 x^{2} - \frac{11}{5} x + 1} \end{array}\right), 0\right)
M = matrix(Kfield, 7, 7, lambda i,j: ZZ[x].random_element(2))
M
\newcommand{\Bold}[1]{\mathbf{#1}}\left(−4x2+x−1−2x+1−x2+33x2+25x−22x2+2x+1−x2−x3x2x2+8x0−5x2−x+1−2x2−x+3x2−x+1−3x+112x2−13x−46x2−xx2+2x−12x2+x−2713x2+x+1x2−x−31−x2+x+1x+2x2+x3x2+14x−14x2−20x2−2x−2x2−xx2−2x+1−x2−11x+23x−42x2+x+2−23x+260x−1−12x2+x−1x2−4x−x2−1−x2−x2−x−82x2−4−x+6−xx2+4x2−8x−x2+16x2−2x−41x2+x+42x2+1−x2+10x+12−4x2+16x−7\begin{array}{rrrrrrr} -4 x^{2} + x - 1 & -2 x + 1 & -x^{2} + 3 & 3 x^{2} + 25 x - 22 & x^{2} + 2 x + 1 & -x^{2} - x & 3 x^{2} \\ x^{2} + 8 x & 0 & -5 x^{2} - x + 1 & -2 x^{2} - x + 3 & x^{2} - x + 1 & -3 x + 1 & 12 x^{2} - 13 x \\ -46 x^{2} - x & x^{2} + 2 x - 1 & 2 x^{2} + x - 27 & 13 x^{2} + x + 1 & x^{2} - x - 31 & -x^{2} + x + 1 & x + 2 \\ x^{2} + x & 3 x^{2} + 14 x - 1 & 4 x^{2} - 20 & x^{2} - 2 x & -2 x^{2} - x & x^{2} - 2 x + 1 & -x^{2} - 11 x + 2 \\ 3 x - 4 & 2 x^{2} + x + 2 & -23 x + 26 & 0 & x - 1 & -12 x^{2} + x - 1 & x^{2} - 4 x \\ -x^{2} - 1 & -x^{2} & -x^{2} - x - 8 & 2 x^{2} - 4 & -x + 6 & -x & x^{2} + 4 \\ x^{2} - 8 x & -x^{2} + 16 & x^{2} - 2 x & -41 x^{2} + x + 4 & 2 x^{2} + 1 & -x^{2} + 10 x + 12 & -4 x^{2} + 16 x - 7 \end{array}\right)
N=escalona(M)[0]
matrix(ZZ,7,7,lambda i,j: max(N[i][j].numerator().degree(),N[i][j].denominator().degree()))
\newcommand{\Bold}[1]{\mathbf{#1}}\left(2122222034444400666650008888000010101000000121200000014\begin{array}{rrrrrrr} 2 & 1 & 2 & 2 & 2 & 2 & 2 \\ 0 & 3 & 4 & 4 & 4 & 4 & 4 \\ 0 & 0 & 6 & 6 & 6 & 6 & 5 \\ 0 & 0 & 0 & 8 & 8 & 8 & 8 \\ 0 & 0 & 0 & 0 & 10 & 10 & 10 \\ 0 & 0 & 0 & 0 & 0 & 12 & 12 \\ 0 & 0 & 0 & 0 & 0 & 0 & 14 \end{array}\right)

Es crecimiento anterior en el grado de los polinomios es lineal. Este crecimiento es mucho más acusado si incluimos denominadores

M = matrix(Kfield, 7, 7, lambda i,j: Kfield.random_element(3))
M
\newcommand{\Bold}[1]{\mathbf{#1}}\left(x3−1−116x3+52x2+711x−1−16x3+x2−12x−5213x2−12x−148−21x3−x2−2x3−x2−13x−x3−43x2+x+11915x3+12x2−x−1−19x3+x2−14x+12−13x3−12x2−5x+36x3+32x2+6x−x3+12x−1x3+34x2−12x+1x3+12x2−13x−2−x3−33x2−2−2x3+6x2+1283−2x3+x2−x+12−3x3−13x2+2x+13−x3+285x2+493x−4x3−9x2−35x−1−14x3+x2−x+132x2−45x−2−x3+3x2+x+3−2x3−2x2−225x3−x−1−x3+176x2−23x−12170x3−516x−13−12x2+1x313x3+124x2−x+5−x3−3101x2−13x−118x3−17x2+5−x3+13x2−x+32x3−4x2+x+16x3+2x2−x815x2+x−15x3+x−122x3−x2−13x3−114x2−x−1−413x3−x2−x−1x2+12x+17−x2+15x−1x3−14x2+12x−136−2x3+12x2−x+114x3+6x2+12x−19x2+16x12x3+2x−33x3+4x−12x3+12x2+12−12x2+12x+123x2x2−x+123x3+17x2+2x+15x3+x2+23x−12x3+3x2−23x−1−25x3−5x2−x+14−3x3+32x2+2x+2−x−5−2x−1x3−2x2−x+112x23x−1−12x2+5x+12−x3−112x−211x3+x2−79x−50x2−1x3−32x2+3x−1x3+12x2+2x+4−12x3−14x2−2419x−1−72x2−x+15−23x3−x2−x−1936x3−14x2+12x+2−x3−13x−14x3−x2+13x+12−x3+32x2−3x+23x3−17x38x3−311x−12−3x3+x2+2x+314−x3−5x+12−7x2+x−111x+2−32x3−23x2−x+135−32x3+x2+3x+972x2+9x+2−13x3+14x+1412x3−x−113x3+112x2−17−12x3274x3−3x2+1096x+111−19x2−13x+3677x3+16x2−2x−3519x3−x3+732x2+1x3+12x2+x−3x3−x2+2x+1312x3+14x2+67x+3\begin{array}{rrrrrrr} \frac{x^{3} - 1}{-\frac{1}{16} x^{3} + \frac{5}{2} x^{2} + \frac{7}{11} x - 1} & \frac{-\frac{1}{6} x^{3} + x^{2}}{-\frac{1}{2} x - \frac{5}{2}} & \frac{\frac{1}{3} x^{2} - \frac{1}{2} x - \frac{1}{48}}{-21 x^{3} - x^{2} - 2} & \frac{x^{3} - x^{2} - \frac{1}{3} x}{-x^{3} - \frac{4}{3} x^{2} + x + \frac{1}{19}} & \frac{\frac{1}{5} x^{3} + \frac{1}{2} x^{2} - x - 1}{-\frac{1}{9} x^{3} + x^{2} - \frac{1}{4} x + \frac{1}{2}} & \frac{-\frac{1}{3} x^{3} - \frac{1}{2} x^{2} - 5 x + 3}{6 x^{3} + \frac{3}{2} x^{2} + 6 x} & \frac{-x^{3} + \frac{1}{2} x - 1}{x^{3} + \frac{3}{4} x^{2} - \frac{1}{2} x + 1} \\ \frac{x^{3} + \frac{1}{2} x^{2} - \frac{1}{3} x - 2}{-x^{3} - 33 x^{2} - 2} & \frac{-2 x^{3} + 6 x^{2} + \frac{1}{283}}{-2 x^{3} + x^{2} - x + \frac{1}{2}} & \frac{-3 x^{3} - \frac{1}{3} x^{2} + 2 x + \frac{1}{3}}{-x^{3} + \frac{28}{5} x^{2} + \frac{49}{3} x} & \frac{-4 x^{3} - 9 x^{2} - \frac{3}{5} x - 1}{-\frac{1}{4} x^{3} + x^{2} - x + 1} & \frac{\frac{3}{2} x^{2} - \frac{4}{5} x - 2}{-x^{3} + 3 x^{2} + x + 3} & \frac{-2 x^{3} - 2 x^{2} - 2}{\frac{2}{5} x^{3} - x - 1} & \frac{-x^{3} + \frac{1}{76} x^{2} - \frac{2}{3} x - \frac{1}{2}}{\frac{1}{70} x^{3} - \frac{5}{16} x - \frac{1}{3}} \\ \frac{-\frac{1}{2} x^{2} + 1}{x^{3}} & \frac{\frac{1}{3} x^{3} + \frac{1}{24} x^{2} - x + 5}{-x^{3} - \frac{3}{101} x^{2} - \frac{1}{3} x - 1} & \frac{\frac{1}{8} x^{3} - \frac{1}{7} x^{2} + 5}{-x^{3} + \frac{1}{3} x^{2} - x + 3} & \frac{2 x^{3} - 4 x^{2} + x + \frac{1}{6}}{x^{3} + 2 x^{2} - x} & \frac{\frac{81}{5} x^{2} + x}{-\frac{1}{5} x^{3} + x - \frac{1}{2}} & \frac{2 x^{3} - x^{2}}{-\frac{1}{3} x^{3} - \frac{1}{14} x^{2} - x - 1} & \frac{-\frac{4}{13} x^{3} - x^{2} - x - 1}{x^{2} + \frac{1}{2} x + \frac{1}{7}} \\ \frac{-x^{2} + \frac{1}{5} x - 1}{x^{3} - \frac{1}{4} x^{2} + \frac{1}{2} x - \frac{1}{36}} & \frac{-2 x^{3} + \frac{1}{2} x^{2} - x + 1}{\frac{1}{4} x^{3} + 6 x^{2} + \frac{1}{2} x} & \frac{-\frac{1}{9} x^{2} + \frac{1}{6} x}{\frac{1}{2} x^{3} + 2 x - 3} & \frac{3 x^{3} + 4 x}{-\frac{1}{2} x^{3} + \frac{1}{2} x^{2} + \frac{1}{2}} & \frac{-\frac{1}{2} x^{2} + \frac{1}{2} x + 1}{\frac{2}{3} x} & \frac{2 x^{2} - x + \frac{1}{2}}{3 x^{3} + \frac{1}{7} x^{2} + 2 x + \frac{1}{5}} & \frac{x^{3} + x^{2} + \frac{2}{3} x}{-\frac{1}{2} x^{3} + 3 x^{2} - \frac{2}{3} x - 1} \\ \frac{-\frac{2}{5} x^{3} - 5 x^{2} - x + \frac{1}{4}}{-3 x^{3} + \frac{3}{2} x^{2} + 2 x + 2} & \frac{-x - 5}{-2 x - 1} & \frac{x^{3} - 2 x^{2} - x + 1}{\frac{1}{2} x^{2}} & \frac{3 x - 1}{-\frac{1}{2} x^{2} + 5 x + \frac{1}{2}} & \frac{-x^{3} - \frac{1}{12} x - 2}{11 x^{3} + x^{2} - \frac{7}{9} x - 50} & \frac{x^{2} - 1}{x^{3} - \frac{3}{2} x^{2} + 3 x - 1} & \frac{x^{3} + \frac{1}{2} x^{2} + 2 x + 4}{-\frac{1}{2} x^{3} - \frac{1}{4} x^{2} - \frac{2}{419} x - 1} \\ \frac{-\frac{7}{2} x^{2} - x + \frac{1}{5}}{-\frac{2}{3} x^{3} - x^{2} - x - \frac{193}{6}} & \frac{x^{3} - \frac{1}{4} x^{2} + \frac{1}{2}}{x + 2} & \frac{-x^{3} - \frac{1}{3} x - 1}{4 x^{3} - x^{2} + \frac{1}{3} x + \frac{1}{2}} & \frac{-x^{3} + \frac{3}{2} x^{2} - 3 x + \frac{2}{3}}{x^{3} - \frac{1}{7} x} & \frac{\frac{3}{8} x^{3} - \frac{3}{11} x - \frac{1}{2}}{-3 x^{3} + x^{2} + 2 x + \frac{3}{14}} & \frac{-x^{3} - 5 x + \frac{1}{2}}{-7 x^{2} + x - 1} & \frac{11 x + 2}{-\frac{3}{2} x^{3} - \frac{2}{3} x^{2} - x + \frac{13}{5}} \\ \frac{-\frac{3}{2} x^{3} + x^{2} + 3 x + \frac{9}{7}}{2 x^{2} + 9 x + 2} & \frac{-\frac{1}{3} x^{3} + \frac{1}{4} x + \frac{1}{4}}{\frac{1}{2} x^{3} - x - 1} & \frac{\frac{1}{3} x^{3} + \frac{1}{12} x^{2} - \frac{1}{7}}{-\frac{1}{2} x^{3}} & \frac{\frac{27}{4} x^{3} - 3 x^{2} + 1096 x + \frac{1}{11}}{-\frac{1}{9} x^{2} - \frac{1}{3} x + 3} & \frac{\frac{6}{77} x^{3} + \frac{1}{6} x^{2} - 2 x - 3}{\frac{5}{19} x^{3}} & \frac{-x^{3} + \frac{7}{32} x^{2} + 1}{x^{3} + \frac{1}{2} x^{2} + x} & \frac{-3 x^{3} - x^{2} + 2 x + \frac{1}{3}}{\frac{1}{2} x^{3} + \frac{1}{4} x^{2} + \frac{6}{7} x + 3} \end{array}\right)

Escalonamos la matriz anterior usando eliminación Gaussiana y escribimos una nueva matriz que contenga en sus entradas el máximo de los grados de numerador y denominador

N=escalona(M)[0]
matrix(ZZ,7,7,lambda i,j: max(N[i][j].numerator().degree(),N[i][j].denominator().degree()))
\newcommand{\Bold}[1]{\mathbf{#1}}\left(33333330121212121212002625272727000464648470000687070000009999000000133\begin{array}{rrrrrrr} 3 & 3 & 3 & 3 & 3 & 3 & 3 \\ 0 & 12 & 12 & 12 & 12 & 12 & 12 \\ 0 & 0 & 26 & 25 & 27 & 27 & 27 \\ 0 & 0 & 0 & 46 & 46 & 48 & 47 \\ 0 & 0 & 0 & 0 & 68 & 70 & 70 \\ 0 & 0 & 0 & 0 & 0 & 99 & 99 \\ 0 & 0 & 0 & 0 & 0 & 0 & 133 \end{array}\right)

En cualquier caso, el método de Gauss permite calcular el determinante de la matriz

def determinante_a_la_Gauss(M):
    M = M.change_ring(K.fraction_field())
    n = M.nrows()
    swap = 0
    for k in range(n-1):
        M, s = haz_ceros(M,k)
        swap = swap + s
    return prod( M.diagonal())*(-1)**swap

M = random_matrix(K, 5)
M.change_ring(Kfield)
M
\newcommand{\Bold}[1]{\mathbf{#1}}\left(−13x2+1−x2−32x−2−4x2+52x−12−x2−12x+218x2+2x−1x2−7160x+11013x2−6x−15−4x2−23x−2x2+12x5x2−x13x2−6x8x2+17x−110x2+12x−13213x2−12x−125−12x−18x2−118x2+263x−1−12x2−x−19x17x2−x−x−1−45x−1−x2−2x−13−3x2−2\begin{array}{rrrrr} -\frac{1}{3} x^{2} + 1 & -x^{2} - \frac{3}{2} x - 2 & -4 x^{2} + \frac{5}{2} x - \frac{1}{2} & -x^{2} - \frac{1}{2} x + 2 & \frac{1}{8} x^{2} + 2 x - 1 \\ x^{2} - \frac{7}{160} x + \frac{1}{10} & \frac{1}{3} x^{2} - 6 x - \frac{1}{5} & -4 x^{2} - \frac{2}{3} x & -2 x^{2} + \frac{1}{2} x & 5 x^{2} - x \\ \frac{1}{3} x^{2} - 6 x & 8 x^{2} + 17 x & -\frac{1}{10} x^{2} + \frac{1}{2} x - \frac{1}{3} & \frac{2}{13} x^{2} - \frac{1}{2} x - \frac{1}{25} & -\frac{1}{2} x \\ -\frac{1}{8} x^{2} - 1 & \frac{1}{8} x^{2} + \frac{26}{3} x - 1 & -\frac{1}{2} x^{2} - x - 1 & 9 x & \frac{1}{7} x^{2} \\ -x & -x - 1 & -45 x - 1 & -x^{2} - 2 x - \frac{1}{3} & -3 x^{2} - 2 \end{array}\right)
determinante_a_la_Gauss(M)
\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1908293}{131040} x^{10} + \frac{16511935163}{11980800} x^{9} + \frac{2392628250307}{251596800} x^{8} + \frac{7027749079727}{628992000} x^{7} + \frac{15439832494061}{628992000} x^{6} - \frac{9746337780659}{235872000} x^{5} + \frac{346857133271}{37739520} x^{4} + \frac{10418523289}{9828000} x^{3} + \frac{965758711}{7862400} x^{2} - \frac{379361}{54000} x - \frac{164}{375}
determinante_a_la_Gauss(M) - M.det()
\newcommand{\Bold}[1]{\mathbf{#1}}0

En la matriz anterior, en cada columna hay elementos de grado, a lo más, 22, luego el determinante solo puede tener grado a lo más 2⋅5=102\cdot 5=10.

Si evaluo la matriz en 1111 elementos r1,…,r11r_1,\ldots,r_{11} de K\mathbb{K} distintos, hay un único polinomio ff de grado a lo más 10 tal que f(ri)=det(M(i))f(r_i) = det(M(i)). Como el determinante es un polinomio que cumple lo anterior, tenemos que el determinante lo podemos calcular con la interpolación de Lagrange/Teorema Chino de los Restos.

def determinante_chino(M):
    n = M.nrows()
    d = [M.row(i).denominator() for i in range(n)]
    M = diagonal_matrix(d)*M
    M = M.change_ring(K)
    cota = sum([max([termino.degree() for termino in M.row(i)]) for i in range(n)])
    puntos = range(cota+1)
    valores = []
    for i in puntos:
        Mi = M(i)
        Mi = Mi.change_ring(QQ)
        DET = determinante_a_la_Gauss(Mi)
        valores.append(DET)
    DET = lagrange_interpolation(puntos, valores)/prod(d)
    return DET

M = random_matrix(K,3)
M
\newcommand{\Bold}[1]{\mathbf{#1}}\left(−12x+1−54x2−219x+12−12x2+12x−134x2+2x−112−12x2−215x18x2+136x−1−23x2+15x+523x2−x+1x+4\begin{array}{rrr} -\frac{1}{2} x + 1 & -\frac{5}{4} x^{2} - \frac{2}{19} x + \frac{1}{2} & -\frac{1}{2} x^{2} + \frac{1}{2} x - \frac{1}{34} \\ x^{2} + 2 x - \frac{1}{12} & -\frac{1}{2} x^{2} - \frac{2}{15} x & 18 x^{2} + \frac{1}{36} x - 1 \\ -\frac{2}{3} x^{2} + \frac{1}{5} x + 5 & \frac{2}{3} x^{2} - x + 1 & x + 4 \end{array}\right)
determinante_chino(M) - M.det()
\newcommand{\Bold}[1]{\mathbf{#1}}0
determinante_a_la_Gauss(M) - M.det()
\newcommand{\Bold}[1]{\mathbf{#1}}0

Comparemos los dos algoritmos para calcular determinantes en matrices aleatorias de distintos tamaños

Caso de entradas polinomios

for i in range(2,13):
    M = random_matrix(K, i)
    M = M.change_ring(Kfield)
    print('matriz de tamaño: '+str(i))
    timeit('determinante_chino(M)',repeat=1,number=1)
    timeit('determinante_a_la_Gauss(M)',repeat=1,number=1)
matriz de tamaño: 2 1 loops, best of 1: 9.3 ms per loop 1 loops, best of 1: 437 µs per loop matriz de tamaño: 3 1 loops, best of 1: 14.1 ms per loop 1 loops, best of 1: 1.07 ms per loop matriz de tamaño: 4 1 loops, best of 1: 23.9 ms per loop 1 loops, best of 1: 2.24 ms per loop matriz de tamaño: 5 1 loops, best of 1: 46.1 ms per loop 1 loops, best of 1: 4.3 ms per loop matriz de tamaño: 6 1 loops, best of 1: 71.7 ms per loop 1 loops, best of 1: 8.56 ms per loop matriz de tamaño: 7 1 loops, best of 1: 123 ms per loop 1 loops, best of 1: 15.3 ms per loop matriz de tamaño: 8 1 loops, best of 1: 186 ms per loop 1 loops, best of 1: 24.4 ms per loop matriz de tamaño: 9 1 loops, best of 1: 272 ms per loop 1 loops, best of 1: 155 ms per loop matriz de tamaño: 10 1 loops, best of 1: 400 ms per loop 1 loops, best of 1: 494 ms per loop matriz de tamaño: 11 1 loops, best of 1: 551 ms per loop 1 loops, best of 1: 2.09 s per loop matriz de tamaño: 12 1 loops, best of 1: 756 ms per loop 1 loops, best of 1: 10.2 s per loop

Observamos que, en esta implementación, para matrices hasta tamaño 9-10 y entradas polinomios de grado 2 con coeficientes pequeños es mejor usar el método de Gauss, sin embargo, a partir de tamaño 11, el coste del algoritmo de Gauss se vuelve prohibitivo, mientras que el método de interpolación todavía es usable.

Vamos a hora a dibujar un gráfico que represente el caso peor del tiempo de ejecución en una matriz de polinomios. Compararemos

  • El método de Gauss
  • El método del Teorema Chino de los restos
  • Los algoritmos internos considerada la matriz como matriz de polinomios
  • Los algoritmos internos considerada la matriz como matriz de funciones racionales
from time import time

tiempo_Gauss1           = [0 for i in range(12)]
tiempo_chino1           = [0 for i in range(20)]
tiempo_sage_polinomio1  = [0 for i in range(20)]
tiempo_sage_racional1   = [0 for i in range(7)]

for j in range(3):
    for i in range(20):
        M = random_matrix(K, i)
        t = time()
        dummy = M.det()
        t = time() -t
        tiempo_sage_polinomio1[i] = max(t, tiempo_sage_polinomio1[i])
        M = M.change_ring(Kfield)
        if i < 12:
            t = time()
            dummy = determinante_a_la_Gauss(M)
            t = time() - t
            tiempo_Gauss1[i] = max(t, tiempo_Gauss1[i])
        if i < 7:
            t = time()
            dummy = M.det()
            t = time() - t
            tiempo_sage_racional1[i] = max(t, tiempo_sage_racional1[i])            
        t = time()
        dummy = determinante_chino(M)
        t = time() - t
        tiempo_chino1[i] = max(t, tiempo_chino1[i])

list_plot(tiempo_Gauss1, color='blue',plotjoined=True) + list_plot(tiempo_chino1, color='red', plotjoined=1) + list_plot(tiempo_sage_polinomio1, color='green', plotjoined=1) + list_plot(tiempo_sage_racional1, color='orange', plotjoined=1)

En azul el tiempo que tarda el método de Gauss

En rojo el tiempo que tarda el método de interpolación

En verde, el tiempo que tarda el comando interno del determinante en Sage, como polinomios

En naranja, el tiempo que tarda el comando interno del determinante en Sage, como funciones racionales

Como conlusión, no es sorprendente ver que el mejor método es el determinante interno cuando las entradas son polinomios. Sin embargo, sorprende lo pésimo que es el método por defecto sobre la misma matriz pero considerada en Q(x)\mathbb{Q}(x) en lugar de Q[x]\mathbb{Q}[x]. Esto es debido a que no hay un método específico para Q(x)\mathbb{Q}(x) y está usando un método libre de divisiones que vale para cualquier cuerpo.

De los dos métodos que hemos programado, para matries de mayor tamaño es mucho mejor la interpolación.

Veamos que ocurre con matrices con entradas funciones racionales.

tiempo_Gauss = [0 for i in range(11)]
tiempo_chino = [0 for i in range(13)]
tiempo_sage  = [0 for i in range(7)]

for j in range(3):
    for i in range(13):
        M = random_matrix(Kfield, i)
        if i < 11:
            t = time()
            dummy = determinante_a_la_Gauss(M)
            t = time() - t
            tiempo_Gauss[i] = max(t, tiempo_Gauss[i])
        if i < 7:
            t = time()
            dummy = M.det()
            t = time() - t
            tiempo_sage[i] = max(t, tiempo_sage[i])            
        t = time()
        dummy = determinante_chino(M)
        t = time() - t
        tiempo_chino[i] = max(t, tiempo_chino[i])

list_plot(tiempo_Gauss, color='blue',plotjoined=True) + list_plot(tiempo_chino, color='red', plotjoined=1) + list_plot(tiempo_sage, color='green', plotjoined=1)

Aquí se puede apreciar que el método interno para calcular determinantes en Q(x)\mathbb{Q}(x) debe ser modificado.