5.1 Regla del Trapecio para Integración numérica con Python

Referencia: Chapra 21.1 p621 pdf645, Burden 4.2 p192 pdf202, Rodríguez 7.1.1 p273

La integración numérica con la regla del trapecio usa un polinomio de primer grado como aproximación de f(x) entre los extremos a y b del intervalo.

Es la primera de las formulas cerradas de Newton-Cotes.

I = \int_a^b f(x) dx \cong \int_a^b f_1 (x) dx

usando el rango entre [a,b] el polinomio se aproxima con una línea recta:

f_1 (x) = f(a) + \frac{f(b)-f(a)}{b-a} (x-a)

el área bajo esta línea recta es una aproximación de la integral de f(x) entre los límites a y b.

El resultado del integral es la regla del trapecio:

I = (b-a) \frac{f(a)+f(b)}{2}

que se interpreta como la multiplicación entre la base y altura promedio de un trapecio. También se llega al resultando sumando las áreas de sus componentes: rectángulo y triángulo.

El error de truncamiento se encuentra como el integral del término que le sigue al polinomio de Taylor en la aproximación, es decir el de grado 2, que al integrarlo tiene un orden de h3. (Rodríguez p275)

error_{truncar} = -\frac{h^3}{12}f''(z)

a < z < b


Algoritmo en Python


Ejemplo de Integración con el método del trapecio

Para integrar la función  en el intervalo [1,3] con 4, 16, 32 ,64 y 128 tramos,

f(x)= \sqrt {(x)} \sin(x)

Tramos = 4

La base de los trapecios (xi-xi+1), si esta base tiene valores equidistantes se sustituye por la variable h .

base =\frac{b-a}{tramos}=\frac{3-1}{4}= 0.5 = h

para cada tramo, el área de un trapecio es:

A_i = base\frac{f(x_i)+f(x_{i+1})}{2}

por lo que cada trapecio en cada subintervalo tendrá un área de:

A_1 = (0.5)\frac{f(1)+f(1.5)}{2} A_2 = (0.5)\frac{f(1.5)+f(2)}{2} A_3 = (0.5)\frac{f(2)+f(2.5)}{2} A_4 = (0.5)\frac{f(2.5)+f(3)}{2}

y el integral de todo el intervalo es la suma de todos los trapecios.

Integral = A_1+ A_2+ A_3+A_4

si se quiere tomar el factor común entre las sumas de áreas de cada tramo, también se tiene la forma de la ecuación con un h equidistante entre cada tramo del intervalo.

Integral = \frac{0.5}{2}\Big(f(1)+2f(1.5) +2f(2)+2f(2.5)+f(3)\Big)

Interpretando la fórmula,  el integral es la suma de:

  • una vez la función evaluada en cada extremo,
  • mas dos veces cada valor de la función evaluada en los puntos intermedios.
  • El resultado de la suma se multiplica por h/2.

En caso que los puntos no sean equidistantes la simplificación NO procede, y se sumaran las áreas de los trapecios de cada tramo.

Siguiendo el procedimiento, para cada cantidad de tramos en el intervalo, los resultados serán:

tramos: 4
Integral: 1.99841708623
tramos:  16
Integral:  2.05019783717
tramos:  32
Integral:  2.05277225085
tramos:  64
Integral:  2.05341563776
tramos: 128
Integral:  2.05357647096

Algoritmo en Python – Regla o Método del trapecio

Algoritmo con fórmula simplificada para varios tramos, con puntos de muestras equidistantes h.

# Integración: Regla de los trapecios
# Usando una función fx()
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: np.sqrt(x)*np.sin(x)

# intervalo de integración
a = 1
b = 3
tramos = 4

# PROCEDIMIENTO
# Regla del Trapecio
# Usando tramos equidistantes en intervalo
h = (b-a)/tramos
xi = a
suma = fx(xi)
for i in range(0,tramos-1,1):
    xi = xi + h
    suma = suma + 2*fx(xi)
suma = suma + fx(b)
area = h*(suma/2)

# SALIDA
print('tramos: ', tramos)
print('Integral: ', area)

La gráfica como ayuda didáctica, si la cantidad de tramos sea «poca», muestra los trapecios, si aumenta la cantidad de tramos, no es necesario poner líneas verticales en blanco para separar los trapecios:

# GRAFICA
# Puntos de muestra
muestras = tramos + 1
xi = np.linspace(a,b,muestras)
fi = fx(xi)
# Linea suave
muestraslinea = tramos*10 + 1
xk = np.linspace(a,b,muestraslinea)
fk = fx(xk)

# Graficando
plt.plot(xk,fk, label ='f(x)')
plt.plot(xi,fi, marker='o',
         color='orange', label ='muestras')

plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Integral: Regla de Trapecios')
plt.legend()

# Trapecios
plt.fill_between(xi,0,fi, color='g')
for i in range(0,muestras,1):
    plt.axvline(xi[i], color='w')

plt.show()

Para rellenar el área bajo la curva se usa la instrucción plt.fill_between(xi,0,fi, color=’g’). La división de los trapecios con una línea blanca (white) se realiza con la instrucción plt.axvline(xi[i], color=’w’).


Algoritmo con muestras

Algoritmo usando las muestras y para tramos que pueden tener distancias arbitrarias. En este caso se usa la formula del área del trapecio para cada tramo.

Usado para ejercicios donde los datos proporcionados son muestras. Se adapta el bloque de ingreso y el procedimiento inicia desde «Regla de Trapecio». La cantidad de tramos será el numero de muestras menos uno tramos = len(xi) - 1

# Integración: Regla de los trapecios
# Usando una muestras xi,fi
import numpy as np
import matplotlib.pyplot as plt

def integratrapecio_fi(xi,fi):
    ''' sobre muestras de fi para cada xi
        integral con método de trapecio
    '''
    n = len(xi)
    suma = 0
    for i in range(0,n-1,1):
        dx = xi[i+1]-xi[i]
        untrapecio = dx*(fi[i+1]+fi[i])/2
        suma = suma + untrapecio
    return(suma)

# PROGRAMA -----------------
# INGRESO
xi = [1. , 1.5, 2. , 2.5, 3. ]
fi = [0.84147098, 1.22167687, 1.28594075,
      0.94626755, 0.24442702]

# PROCEDIMIENTO
Area = integratrapecio_fi(xi,fi)

# SALIDA
print('tramos: ',len(xi)-1)
print('Integral con trapecio: ',Area)

4.5 Interpolación paramétrica

Referencia: Rodriguez 6.9.2 p236, Burden 9Ed 3.6 p164

En algunos casos, los datos (x,y) no tienen una relación de tipo funcional y(x), entonces no se pueden aplicar directamente los métodos de interpolación revisados.

Por ejemplo, en la trayectoria del balón en el «gol imposible», la gráfica de la trayectoria en el espacio o sus proyecciones en los planos dependen del parámetro tiempo «t»en lugar de una relación de x,y,z

Referencia: 1Eva_IT2018_T4 El gol imposible

Sin embargo si las coordenadas (x,y) se expresan como funciones de otra variable t denomiada parámetro, entonces los puntos x(t), y(t) tienen relación funcional, y se pueden construir polinomios de interpolación.


Ejemplo

Las coordinadas x(t) y y(t) del recorrido de un cohete registradas en los instantes t fueron:

ti  = [0,1,2,3]
xti = [2,1,3,4]
yti = [0,4,5,0]

Usaremos un algoritmo en Python para mostrar la trayectoria x,y para el problema planteado.

Al realizar la interpolación de los puntos para obtener polinomios que dependen de «t» se obtiene:

px = lambda t: (-2/3)*(t**3) + (7/2)*(t**2) + (-23/6)*t + 2
py = lambda t: (-1/2)*(t**3) + (9/2)*t

polinomios con los que se puede realizar la gráfica px(t), py(t) en forma separada. Pero para comprender mejor la trayectoria del cohete, se utiliza la gráfica px(t) vs py(t) en el intervalo t entre[0,3]

Las intrucciones para mostrar el resultado son:

# interpolación paramétrica
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
ti = [0,1,2,3]
xti = [2,1,3,4]
yti = [0,4,5,0]

# PROCEDIMIENTO
# interpolando con lagrange
px = lambda t: (-2/3)*(t**3) + (7/2)*(t**2) + (-23/6)*t + 2
py = lambda t: (-1/2)*(t**3) + (9/2)*t

t = np.arange(0,3,0.01)
puntosx = px(t)
puntosy = py(t)

# Salida
plt.plot(puntosx,puntosy)
plt.show()

4.4.1 Trazadores Cúbicos Natural o libre

Referencia:  Burden 9Ed 3.5 p147 pdf.151, Chapra 18.6.3 p532 pdf556, Rodriguez 6.11.1 p244

Tiene como objetivo incorporar condiciones adicionales para la función en los extremos del intervalo donde se encuentran los puntos o «nodos».

Por ejemplo, si los datos son los de una trayectoria en un experimento de física, se podría disponer de la aceleración el punto inicial de las muestras (izquierda) y de salida (derecha) del intervalo de las muestras.

El método busca obtener un polinómio de tercer grado para cada sub-intervalo o tramo entre «nodos» consecutivos (j,  j+1), de la forma:

S_j(x_j) = a_j + b_j (x-x_j) + c_j(x-xj)^2 + d_j(x-x_j)^3

para cada j = 0, 1, …, n-1

Para n datos existen n-1 tramos, cuatro incógnitas (coeficientes) que evaluar por cada tramo, como resultado 4(n-1) incógnitas para todo el intervalo .

Para los términos (xj+1xj) de los tramos que se usan varias veces en el desarrollo, se simplifican como hj:

h_j = x_{j+1} - x_{j} S_j(x_j) = a_j + b_j h_j + c_j h_j^2 + d_jh_j^3

Para generar el sistema de ecuaciones, se siguen los siguientes criteros:

1.  Los valores de la función deben ser iguales en los nodos interiores

S_j(x_{j+1}) = f(x_{j+1}) = S_{j+1}(x_{j+1})

2. Las primeras derivadas en los nodos interiores deben ser iguales

S'_j(x_{j+1}) = S'_{j+1}(x_{j+1})

3. Las segundas derivadas en los nodos interiores deben ser iguales

S''_j(x_{j+1}) = S''_{j+1}(x_{j+1})

4. El primer y último polinomio deben pasar a través de los puntos extremos

f(x_0) = S_0(x_0) = a_0 f(x_n) = S_n(x_n) = a_n

5. Una de las condiciones de frontera se satisface:

5a. frontera libre o natural: Las segundas derivadas en los nodos extremos son cero

S''(x_0) = S''(x_n) = 0

5b. frontera sujeta: las primeras derivadas en los nodos extremos son conocidas

S'(x_0) = f'(x_0)

S'(x_n) = f'(x_n)

Tarea: Revisar en los textos el planteamiento de las ecuaciones para resolver el sistema que se genera al plantear el polinomio.

Ubicar en el texto las ecuaciones resultantes, que son las que se aplicarán en el algoritmo.


Algoritmo en Python

El algoritmo parte de lo desarrollado para «trazadores lineales», donde se presentaron los bloques para:

  • construir el trazador en una tabla de polinomios por tramos,
  • graficar los trazadores en cada tramo, requientdo,
  • evaluar cada uno de ellos en cada tramo con muestras suficientes para presentar una buena precisión en la gráfica.

Del algoritmo básico se modifica entonces el bloque del cálculo de los polinomios de acuerdo al planteamiento de formulas y procedimientos para trazadores cúbicos naturales (enviado a revisar como tarea).

# Trazador cúbico natural
# Condición: S''(x_0) = S''(x_n) = 0
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

def traza3natural(xi,yi):
    n = len(xi)
    
    # Valores h
    h = np.zeros(n-1, dtype = float)
    for j in range(0,n-1,1):
        h[j] = xi[j+1] - xi[j]
    
    # Sistema de ecuaciones
    A = np.zeros(shape=(n-2,n-2), dtype = float)
    B = np.zeros(n-2, dtype = float)
    S = np.zeros(n, dtype = float)

    A[0,0] = 2*(h[0]+h[1])
    A[0,1] = h[1]
    B[0] = 6*((yi[2]-yi[1])/h[1] - (yi[1]-yi[0])/h[0])

    for i in range(1,n-3,1):
        A[i,i-1] = h[i]
        A[i,i] = 2*(h[i]+h[i+1])
        A[i,i+1] = h[i+1]
        factor21 = (yi[i+2]-yi[i+1])/h[i+1]
        factor10 = (yi[i+1]-yi[i])/h[i]
        B[i] = 6*(factor21 - factor10)
        
    A[n-3,n-4] = h[n-3]
    A[n-3,n-3] = 2*(h[n-3]+h[n-2])
    factor12 = (yi[n-1]-yi[n-2])/h[n-2]
    factor23 = (yi[n-2]-yi[n-3])/h[n-3]
    B[n-3] = 6*(factor12 - factor23)
    
    # Resolver sistema de ecuaciones S
    r = np.linalg.solve(A,B)
    for j in range(1,n-1,1):
        S[j] = r[j-1]
    S[0] = 0
    S[n-1] = 0
    
    # Coeficientes
    a = np.zeros(n-1, dtype = float)
    b = np.zeros(n-1, dtype = float)
    c = np.zeros(n-1, dtype = float)
    d = np.zeros(n-1, dtype = float)
    for j in range(0,n-1,1):
        a[j] = (S[j+1]-S[j])/(6*h[j])
        b[j] = S[j]/2
        factor10 = (yi[j+1]-yi[j])/h[j]
        c[j] = factor10 - (2*h[j]*S[j]+h[j]*S[j+1])/6
        d[j] = yi[j]
    
    # Polinomio trazador
    x = sym.Symbol('x')
    px_tabla = []
    for j in range(0,n-1,1):

        pxtramo = a[j]*(x-xi[j])**3 + b[j]*(x-xi[j])**2
        pxtramo = pxtramo + c[j]*(x-xi[j])+ d[j]
        
        pxtramo = pxtramo.expand()
        px_tabla.append(pxtramo)
    
    return(px_tabla)

# PROGRAMA -----------------------
# INGRESO , Datos de prueba
xi = np.array([0.1 , 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.8, 1.7, 2.0])
muestras = 10 # entre cada par de puntos

# PROCEDIMIENTO
# Tabla de polinomios por tramos
n = len(xi)
px_tabla = traza3natural(xi,fi)

# SALIDA
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
    print(' x = ['+str(xi[tramo-1])
          +','+str(xi[tramo])+']')
    print(str(px_tabla[tramo-1]))

con lo que se obtiene el resultado por cada tramo

Polinomios por tramos: 
 x = [0.1,0.2]
-146.666666666667*x**3 + 44.0*x**2 + 0.566666666666666*x + 1.1
 x = [0.2,0.3]
283.333333333333*x**3 - 214.0*x**2 + 52.1666666666667*x - 2.34
 x = [0.3,0.4]
-136.666666666667*x**3 + 164.0*x**2 - 61.2333333333333*x + 9.0
>>>

Para el trazado de la gráfica se añade al final del algoritmo:

# GRAFICA
# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,muestras)
    
    # evalua polinomio del tramo
    pxtramo = px_tabla[tramo-1]
    pxt = sym.lambdify('x',pxtramo)
    ytramo = pxt(xtramo)

    # vectores de trazador en x,y
    xtraza = np.concatenate((xtraza,xtramo))
    ytraza = np.concatenate((ytraza,ytramo))
    tramo = tramo + 1

# Gráfica
plt.plot(xi,fi,'ro', label='puntos')
plt.plot(xtraza,ytraza, label='trazador'
         , color='blue')
plt.title('Trazadores Cúbicos Naturales')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()

Si los polinomios no se igualan entre los nodos, tampoco sus velocidades y aceleraciones; podría considerar la experiencia semejante a variaciones de velocidad y aceleración en una trayectoria:

Car sales woman scares customers. maxman.tv. 4 enero 2016

Pilot shows off the impact of flying at 9.5G.@loteofficial. 2 Marzo2023

https://www.youtube.com/shorts/i7KU7bqLSKQ

4.4 Trazadores lineales (Splines) grado1

Referencia: Chapra 18.6.1 p525 pdf549

El concepto de trazador se originó en la técnica de dibujo que usa una cinta delgada y flexible (spline) para dibujar curvas suaves a través de un conjunto de puntos.

La unión más simple entre dos puntos es una línea recta. El método crea un polinomio para cada par de puntos consecutivos en el intervalo, por lo que el resultado será una tabla de polinomios.

Los trazadores de primer grado para un grupo de datos ordenados pueden definirse como un conjunto de funciones lineales.

f(x) = f(x_0) + m_0(x-x_0), x_0\leq x\leq x_1 f(x) = f(x_1) + m_1(x-x_1), x_1\leq x\leq x_2

f(x) = f(x_{n-1}) + m_{n-1}(x-x_{n-1}), x_{n-1}\leq x\leq x_n

donde

m_i = \frac{f(x_{i+1}) - f(x_i)}{(x_{i+1}-x_i)}

Observe que la expresión de f(x) para un tramo entre dos puntos es el polinomio de grado 1 realizado con diferencia finita avanzadas  o las diferencias divididas.

Las ecuaciones se pueden usar para evaluar la función en cualquier punto entre x0 y xn. Al localizar primero el intervalo dentro del cual está el punto, puede seleccionar el polinomio que corresponde a ese tramo.

Polinomios por tramos: 
 x = [0.1,0.2]    p(x) = 3.5*x + 1.1
 x = [0.2,0.3]    p(x) = -1.0*x + 2.0
 x = [0.3,0.4]    p(x) = 3.0*x + 0.8

Algoritmo en Python

Datos de los puntos como ejemplo para el algoritmo

xi = [0.1 , 0.2, 0.3, 0.4]
fi = [1.45, 1.8, 1.7, 2.0]

El método con trazadores lineales, permite plantear los bloques necesarios para manejar varios polinomios, uno por cada tramo entre dos puntos dentro del intervalo del problema

El método permitirá disponer de un punto de partida para trazadores de mayor grado, por ejemplo los cúbicos.

Los polinomios de cada tramo se almacenan en una tabla, cada uno puede ser utilizado individualmente en su respectivo tramo, por ejemplo para realizar la gráfica de línea entre tramos.

# Trazador (spline) lineal, grado 1
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

def trazalineal(xi,fi):
    n = len(xi)
    x = sym.Symbol('x')
    px_tabla = []
    
    tramo = 1
    while not(tramo>=n):
        # con 1ra diferencia finita avanzada 
        numerador = fi[tramo]-fi[tramo-1]
        denominador = xi[tramo]-xi[tramo-1]
        m = numerador/denominador
        pxtramo = fi[tramo-1] + m*(x-xi[tramo-1])
        px_tabla.append(pxtramo)
        tramo = tramo + 1
        
    return(px_tabla)

# PROGRAMA
# INGRESO , Datos de prueba
xi = [0.1 , 0.2, 0.3, 0.4]
fi = [1.45, 1.8, 1.7, 2.0]
muestras = 10 # entre cada par de puntos

# PROCEDIMIENTO
# Tabla de polinomios por tramos
n = len(xi)
px_tabla = trazalineal(xi,fi)

# SALIDA
print('Polinomios por tramos: ')
for tramo in range(1,n,1):
    print('  x = ['+str(xi[tramo-1])
          +','+str(xi[tramo])+']')
    print(str(px_tabla[tramo-1]))

Se obtiene como resultado:

Polinomios por tramos: 
  x = [0.1,0.2]
3.5*x + 1.1
  x = [0.2,0.3]
-1.0*x + 2.0
  x = [0.3,0.4]
3.0*x + 0.8
>>> 

Para añadir la gráfica se añaden las instrucciones para:

  • evaluar el polinomio en cada tramo
  • concatenar los resultados de todos los tramos en los vectores xtraza, ytraza.
  • poner en la gráfica los puntos del problema y las líneas que genera cada polinomio
# GRAFICA
# Puntos para graficar cada tramo
xtraza = np.array([])
ytraza = np.array([])
tramo = 1
while not(tramo>=n):
    a = xi[tramo-1]
    b = xi[tramo]
    xtramo = np.linspace(a,b,muestras)

    # evalua polinomio del tramo
    pxtramo = px_tabla[tramo-1]
    pxt = sym.lambdify('x',pxtramo)
    ytramo = pxt(xtramo)

    # vectores de trazador en x,y
    xtraza = np.concatenate((xtraza,xtramo))
    ytraza = np.concatenate((ytraza,ytramo))
    tramo = tramo + 1

# Gráfica
plt.plot(xi,fi,'o', label='puntos')
plt.plot(xtraza,ytraza, label='trazador')
plt.title('Trazadores lineales (splines)')
plt.xlabel('xi')
plt.ylabel('px(xi)')
plt.legend()
plt.show()

4.3 Interpolación polinómica de Lagrange con Python

Referencia: Chapra 18.2 pdf 540, Burden 9Ed 3.1 p108, Rodriguez 6.2 p195

El polinomio de interpolación de Lagrange es una reformulación del polinomio de interpolación de Newton que el método evita el cálculo de las diferencias divididas. El método tolera las diferencias entre distancias x de los puntos de muestra.

El polinomio de Lagrange se construye a partir de las fórmulas:

f_{n} (x) = \sum_{i=0}^{n} L_{i} (x) f(x_{i}) L_{i} (x) = \prod_{j=0, j \neq i}^{n} \frac{x-x_j}{x_i - x_j}

Donde una vez que se han seleccionado los puntos a usar que generan la misma cantidad de términos que puntos.

Ejemplo: Polinomio Interpolante de Lagrange con grado 3

Dados los 4 puntos en la tabla se requiere generar un polinomio de grado 3 de la forma:

p(x) = a_0 x^3 + a_1 x^2 + a_2 x^1 + a_3
xi 0 0.2 0.3 0.4
fi 1 1.6 1.7 2.0

Revisando la aplicación de la fórmula se tiene que calcular cuatro términos,

término 1

L_{0} (x) = \frac{(x-0.2)(x-0.3)(x-0.4)}{(0-0.2)(0-0.3)(0-0.4)}

término 2

L_{1} (x) = \frac{(x-0)(x-0.3)(x-0.4)}{(0.2-0)(0.2-0.3)(0.2-0.4)}

término 3

L_{2} (x) = \frac{(x-0)(x-0.2)(x-0.4)}{(0.3-0)(0.3-0.2)(0.3-0.4)}

término 4

L_{3} (x) = \frac{(x-0)(x-0.2)(x-0.3)}{(0.4-0)(0.4-0.2)(0.4-0.3)}

se construye el polinomio usando la fórmula para fn(x) para cada valor fi,

p_3(x) = 1 L_{0} (x) + 1.6 L_{1} (x) + 1.7 L_{2} (x) + 2 L_{3} (x) p_3(x) = 1 \frac{(x-0.2)(x-0.3)(x-0.4)}{(0-0.2)(0-0.3)(0-0.4)} + + 1.6 \frac{(x-0)(x-0.3)(x-0.4)}{(0.2-0)(0.2-0.3)(0.2-0.4)} + 1.7 \frac{(x-0)(x-0.2)(x-0.4)}{(0.3-0)(0.3-0.2)(0.3-0.4)} + 2 \frac{(x-0)(x-0.2)(x-0.3)}{(0.4-0)(0.4-0.2)(0.4-0.3)}

La simplificación de la expresión del polinomio se puede dejar como tarea.

Una forma de verificar que el polinomio es correcto es usar un punto original xi[i] y comprobar que la evaluación tiene como resultado fi[i].

xi = np.array([0, 0.2, 0.3, 0.4])
fi = np.array([1, 1.6, 1.7, 2.0])

Algoritmo en Python

En el algoritmo en Python, para construir las expresiones de cada término se usa la forma simbólica con Sympy.

En cada término Li(x) se usan todos los elementos i , excepto el mismo elemento i, en el numerador y denominador de la expresión.

En polinomio se agrupan todos los términos multiplicados por su respectivo valor fi[i].

    valores de fi:  [1.  1.6 1.7 2. ]
divisores en L(i):  [-0.024  0.004 -0.003  0.008]

Polinomio de Lagrange, expresiones
400.0*x*(x - 0.4)*(x - 0.3) 
  - 566.666666666667*x*(x - 0.4)*(x - 0.2) 
  + 250.0*x*(x - 0.3)*(x - 0.2) 
  - 41.6666666666667*(x - 0.4)*(x - 0.3)*(x - 0.2)

Polinomio de Lagrange: 
41.666666666667*x**3 - 27.5*x**2 + 6.8333333333334*x + 1.0

Las expresiones del polinomio contiene los binomios en su forma básica, para resolver y simplificar las ecuaciones se usa polinomio.expand().

Para realizar la gráfica del polinomio es conveniente convertirlo a la forma lambda con Numpy, de esta forma se evalúa en una sola linea todos los puntos para el intervalo [a,b] en x.

# Interpolacion de Lagrange
# divisoresL solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([0, 0.2, 0.3, 0.4])
fi = np.array([1, 1.6, 1.7, 2.0])

# PROCEDIMIENTO
# Polinomio de Lagrange
n = len(xi)
x = sym.Symbol('x')
polinomio = 0
divisorL = np.zeros(n, dtype = float)
for i in range(0,n,1):
    
    # Termino de Lagrange
    numerador = 1
    denominador = 1
    for j  in range(0,n,1):
        if (j!=i):
            numerador = numerador*(x-xi[j])
            denominador = denominador*(xi[i]-xi[j])
    terminoLi = numerador/denominador

    polinomio = polinomio + terminoLi*fi[i]
    divisorL[i] = denominador

# simplifica el polinomio
polisimple = polinomio.expand()

# para evaluación numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('    valores de fi: ',fi)
print('divisores en L(i): ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.show()

4.2.2 Diferencias divididas de Newton

Referencia: Chapra 18.1.3 p508 pdf532, Rodriguez 6.7 p223, Burden 9Ed 3.3 p124.

El método se usa en el caso que los puntos en el «eje x» se encuentran espaciados de forma arbitraria y provienen de una función desconocida pero supuestamente diferenciable.

La n-ésima diferencia dividida finita es:

f[x_{n}, x_{n-1}, \text{...}, x_{1}, x_{0}] = \frac{f[x_{n}, x_{n-1}, \text{...}, x_{1}]- f[x_{n-1}, x_{n-2}, \text{...}, x_{0}]}{x_{n}-x_{0}}

Para lo cual se debe interpretar la tabla de diferencias divididas, también como una aproximación a una derivada:

i xi f[xi] Primero Segundo Tercero
0 x0 f[x0] f[x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0]
1 x1 f[x1] f[x2,x1] f[x3,x2,x1]
2 x2 f[x2] f[x3,x2]
3 x2 f[x3]

En la fórmula del polinomio, las diferencias divididas sirven para evaluar los coeficientes de cada término adicional para aumentar el grado. Semejante al proceso realizado para  «diferencias finitas divididas»:

f_n(x) = f(x_0)+(x-x_0)f[x_1,x_0] + + (x-x_0)(x-x_1)f[x_2,x_1,x_0] + \text{...}+ + (x-x_0)(x-x_1)\text{...}(x-x_{n-1})f[x_n,x_{n-1},\text{...},x_0]

Este polinomio de interpolación se conoce como de Newton en diferencias divididas.


Ejercicio

El ejercicio perminte mostrar con mayor detalle las operaciones del procedimiento. Para el ejercicio, se toma como referencia el descrito en el enlace, observando que las distancias entre puntos para el «eje x» NO son equidistantes,

Referencia: 1Eva_IIT2008_T3_MN Ganancia en inversión

Los datos del problema muestran varios valores de inversión xi que generan cada uno ganancias fi.

Considere que los valores invertidos en materia prima para producción, dependen de la demanda de del producto en el mercado, motivo por el que los valores de inversión no guardan distancias equidistante entre si.

Siguiento del formato usado en los algoritmos, los datos son:

xi = np.array([3.2 , 3.8 , 4.2 , 4.5 ])
fi = np.array([5.12, 6.42, 7.25, 6.85])

Con los datos se llena la tabla de  diferencias divididas, donde por simplicidad, se escriben las operaciones en cada casilla. La última de deja como tarea.:

i xi f[xi] Primero Segundo Tercero
0 3.2 5.12 \frac{6.42-5.12}{3.8-3.2} =2.1667 \frac{2.075-2.1667}{4.2-3.2} =-0.0917 f[x3,x2,x1,x0]
1 3.8 6.42 \frac{7.25-6.42}{4.2-3.8} =2.075 \frac{-1.3333-2.075}{4.5-3.8} =-4.869
2 4.2 7.25 \frac{6.85-7.25}{4.5-4.2} =-1.3333
3 4.5 6.85

Las diferencias divididas de la primera fila son los valores usados para la expresión del polinomio de interpolacion:

ddividida = tabla[0,3:] = [ 2.1667, -0.0917, -3.6749, 0. ]

La expresión del polinomio inicia con el primer valor de la función f(x0).

p_3(x) = f_0 = 5.12

se calcula el primer término usando el factor de la primera diferencia dividida y se añade a la expresión del polinomio.

término = (x-x_0)f[x1,x0] = (x-3.2)2.1667 p_3(x) = 5.12 + 2.1667(x-3.2)

con el siguiente valor de ddividida[] se procede de manera semejante:

término = (x-x_0)(x-x_1)f[x1,x0] = = (x-3.2)(x-3.8)(-0.0917) p_3(x) = 5.12 + 2.1667*(x-3.2)+ + (-0.0917)(x-3.2)(x-3.8)

Realizando el cálculo del tercer y último termino con diferencias divididas,  el polinomio obtenido es:

p_3(x) = 5.12 + 2.1667(x-3.2) + + (-0.0917)(x-3.2)(x-3.8) + + (-3.6749)(x - 3.2)(x - 3.8)(x - 4.2)

que simplificando al multiplicar entre los términos (x-xi):

p_3(x) = 184.7569 - 149.9208 x + 41.0673 x^2 - 3.6749 x^3

El polinomio graficado en el intervalo de xi muestra que pasa por todos los puntos.


Algoritmo en Python

El algoritmo para  interpolación de diferencias divididas o Newton, considera reutilizar el procedimiento de cálculo de diferencias finitas incorporando la parte del denominador.

La creación de la expresión del polinomio también es semejante a la usada para diferencias finitas avanzadas.

Se incorpora la parte gráfica para observar los resultados en el intervalo xi, con el número de muestras = 101 para tener una buena resolución de la línea del polinomo.

# Polinomio interpolación
# Diferencias Divididas de Newton
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([3.2, 3.8, 4.2, 4.5])
fi = np.array([5.12, 6.42, 7.25, 6.85])

# PROCEDIMIENTO

# Tabla de Diferencias Divididas Avanzadas
titulo = ['i   ','xi  ','fi  ']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias divididas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('F['+str(j-2)+']')

    # cada fila de columna
    i = 0
    paso = j-2 # inicia en 1
    while (i < diagonal):
        denominador = (xi[i+paso]-xi[i])
        numerador = tabla[i+1,j-1]-tabla[i,j-1]
        tabla[i,j] = numerador/denominador
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Divididas
# caso: puntos equidistantes en eje x
dDividida = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    factor = dDividida[j-1]
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
np.set_printoptions(precision = 4)
print('Tabla Diferencia Dividida')
print([titulo])
print(tabla)
print('dDividida: ')
print(dDividida)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencias Divididas - Newton')
plt.show()

El resultado del algoritmo se muestra a continuación:

Tabla Diferencia Dividida
[['i   ', 'xi  ', 'fi  ', 'F[1]', 'F[2]', 'F[3]', 'F[4]']]
[[ 0.      3.2     5.12    2.1667 -0.0917 -3.6749  0.    ]
 [ 1.      3.8     6.42    2.075  -4.869   0.      0.    ]
 [ 2.      4.2     7.25   -1.3333  0.      0.      0.    ]
 [ 3.      4.5     6.85    0.      0.      0.      0.    ]]
dDividida: 
[ 2.1667 -0.0917 -3.6749  0.    ]
polinomio: 
2.16666666666667*x - 3.67490842490842*(x-4.2)*(x-3.8)*(x-3.2)
 - 0.0916666666666694*(x - 3.8)*(x - 3.2) - 1.81333333333334
polinomio simplificado: 
-3.67490842490842*x**3 + 41.0673076923077*x**2 
- 149.920860805861*x + 184.756923076923
>>> 

4.2.1 Diferencias finitas avanzadas – polinomio de interpolación

Referencia: Rodriguez 6.6.4 Pdf.221, Burden 9Ed p129

Se usa en interpolación cuando los puntos en el «eje x» se encuentran igualmente espaciados, la diferencia entre puntos consecutivos xi es una constante denominada h.

h = xi+1 – xi


Una relación entre derivadas y diferencias finitas se establece mediante:

f^{(n)}(z) = \frac{\Delta ^{n} f_{0}}{h^{n}} para algún k en el intervalo [x0,xn]

\frac{\Delta ^{n} f_{0}}{h^{n}} es una aproximación para la n-ésima derivada f(n) en el intervalo [x0,xn]

El polinomio de interpolación se puede contruir por medio de diferencias finitas avanzadas con las siguiente fórmula:

p_n (x) = f_0 + \frac{\Delta f_0}{h} (x - x_0) + + \frac{\Delta^2 f_0}{2!h^2} (x - x_0)(x - x_1) + + \frac{\Delta^3 f_0}{3!h^3} (x - x_0)(x - x_1)(x - x_2) + \text{...} + \frac{\Delta^n f_0}{n!h^n} (x - x_0)(x - x_1) \text{...} (x - x_{n-1})

Observe que en la medida que se toman más puntos de muestra, el grado del polinomio puede ser mayor.


Ejercicio

Se toman los datos del ejercicio de «diferencias finitas» y se grafican los resultados. Se resalta en la gráfica la distancia equidistante h entre los puntos consecutivos xi como los mostrados en la gráfica anterior.

xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

y se usan los resulados previos de diferencias finitas:

Tabla Diferencia Finita: 
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]

valores con los que se inicia el cálculo de la distancia constante entre puntos xi:

h = xi[1] – xi[0] = 0.2-0.1 = 0.1

Los valores de diferencias finitas a usar en la fórmula son de la primera fila desde la columna 3 en adelante.

dfinita = tabla[0,3:] = [0.15, -0.05, 0.25, 0.]

La construción de la expresión del polinomio empieza con el valor del primer punto de la función f(0.1)

p_3(x) = f_0 = 1.45

para añadirle el primer término j=1, hay que calcular el factor:

factor = \frac{\Delta^1 f_0}{1!h^1} = \frac{0.15}{0.1} = 1.5

con lo que el término se completa como:

término = factor(x-x_0) = 1.5(x-0.1) p_3(x) = 1.45 + 1.5(x-0.1)

Para el segundo término j=2, se repite el proceso:

factor = \frac{\Delta^2 f_0}{2!h^2} = \frac{-0.05}{2!(0.1)^2} = -2.5 término = factor(x-x_0) = -2.5(x-0.1)(x-02) p_3(x)= 1.45 + 1.5(x-0.1) +(-2.5)(x-0.1)(x-02)

Finalmente, añadiendo el tercer término j=3 cuyo cálculo es semejante a los anteriores, se deja como tarea.

El resultado del método es:

p_3(x)= 1.45 + 1.5(x-0.1) -2.5(x-0.1)(x-02) + + 41.667(x - 0.3)(x - 0.2)(x - 0.1)

Se puede seguir simpificando la respuesta, por ejemplo usando solo el término de grado con 1.5(x-0-1) se tiene que:

p_3(x) = 1.3 + 1.5x -2.5(x-0.1)(x-02)+ + 41.667(x - 0.3)(x - 0.2)(x - 0.1)

Seguir simplificando la expresión en papel, se deja como tarea.

La gráfica del polinomio obtenido comprueba que pasa por cada uno de los puntos dados para el ejercicio.


Algoritmo en Python

El polinomio se contruye usando el ejercicio de diferencias finitas.

Para contruir la expresión del polinomio añadiendo los términos de la fórmula, se define la variable simbólica x con sympy.

Para simplificar el polinomio resultante con las expresiones de multiplicación, se utiliza la instrucción sym.expand().

En caso de requerir evaluar la fórmula con un vector de datos se la convierte a la forma lambda.

Se añaden las instrucciones para realizar la gráfica en el intervalo [a,b] dado por los valores xi, definiendo el número de muestras = 101 que son suficientes para que la gráfica se observe sin distorsión.

# Polinomio interpolación
# Diferencias finitas avanzadas
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

# PROCEDIMIENTO

# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('df'+str(j-2))
    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Finitas avanzadas
# caso: puntos equidistantes en eje x
h = xi[1] - xi[0]
dfinita = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    denominador = np.math.factorial(j)*(h**j)
    factor = dfinita[j-1]/denominador
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('Tabla Diferencia Finita')
print([titulo])
print(tabla)
print('dfinita: ')
print(dfinita)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')

plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación polinómica')
plt.show()

teniendo como resultado:

Tabla Diferencia Finita
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]
dfinita: 
[ 0.15 -0.05  0.25  0.  ]
polinomio: 
1.5*x + 41.6666666666667*(x - 0.3)*(x - 0.2)*(x - 0.1) 
- 2.50000000000001*(x - 0.2)*(x - 0.1) + 1.3
polinomio simplificado: 
41.6666666666667*x**3 - 27.5*x**2 
+ 6.83333333333335*x + 0.999999999999999

el polinomio de puede evaluar como px(valor) una vez que se convirtió a lambda:

>>> px(0.1)
1.4500000000000004
>>> px(0.2)
1.6000000000000025
>>> px0.3)
1.7000000000000042
>>> 

Como tarea se recomienda realizar las gráficas comparativas entre métodos:

4.2 Diferencias finitas con Python

Referencia: Rodriguez 6.5 Pdf.211, Chapra 18.1.3 pdf.509

Establece relaciones simples entre los puntos dados que describen la funcion f(x). Las tabla de diferencias finitas es un elemento base para varios métodos de interpolación, por lo que se trata como un tema inicial.

Para un ejemplo se toman los siguientes puntos:

xi 0.10 0.2 0.3 0.4
fi 1.45 1.6 1.7 2.0

La tabla de diferencias finitas se construye tomando los datos y sus índices como parte de las primeras columnas.

i xi fi Δ1fi Δ2fi Δ3fi Δ4fi
0 0.1 1.45 0.15 -0.05 0.25 0.
1 0.2 1.6 0.1 0.2 0. 0.
2 0.3 1.7 0.3 0. 0. 0.
3 0.4 2.0 0. 0. 0. 0.

Cada casila de diferencia finita Δjfi se obtiene restando los dos valores consecutivos de la columna anterior. Por ejemplo, para la primera columna:

\Delta ^1 f_0 = 1.6-1.45 = 0.15 \Delta ^1 f_1 = 1.7-1.6 = 0.1 \Delta ^1 f_2 = 2.0-1.7 = 0.3

y la siguiente Δ1f3 no es posible calcular.

Siguiendo el mismo procedimiento se calculan los valores de la siguiente columna Δ2fi como las diferencias de la columna anterior, así sucesivamente.

Si la función f(x) de donde provienen los datos es un polinonio de grado n, entonces la n-ésima diferencia finita será una constante, y las siguientes diferencias se anularán.


Algoritmo en Python

Para crear la tabla de diferencias finitas, las primeras columnas requieren concatenar los valores de los índice ixi y fi.

xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

Los índices i se crean en un vector ki, pues la variable i es usada como fila en matrices, por lo que evitamos confundirlas al usar la variable.

A la matriz con las tres columnas descritas, se le añade a la derecha una matriz de nxn para calcular las diferencias.

Se calculan las diferencias para cada columna, realizando la operación entre los valores de cada fila. Considere que no se realizan cálculos desde la diagonal hacia abajo en la tabla, los valores quedan como cero.

Al final se muestra el título y el resultado de la tabla.

# Tabla de Diferencias finitas
# resultado en: [título,tabla]
# Tarea: verificar tamaño de vectores
import numpy as np

# INGRESO, Datos de prueba
xi = np.array([0.10, 0.2, 0.3, 0.4])
fi = np.array([1.45, 1.6, 1.7, 2.0])

# PROCEDIMIENTO

# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n  = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla   = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('df'+str(j-2))

    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i + 1

    diagonal = diagonal - 1
    j = j + 1

# SALIDA
print('Tabla Diferencia Finita: ')
print([titulo])
print(tabla)

el resultado de aplicar el algoritmo es:

Tabla Diferencia Finita: 
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[ 0.    0.1   1.45  0.15 -0.05  0.25  0.  ]
 [ 1.    0.2   1.6   0.1   0.2   0.    0.  ]
 [ 2.    0.3   1.7   0.3   0.    0.    0.  ]
 [ 3.    0.4   2.    0.    0.    0.    0.  ]]

>>> 

Gráfica de puntos con líneas horizontales

Para tener una referencia visual sobre las primeras diferencias finitas, en una gráfica se trazan las líneas horizontales que pasan por cada punto. Para las segundas diferencias se debe graficar las primeras diferencias finitas vs xi repitiendo el proceso. Las líneas de distancia se añadieron con un editor de imágenes.

las intrucciones adicionales al algoritmo anterior para añadir la gráfica son:

# Gráfica
import matplotlib.pyplot as plt

for i in range(0,n,1):
    plt.axhline(fi[i],ls='--', color='yellow')

plt.plot(xi,fi,'o', label = 'Puntos')

plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencia Finita')
plt.show()

4.1 El polinomio de interpolación con Sympy-Python

Referencia: Rodriguez 6.1 p191, Chapra 18.3 p520 pdf544

Unir n puntos en el plano usando un polinomios de interpolación se puede realizar con un polinomio de grado (n-1) y la misma cantidad de puntos en el plano.

Por ejemplo, dados los 4 puntos en la tabla se requiere generar un polinomio de grado 3 de la forma:

p(x) = a_0 x^3 + a_1 x^2 + a_2 x^1 + a_3
xi 0 0.2 0.3 0.4
fi 1 1.6 1.7 2.0

Es posible generar un sistema de ecuaciones en base p(x) que pase por cada uno de los puntos. Por ejemplo si se toma el primer punto con xi=0 y fi=1 se establece la ecuación:

p(0) = 1 a_0 (0)^3 + a_1 (0)^2 + a_2 (0)^1 + a_3 = 1

Note que ahora las incógnitas son los coeficientes ai. Luego se continúa con los otros puntos seleccionados hasta completar las ecuaciones necesarias para el grado del polinomio seleccionado.

a_0 (0.0)^3 + a_1 (0.0)^2 + a_2 (0.0)^1 + a_3 = 1.0 a_0 (0.2)^3 + a_1 (0.2)^2 + a_2 (0.2)^1 + a_3 = 1.6 a_0 (0.3)^3 + a_1 (0.3)^2 + a_2 (0.3)^1 + a_3 = 1.7 a_0 (0.4)^3 + a_1 (0.4)^2 + a_2 (0.4)^1 + a_3 = 2.0

El sistema obtenido se resuelve con alguno de los métodos conocidos.

Los métodos requieren usar la matriz A y vector B:

\begin{pmatrix} 0.0^3 & 0.0^2 & 0.0^1 & 1 \\ 0.2^3 & 0.2^2 & 0.2^1 & 1 \\ 0.3^3 & 0.3^2 & 0.3^1 & 1 \\ 0.4^3 & 0.4^2 & 0.4^1 & 1 \end{pmatrix} . \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{pmatrix} = \begin{pmatrix} 1.0 \\ 1.6 \\ 1.7 \\ 2.0 \end{pmatrix}

La matriz A también se conoce como Matriz Vandermonde D, se construye observando que los coeficientes se elevan al exponente referenciado al índice columna pero de derecha a izquierda. La última columa da 1 por tener como coeficiente el valor de xi0.

Por los métodos para resolver sistemas de ecuaciones, se observa que hay que realizar la matriz aumentada, pivotear por filas, etc.

Para la solución se propone un algoritmo realizado en Python, que para enfocarnos en la interpolación, se obtiene el siguente resultado:

los coeficientes del polinomio: 
[[ 41.66666667]
 [-27.5       ]
 [  6.83333333]
 [  1.        ]]

con lo que se contruye el polinomio con los valores ai.

Polinomio de interpolación p(x): 
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333333*x + 1.0

que re-escrito es

p(x) = 41.6666 x^3 - 27.5 x^2 + 6.8333 x + 1

Polinomio que se grafica dentro del intervalo de los datos del problema y se obtiene la línea que pasa por los puntos.


Algoritmo en Python

Para desarrollar el ejercicio, se realiza un bloque para construir la matriz A y vector B, usando los vectores que representan los puntos de muestra de la función o experimento.

xi = [0,0.2,0.3,0.4]
fi = [1,1.6,1.7,2.0]

Observe que en éste caso los datos se dan en forma de lista y se deben convertir hacia arreglos.

La matriz A o Matriz Vandermonde D, se construye observando que los coeficientes del vector xi se elevan al exponente referenciado al índice columna pero de derecha a izquierda.

Construir el polinomio consiste en resolver el sistema de ecuaciones indicado en la sección anterior. Para simplificar la solución del sistema, se usa numpy para algebra lineal, que entrega el vector solución que representan los coeficientes del polinomio de interpolación.

Para contruir la expresión del polinomio, se usa la forma simbólica con Sympy, de forma semejante a la usada para construir el polinomio de Taylor.

Para facilitar la evación del polinomio de forma numérica, se convierte el polinomio a la forma Lambda. Con un número de muestras suficientes para suavizar la curva dentro del intervalo [a,b] del problema se puede realizar la gráfica del polinomio.

Finalmente, con la gráfica se comprueba que la curva pase por todos los puntos del intervalo.

# El polinomio de interpolación
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
xi = [0,0.2,0.3,0.4]
fi = [1,1.6,1.7,2.0]
# muestras = tramos+1
muestras = 101

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xi)
B = np.array(fi)
n = len(xi)

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype =float)
for i in range(0,n,1):
    for j in range(0,n,1):
        potencia = (n-1)-j # Derecha a izquierda
        D[i,j] = xi[i]**potencia

# Aplicar métodos Unidad03. Tarea
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)

# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i in range(0,n,1):
    potencia = (n-1)-i   # Derecha a izquierda
    termino = coeficiente[i]*(x**potencia)
    polinomio = polinomio + termino

# Polinomio a forma Lambda
# para evaluación con vectores de datos xin
px = sym.lambdify(x,polinomio)

# Para graficar el polinomio en [a,b]
a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)

# Usando evaluación simbólica
##yin = np.zeros(muestras,dtype=float)
##for j in range(0,muestras,1):
##    yin[j] = polinomio.subs(x,xin[j])
    
# SALIDA
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)

# Grafica
plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xin,yin, label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()

con lo que se obtiene el siguiente resultado:

Matriz Vandermonde: 
[[0.    0.    0.    1.   ]
 [0.008 0.04  0.2   1.   ]
 [0.027 0.09  0.3   1.   ]
 [0.064 0.16  0.4   1.   ]]
los coeficientes del polinomio: 
[[ 41.66666667]
 [-27.5       ]
 [  6.83333333]
 [  1.        ]]
Polinomio de interpolación: 
41.6666666666667*x**3 - 27.5*x**2 + 6.83333333333333*x + 1.0

formato pprint
                  3         2                           
41.6666666666667*x  - 27.5*x  + 6.83333333333333*x + 1.0

Para mostrar el polinomio de una manera más fácil de interpretar se usa la instrucción sym.pprint(), usada al final del algoritmo.

3.6.1 Método de Gauss-Seidel con Python

Referencia: Ejemplo 01 Chapra Ejemplo 11.3  p311 pdf335

\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}

El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fué implementado el procedimiento.


Planteamiento

Para el desarrollo con lápiz y papel, se despeja una variable de cada ecuación, se empieza con la primera expresión para obtener x0

3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 3 x_0 = 7.85 + 0.1 x_1 + 0.2 x_2 x_0 = \frac{7.85 + 0.1 x_1 + 0.2 x_2}{3}

con la segunda ecuación se obtiene x1

0.1 x_0 +7x_1 -0.3 x_2 = -19.3 7x_1 = -19.3 - 0.1 x_0 +0.3 x_2 x_1 = \frac{ -19.3 - 0.1 x_0 +0.3 x_2}{7}

usando la tercera ecuación se obtiene x2

0.3 x_0 - 0.2 x_1 + 10 x_2 = 71.4 10 x_2 = 71.4 - 0.3 x_0 + 0.2 x_1 x_2 = \frac{71.4 - 0.3 x_0 + 0.2 x_1}{10}

Vector inicial

Como el vector inicial no se indica en el enunciado, se considera usar el vector de ceros para iniciar  las iteraciones.

X = [0,0,0]

con tolerancia de 0.00001


Iteraciones

Con las ecuaciones obtenidas en el planteamiento, se desarrolla usando el vector inicial presentado, hasta que el |error|<tolera

itera = 0

x_0 = \frac{7.85 + 0.1 (0) + 0.2 (0)}{3} = 2.61 x_1 = \frac{ -19.3 - 0.1 (2.61) +0.3 (0)}{7} = -2.79 x_2 = \frac{71.4 - 0.3 (2.61) + 0.2 (-2.79)}{10} = 7.00

X = [2.61, -2.79, 7.00]

diferencias = [2.61, -2.79, 7.00] – [0,0,0]
diferencias = [2.61, -2.79, 7.00]
error = ||diferencias|| , se usa la norma de max(diferencias)
error = 7.00

itera = 1

X = [2.61, -2.79, 7.00]

x_0 = \frac{7.85 + 0.1 (-2.79) + 0.2 (7.00)}{3} = 2.99 x_1 = \frac{ -19.3 - 0.1 (2.99) +0.3 (7.00)}{7} = -2.49 x_2 = \frac{71.4 - 0.3 (2.99) + 0.2 (-2.49)}{10} = 7.00

X = [2.99, -2.49, 7.00]

diferencias = [2.99, -2.49, 7.00] – [2.61, -2.79, 7.00]

diferencias = [0.38, 0.3 , 0. ]

error = ||diferencias||, se usa la norma de max(diferencias)

error = 0.38

itera = 2

X = [2.99, -2.49, 7.00]

x_0 = \frac{7.85 + 0.1 (-2.49) + 0.2 (7.00)}{3} = 3.00 x_1 = \frac{ -19.3 - 0.1 (3.00) +0.3 (7.00)}{7} = -2.5 x_2 = \frac{71.4 - 0.3 (3.00) + 0.2 (-2.5)}{10} = 7.00

X = [3.00, -2.50, 7.00]

diferencias = [3.00, -2.50, 7.00] – [2.99, -2.49, 7.00]

diferencias = [ 0.01, -0.01, 0. ]

error = ||diferencias||, se usa la norma de max(diferencias)

error = 0.01

El error aún es mayor que tolera, por lo que es necesario continuar con las iteraciones.

Observaciones

El error disminuye en cada iteración, por lo que el método converge hacia la raiz.


Algoritmo con Python

Con la descripción dada para el método de Gauss-Seidel, se muestra una forma básica de implementar el algoritmo.

Ejemplo 01 Chapra Ejemplo 11.3 p311 pdf335

\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}

El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fué implementado el procedimiento. Para generalizar el algoritmo, incluya como tarea aumentar el procedimiento de pivoteo por filas.

# Método de Gauss-Seidel
# solución de sistemas de ecuaciones
# por métodos iterativos

import numpy as np

# INGRESO
A = np.array([[3. , -0.1, -0.2],
              [0.1,  7  , -0.3],
              [0.3, -0.2, 10  ]])

B = np.array([7.85,-19.3,71.4])

X0  = np.array([0.,0.,0.])

tolera = 0.00001
iteramax = 100

# PROCEDIMIENTO

# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
#  valores iniciales
X = np.copy(X0)
diferencia = np.ones(n, dtype=float)
errado = 2*tolera

itera = 0
while not(errado<=tolera or itera>iteramax):
    # por fila
    for i in range(0,n,1):
        # por columna
        suma = 0 
        for j in range(0,m,1):
            # excepto diagonal de A
            if (i!=j): 
                suma = suma-A[i,j]*X[j]
        
        nuevo = (B[i]+suma)/A[i,i]
        diferencia[i] = np.abs(nuevo-X[i])
        X[i] = nuevo
    errado = np.max(diferencia)
    itera = itera + 1

# Respuesta X en columna
X = np.transpose([X])

# revisa si NO converge
if (itera>iteramax):
    X=0
# revisa respuesta
verifica = np.dot(A,X)

# SALIDA
print('respuesta X: ')
print(X)
print('verificar A.X=B: ')
print(verifica)

que da como resultado:

respuesta X: 
[[ 3. ]
 [-2.5]
 [ 7. ]]
verificar A.X=B: 
[[  7.84999999]
 [-19.3       ]
 [ 71.4       ]]
>>> 

que es la respuesta del problema obtenida durante el desarrollo del ejemplo con otros métodos.


Tarea

Completar el algoritmo para usar una matriz diagonal dominante, usando al menos el pivoteo parcial por filas.