Publicado en por Edison Del Rosario

s1Eva_IT2019_T1 Oxígeno y temperatura en agua

Ejercicio: 1Eva_IT2019_T1 Oxígeno y temperatura en agua

Literal a

Se requiere un polinomio de grado 3 siendo el eje x correspondiente a temperatura. Son necesarios 4 puntos de referencia alrededor de 15 grados, dos a la izquierda y dos a la derecha.

Se observa que los datos en el eje x son equidistantes, h=8, y ordenados en forma ascendente, se cumple con los requisitos para usar diferencias finitas avanzadas. que tiene la forma de:

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{...}

Tabla

xi f[xi] f[x1,x0] f[x2,x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0]
8 11.5 9.9-11.5=
-1.6		 -1.5-(-1.6) =
0.1		 0.4-0.1=
0.3		 —
16 9.9 8.4-9.9=
-1.5		 -1.1-(1.5)=
0.4		 — —
24 8.4 7.3-8.4=
-1.1		 — — —
32 7.3 — — — —

Con lo que el polinomio buscado es:

p_3 (x) = 11.5 + \frac{-1.6}{8} (x - 8) + + \frac{0.1}{2!8^2} (x - 8)(x - 16) + \frac{0.3}{3!8^3} (x - 8)(x - 16)(x - 24)

Resolviendo y simplificando el polinomio, se puede observar que al aumentar el grado, la constante del término disminuye.

p_3(x)=12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3

para el cálculo del error se puede usar un término adicional del polinomio, añadiendo un punto más a la tabla de diferencia finitas. Se evalúa éste término y se estima el error que dado que el término de grado 3 es del orden de 10-5, el error será menor. (Tarea)

p_3(15)=12.9- 0.15 (15) - 0.00390625 (15)^2 + 0.00009765625 (15)^3

Evaluando el polinomio en temperatura = 15:

p3(15) = 10.1006835937500

literal b

se deriva el polinomio del literal anterior y se evalua en 16:

p'_3(x)=0- 0.15 - 0.00390625 (2) x + 0.00009765625 (3)x^2 p'_3(16)=0- 0.15 - 0.00390625 (2)(16) + 0.00009765625 (3)(16)^2

p’3(16) = -0.20

literal c

El valor de oxígeno usado como referencia es 9, cuyos valores de temperatura se encuentran entre 16 y 24 que se toman como rango inicial de búsqueda [a,b]. Por lo que el polinomio se iguala a 9 y se crea la forma estandarizada del problema:

p_3(x)=9 9 = 12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3 12.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3 -9 = 0 f(x) = 3.9- 0.15 x - 0.00390625 x^2 + 0.00009765625 x^3

Para mostrar el procedimiento se realizan solo tres iteraciones,

1ra Iteración
a=16 , b = 24, c = (16+24)/2 = 20
f(a) = 0.9, f(b) = -0.6, f(c) = 0.011
error = |24-16| = 8
como f(c) es positivo, se mueve el extremo f(x) del mismo signo, es decir a

2da Iteración
a=20 , b = 24, c = (20+24)/2 = 22
f(a) = 0.119, f(b) = -0.6, f(c) = -0.251
error = |24-20|= 4
como f(c) es negativo, se mueve el extremo f(x) del mismo signo, b

3ra Iteración
a=20 , b = 22, c = (20+22)/2 = 21
f(a) = 0.119, f(b) = -0.251, f(c) = -0.068
error = |22-20| = 2
como f(c) es negativo, se mueve el extremo f(x) del mismo signo, b
y así sucesivamente hasta que error< que 10-3

Usando el algoritmo en python se obtendrá la raiz en 20.632 con la tolerancia requerida.

Revisión de resultados

Usando como base los algoritmos desarrollados en clase:

Literal a
[[ 1.   8.  11.5 -1.6  0.1  0.3  0. ]
 [ 2.  16.   9.9 -1.5  0.4  0.   0. ]
 [ 3.  24.   8.4 -1.1  0.   0.   0. ]
 [ 4.  32.   7.3  0.   0.   0.   0. ]]
9.765625e-5*x**3 - 0.00390625*x**2 - 0.15*x + 12.9
p(15) =  10.1006835937500
Literal b
0.00029296875*x**2 - 0.0078125*x - 0.15
dp(16) = -0.200000000000000

literal c
[ i, a, c, b, f(a), f(c), f(b), paso]
1 16.000 20.000 24.000 0.900 0.119 -0.600 8.000 
2 20.000 22.000 24.000 0.119 -0.251 -0.600 4.000 
3 20.000 21.000 22.000 0.119 -0.068 -0.251 2.000 
4 20.000 20.500 21.000 0.119 0.025 -0.068 1.000 
5 20.500 20.750 21.000 0.025 -0.022 -0.068 0.500 
6 20.500 20.625 20.750 0.025 0.001 -0.022 0.250 
7 20.625 20.688 20.750 0.001 -0.010 -0.022 0.125 
8 20.625 20.656 20.688 0.001 -0.004 -0.010 0.062 
9 20.625 20.641 20.656 0.001 -0.002 -0.004 0.031 
10 20.625 20.633 20.641 0.001 -0.000 -0.002 0.016 
11 20.625 20.629 20.633 0.001 0.001 -0.000 0.008 
12 20.629 20.631 20.633 0.001 0.000 -0.000 0.004 
13 20.631 20.632 20.633 0.000 0.000 -0.000 0.002 
14 20.632 20.632 20.633 0.000 0.000 -0.000 0.001 
raiz:  20.63232421875

Algoritmos Python usando la funcion de interpolación y un procedimiento encontrado en:

http://blog.espol.edu.ec/analisisnumerico/5-1-1-diferencias-finitas-avanzadas-polinomio/

http://blog.espol.edu.ec/analisisnumerico/2-1-1-biseccion-ejemplo01/

# 1Eva_IT2019_T1 Oxígeno y temperatura en mar
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

def interpola_dfinitas(xi,fi):
    '''
    Interpolación de diferencias finitas avanzadas
    resultado: polinomio en forma simbólica
    '''
    # Tabla de diferencias finitas
    titulo = ['i','xi','fi']
    n = len(xi)
    # cambia a forma de columnas
    i = np.arange(1,n+1,1)
    i = np.transpose([i])
    xi = np.transpose([xi])
    fi = np.transpose([fi])
    # Añade matriz de diferencias
    dfinita = np.zeros(shape=(n,n),dtype=float)
    tabla = np.concatenate((i,xi,fi,dfinita), axis=1)
    # Sobre matriz de diferencias, por columnas
    [n,m] = np.shape(tabla)
    c = 3
    diagonal = n-1
    while (c<m):
        # Aumenta el título para cada columna
        titulo.append('df'+str(c-2))
        # calcula cada diferencia por fila
        f = 0
        while (f < diagonal):
            tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
            f = f+1
        
        diagonal = diagonal - 1
        c = c+1

    # POLINOMIO con diferencias finitas
    # caso: puntos en eje x equidistantes
    dfinita = tabla[:,3:]
    n = len(dfinita)
    x = sym.Symbol('x')
    h = xi[1,0]-xi[0,0]
    polinomio = fi[0,0]
    for c in range(1,n,1):
        denominador = np.math.factorial(c)*(h**c)
        factor = dfinita[0,c-1]/denominador
        termino=1
        for f  in range(0,c,1):
            termino = termino*(x-xi[f])
        polinomio = polinomio + termino*factor
    # simplifica polinomio, multiplica los (x-xi)
    polinomio = polinomio.expand()
    
    return(tabla, polinomio)

# INGRESO
tm = [0.,8,16,24,32,40]
ox = [14.6,11.5,9.9,8.4,7.3,6.4]

xi = [8,16,24,32]
fi = [11.5,9.9,8.4,7.3]

# PROCEDIMIENTO
x = sym.Symbol('x')
# literal a
tabla, polinomio = interpola_dfinitas(xi,fi)
p15 = polinomio.subs(x,15)
# literal b
deriva = polinomio.diff(x,1)
dp16 = deriva.subs(x,16)

px =  sym.lambdify(x,polinomio)
xk = np.linspace(np.min(xi),np.max(xi))
pk = px(xk)

# SALIDA
print('Literal a')
print(tabla)
print(polinomio)
print('p(15) = ',p15)
print('Literal b')
print(deriva)
print('dp(16) =', dp16)

# gráfica
plt.plot(tm,ox,'ro')
plt.plot(xk,pk)
plt.axhline(9,color="green")
plt.xlabel('temperatura')
plt.ylabel('concentracion de oxigeno')
plt.grid()
plt.show()

# --------literal c ------------

# Algoritmo de Bisección
# [a,b] se escogen de la gráfica de la función
# error = tolera

# se convierte forma de símbolos a numéricos
buscar = polinomio-9
fx = sym.lambdify(x,buscar)

# INGRESO
a = 16
b = 24
tolera = 0.001

# PROCEDIMIENTO
tabla = []
tramo = b-a

fa = fx(a)
fb = fx(b)
i = 1
while (tramo > tolera):
    c = (a+b)/2
    fc = fx(c)
    tabla.append([i,a,c,b,fa,fc,fb,tramo])
    i = i+1
                 
    cambia = np.sign(fa)*np.sign(fc)
    if (cambia<0):
        b = c
        fb = fc
    else:
        a=c
        fa = fc
    tramo = b-a
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
tabla = np.array(tabla)

raiz = c

# SALIDA
print('\n literal c')
np.set_printoptions(precision = 4)
print('[ i, a, c, b, f(a), f(c), f(b), paso]')
# print(tabla)

# Tabla con formato
n=len(tabla)
for i in range(0,n,1):
    unafila = tabla[i]
    formato = '{:.0f}'+' '+(len(unafila)-1)*'{:.3f} '
    unafila = formato.format(*unafila)
    print(unafila)
    
print('raiz: ',raiz)