s1Eva_2022PAOI_T2 Capacidad de alimentos para pacientes internos

Ejercicio: 1Eva_2022PAOI_T2 Capacidad de alimentos para pacientes internos

a) Realice el planteamiento del sistema de ecuaciones que permita determinar la cantidad máxima de pacientes de cada grupo que podrían ser atendidos usando todos los productos disponibles. Una vez planteadas las ecuaciones, se le indica que la capacidad de atención para emergencia sea fija en K = 10 pacientes (variable libre).

Producto\ Paciente Maternidad Pos – operatorio Covid_19 emergencia Suministro diario
Producto A 0.2 0.1 1.7 0.25 135
Producto B 0.5 2 0.05 0.4 320
Producto C 1.5 0.2 0.75 1.4 410

Sistema de ecuaciones usando la columna de emergencias como variable libre y valor constante para la variable igual a K=10

0.2 x_1 + 0.1 x_2 +1.7 x_3 = 135-0.25 K 0.5 x_1 + 2 x_2 +0.05 x_3 = 320-0.4 K 1.5 x_1 + 0.2 x_2 +0.75 x_3 = 410-1.4 K

b) Muestre los pasos detallados para la matriz aumentada y pivoteo parcial por filas.

matriz aumentada

\begin{pmatrix} 0.2 & 0.1 & 1.7 & 135-0.25*10 \\ 0.5 & 2 &0.05 & 320-0.4*10 \\ 1.5 & 0.2 & 0.75 &410-1.4*10 \end{pmatrix}

pivoteo parcial por filas

\begin{pmatrix} 1.5 & 0.2 & 0.75 & 396 \\ 0.5 & 2 &0.05 & 316 \\ 0.2 & 0.1 & 1.7 & 132.5 \end{pmatrix}

c) Desarrolle al menos 3 iteraciones para el método requerido, con expresiones completas.
1.5 x_1 + 0.2 x_2 +0.75 x_3 = 396

0.5 x_1 + 2 x_2 +0.05 x_3 = 316 0.2 x_1 + 0.1 x_2 +1.7 x_3 = 132.5

ecuaciones a usar

x_1 = \frac{1}{1.5}(396 - 0.2 x_2 - 0.75 x_3) x_2 = \frac{1}{2}(316 - 0.5 x_1 - 0.05 x_3) x_3 = \frac{1}{1.7}(132.5 - 0.2 x_1 - 0.1 x_2)

Dado que el valor de la variable libre se establecía con K=10, el vector inicial podría ser el doble de este valor

X_0 = [20,20,20]

itera=1

X_0 = [20,20,20] x_1 = \frac{1}{1.5}(396 - 0.2(20) - 0.75(20)) =251.33 x_2 = \frac{1}{2}(316 - 0.5(20) - 0.05 (20)) = 152.5 x_3 = \frac{1}{1.7}(132.5 - 0.2 (20) - 0.1 (20)) =74.41 diferencia = [231.33, 132.5, 54.41] errado = max |[231.33, 132.5, 54.41]| = 231.33

itera=2

X_1 = [251.33, 152.5, 74.41] x_1 = \frac{1}{1.5}(396 - 0.2(152.5) - 0.75(74.41)) =206.46 x_2 = \frac{1}{2}(316 - 0.5(251.33) - 0.05 (74.41)) = 96.30 x_3 = \frac{1}{1.7}(132.5 - 0.2 (251.33) - 0.1 (152.5)) =39.40 diferencia = [-44.86, -59.27, -211.92] errado = max |[-44.86, -59.27, -211.92]| = 211.92

itera=3

X_2 = [231.85, 105.39, 48.16] x_1 = \frac{1}{1.5}(396 - 0.2(105.39) - 0.75(48.16)) =231.85 x_2 = \frac{1}{2}(316 - 0.5(231.85) - 0.05 (48.16)) = 105.39 x_3 = \frac{1}{1.7}(132.5 - 0.2 ( 231.85) - 0.1 (105.39)) =48.16 diferencia = [25.39, 121.09, -158.29] errado = max |[25.39, 121.09, -158.29]| = 158.29 X_3 = [231.85, 105.39, 48.16]

d) Realice las observaciones necesarias sobre los errores entre iteraciones y la convergencia.

El error entre iteraciones disminuye, el método converge a:

respuesta X: 
[[228.22]
 [ 99.81]
 [ 44.92]]

e) Si se decide no atender a los pacientes del grupo emergencias, ¿Qué aumento individual de cada una de otros grupos de pacientes podría soportarse con la cantidad diaria de alimento disponible? (use el algoritmo.py).

Se interpreta como K=0 y usando el algoritmo se obtiene:

respuesta X: 
[[237.23]
 [ 99.54]
 [ 45.64]]

el aumento de pacientes entre grupos es

diferencia = [9.01, -0.27, 0.72]

en términos de pacientes, como número entero, solo se gana 9 pacientes para el primer grupo. Se descarta la parte decimal si los pacientes no se cuentan en miles, por lo que las otras variaciones se interpretan como cero.

Algoritmo en Python

Datos tomados luego de pivoteo parcial por filas

Resultados:

0 [20. 20. 20.]
0 [251.33333333 152.5         74.11764706]
1 [206.60784314  93.31372549  39.10784314]
2 [232.00424837 105.37034314  47.85121107]
3 [226.02501538  98.80265763  44.15418589]
4 [228.74921937 100.38989151  45.24395951]
5 [227.99270138  99.68159617  44.83009822]
6 [228.2940714   99.8810722   44.96076477]
7 [228.20214132  99.80246303  44.91357559]
8 [228.23621714  99.82662528  44.92901496]
9 [228.22527582  99.81772034  44.92358473]
10 [228.22917825  99.82059143  44.92539577]
11 [228.22788993  99.81957054  44.92476777]
12 [228.22834004  99.81990832  44.92497939]
13 [228.2281892   99.8197905   44.92490656]
14 [228.22824132  99.81983004  44.92493124]
numero de condicion: 2.166985328561448
respuesta con Jacobi
[[228.22824132]
 [ 99.81983004]
 [ 44.92493124]]
verificando:
[[396.00002641]
 [316.00002729]
 [132.00001438]]

Algoritmo en Python

# 1Eva_2022PAOI_T2 Capacidad de alimentos para pacientes internos
import numpy as np

def jacobi(A,B,tolera,X,iteramax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    xnuevo = np.copy(X)

    itera = 0
    print(itera, X)
    while not(errado<=tolera or itera>iteramax):
        
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            xnuevo[i] = nuevo
        diferencia = np.abs(xnuevo-X)
        errado = np.max(diferencia)
        print(itera, xnuevo)
        X = np.copy(xnuevo)
        itera = itera + 1
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=itera
    return(X)


# INGRESO
A = np.array([[1.5, 0.2, 0.75],
              [0.5, 2.0, 0.05],
              [0.2, 0.1, 1.7 ]],
             dtype=float)

B = np.array([396.,316.,132.],
             dtype=float)
tolera = 1e-4

X = np.array([20.,20.,20.],
             dtype=float)

# PROCEDIMIENTO

# numero de condicion
ncond = np.linalg.cond(A)

respuesta = jacobi(A,B,tolera,X)

verifica = np.dot(A,respuesta)

# SALIDA
print('numero de condicion:', ncond)
print('respuesta con Jacobi')
print(respuesta)
print('verificando:')
print(verifica)

s1Eva_2022PAOI_T1 Impacto en trayectoria del drone

Ejercicio: 1Eva_2022PAOI_T1 Impacto en trayectoria del drone

Desarrollo analítico

a) Realice el planteamiento del problema usando inicialmente las trayectorias en el eje x, donde para el intervalo de operación del misil antidrone, se observa más de un impacto.

x1(t) = x2(t)

f(t) = cos(t) –  sin(0.75 t) =0

y1(t) = y2(t)

sin(2 t) =kt

k = \frac{sin(2 t)}{t}

b) Usando el método de Newton-Raphson encuentre el valor de t en el cual se pretende realizar el impacto al drone. Realice al menos 3 iteraciones de forma analítica, use tolerancia de 10-4,

f(t) = cos(t) - sin(0.75 t) f'(t) = - sin(t) - 0.75 cos(0.75 t)

Como punto inicial para encontrar la raíz de f(t) podría ser t0=4 para el punto marcado en rojo. Para el método de Newton-Raphson se tiene que

t_{i+1} = t_i - \frac{f(t_i)}{f'(t_i)} error = |t_{i+1} - t_i|

iteración 1 t0=4

t_1 = 4 - \frac{cos(4) - sin(0.75*4)}{- sin(4) - 0.75 cos(0.75*4)} t_1 = 4 - \frac{-0.7947}{1.4992} = 4.5300 tramo = |4.5300-4| = 0.53

iteración 2 t1=4.53

t_{2} = 4.53 - \frac{cos(4.53) - sin(0.75*4.53)}{- sin(4.53) - 0.75 cos(0.75*4.53)} t_2 = 4.53 - \frac{-0.0717}{1.7089} = 4.4880 tramo= |4.4880-4.53| = 0.042

iteración 3 t2=4.4880

t_{3} = 4.4880 - \frac{cos(4.4880) - sin(0.75*4.4880)}{ - sin(4.4880) - 0.75 cos(0.75*4.4880)} t_3 = 4.4880 - \frac{-0.0000179}{1.7061} = 4.4879 tramo= |4.4879-4.4880| = 0.0001

c) Realice el análisis de la convergencia del método.

El error disminuye, el método converge. La raiz se encuentra en t=4.487989505154422

d) Con el resultado de t anterior, determine el valor de la constante k para la expresión de y2(t) que asegura el impacto contra el drone.

y_1(t) = y_2(t) sin(2 t) =kt k = \frac{sin(2 t)}{t} = \frac{sin(2*4.487989505154422)}{4.487989505154422} = 0.096714

Desarrollo con Algoritmo

Resultados

['xi', 'xnuevo', 'tramo']
[[4.0000e+00 4.5301e+00 5.3009e-01]
 [4.5301e+00 4.4880e+00 4.2071e-02]
 [4.4880e+00 4.4880e+00 3.0277e-05]]
raiz en:  4.487989505154422
con error de:  3.0276981949128867e-05

Algoritmo en Python

# Método de Newton-Raphson
import numpy as np

# INGRESO
fx  = lambda t: np.cos(1*t) - np.sin(0.75*t)
dfx = lambda t: -np.sin(t) - np.cos(0.75*t)*0.75

x0 = 4
tolera = 0.0001

# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo  = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

# convierte la lista a un arreglo.
tabla = np.array(tabla)
n = len(tabla)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)

s1Eva_2021PAOII_T3 Nutrientes en asociación de cultivos

Ejercicio: 1Eva_2021PAOII_T3 Nutrientes en asociación de cultivos

a. Realice el planteamiento del sistema de ecuaciones, presente en la forma Ax=B.

Las variables del ejercicio corresponden al consumo de nutrientes:

a = consumo de nutrientes planta de plátano
b = consumo de nutrientes planta de café
c = consumo de nutrientes planta de cacao

que generan las siguientes expresiones de consumo segun la tabla del ejercicio:

\begin{cases} 40a+110b+310c = 750 \\ 400a+15b+25c = 445 \\ 200a+560b+310c = 10 \end{cases}

que en forma A.x=B:

\begin{bmatrix} 40 && 110 && 310 \\ 400 && 15 && 25 \\ 200 && 560 && 310 \end{bmatrix}.\begin{bmatrix} a \\ b \\ c \end{bmatrix}=\begin{bmatrix} 750 \\ 445 \\ 10 \end{bmatrix}

b. De ser necesario, realice operaciones con la matriz aumentada para mejorar la convergencia con un método iterativo.

\begin{bmatrix} 40 && 110 && 310 && 750\\ 400 && 15 && 25 && 445\\ 200 && 560 && 310 &&10\end{bmatrix}

aplicando pivoteo parcial por filas

\begin{bmatrix} 400 && 15 && 25 && 445 \\ 200 && 560 && 310 &&10\\ 40 && 110 && 310 && 750\end{bmatrix}

c. En el contexto del problema, proponga un vector inicial y tolerancia.

el vector inicial no se usa con ceros, un suelo sin nutrientes no es útil para un cultivo, el vector de unos puede considerarse una opción. Otra opción es consultar con los agricultores para disponer de valores iniciales.

[1,1,1]

la toreleancia podría ser de 0.001 si consideramos que los nutrientes estan dados en kilogramos, la tolerancia sería en gramos.

d. Realice 3 iteraciones con el método de Gauss-Seidel y estime el error (papel y lápiz)

las ecuaciones para el método iterativo serán

\begin{bmatrix} 400 && 15 && 25 && 445 \\ 200 && 560 && 310 &&10\\ 40 && 110 && 310 && 750\end{bmatrix} a = \frac{445- 15 b - 25 c}{400} b = \frac{10- 200 a - 310 c}{560} c = \frac{750- 40 a -110 c}{400}

iteración 1

a = \frac{445- 15 (1) - 25 (1)}{400} = 1.0125 b = \frac{10- 200 (1.0125) - 310 (1)}{560}= -0.8973 c = \frac{750- 40 (1.0125) -110 (-0.8973)}{400} = 2.6071 error_a = |1.0125-1| = 0.0125 error_b = |-0.89736-1| = 1.8973 error_c =|2.6071-1| = 1.6071

error_iteración = max| [0.0125, 1.8973, 1.6071]| =1.8973

iteración 2 y 3 se deja como tarea:

Resultados:

respuesta X: 
[[ 0.99999468]
 [-1.99993467]
 [ 2.9999775 ]]
verificar A.X=B: 
[[444.99829063]
 [ 10.02854772]
 [750.        ]]
>>> 

e. Describa y justifique su observación sobre la convergencia del método y estime una descripción de los resultados.
El resultado muestra que el consumo de la planta b (café) es negativo, puede interpretarse como un error de datos por revisar o se iterpreta como si la planta aporta nutrientes al suelo.
Instrucciones en Python

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

import numpy as np

# INGRESO
A = np.array([[40,110,310],
              [400,15,25],
              [200,560,310]],dtype=float)
B = np.array([[750],[445],[10]],dtype=float)

X0  = np.array([1.,1.,1.],dtype=float)

tolera = 0.001
iteramax = 50

# PROCEDIMIENTO

# Matriz aumentada
AB  = np.concatenate((A,B),axis=1)
AB0 = np.copy(AB)

# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

# Para cada fila en AB
for i in range(0,n-1,1):
    # columna desde diagonal i en adelante
    columna = abs(AB[i:,i])
    dondemax = np.argmax(columna)
    
    # dondemax no está en diagonal
    if (dondemax !=0):
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal

# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]

A = np.copy(AB[:n,:m])
B = np.copy(AB[:,m])

#  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
        print(diferencia,X)
    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)

s1Eva_2021PAOII_T2 Intersección de funciones – Obstrucción Radioenlace

Ejercicio: 1Eva_2021PAOII_T2 Intersección de funciones – Obstrucción Radioenlace

Planteamiento del ejercicio

Para analizar la obstrucción se usan las expresiones del polinomio encontrado en el tema anterior y la función descrita en el enunciado. Se debe encontrar las distrancias en las que ambas expresiones son iguales, o también:

p(d)- f(d) = 0

a. Establezca un intervalo de análisis para cada raíz.

Segun la gráfica presentada en el Tema 1 se tendrán que encontrar dos raíces:

[0,700]  y [700,1300]

La expresión a usar en el ejercicio se obtiene como:

p_{3}(d) = 85.0+ 0.05941 d - 0.0001015 d^2 +3.842x10^{-8}d^3

y para el radio de fresnel:

f(d) = h_{antena} - r r = \sqrt{\frac{n \lambda d_1 d_2}{d_1+d_2}} d_2 = d_{enlace} - d_1

reemplazando los valores en las expresiones

f(d) = 100 - \sqrt{\frac{1 (0.3278 d_1 (3700-d_1}{3700}}

La expresión para el ejercicio en el intervalo de obstrucción según la gráfica del tema 1 es:

g(d) = p(d)- f(d) = 0 g(d) = 85.0+ 0.05941 d - 0.0001015 d^2 +3.842x10^{-8}d^3 - - 100 - \sqrt{\frac{1 (0.3278 d_1 (3700-d_1}{3700}}

Dado que es un ejercicio que se utiliza muchas veces en el análisis de radioenlaces, se utilizaría un algoritmo para aplicarlo en diferentes situaciones.


b. Realice al menos 3 iteraciones con el método de la Bisección para encontrar la primera raíz (izquierda)

a = 0, b = 700

con a = 0

g(0) = p(0)- f(0) = -15 p(0) = 85.0+ 0.05941 (0) - 0.0001015 (0)^2 +3.842x10^{-8}(0)^3 = 85 f(d) = 100 - \sqrt{\frac{1 (0.3278 d_1 (3700-d_1}{3700}} =100

con b = 700

g(700) = p(700)- f(700) = 3.67 p(700) = 85.0+ 0.05941 (700) - 0.0001015 (700)^2 +3.842x10^{-8}(700)^3 f(700) = 100 - \sqrt{\frac{1 (0.3278 (700) (3700-(700)}{3700}}

como hay cambio de signo entre g(a) y g(b), debe existir una raiz en el intervalo

iteración  1

c= \frac{a+b}{2} = \frac{700-0}{2} = 350 p(350) = 85.0+ 0.05941 (350) - 0.0001015 (350)^2 +3.842x10^{-8}(350)^3 f(350) = 100 - \sqrt{\frac{1 (0.3278 (350) (3700-(350}{3700}} g(350) = p(350)- f(350) = 5.199

como la referencia es

g(0) = -15 g(700) = 3.67

hay cambio de signo entre [0,350]

error estimado = 350-0 = 350

iteración  2

c= \frac{a+b}{2} = \frac{350-0}{2} = 175 p(175) = 85.0+ 0.05941 (175) - 0.0001015 (175)^2 +3.842x10^{-8}(175)^3 f(175) = 100 - \sqrt{\frac{1 (0.3278 (175) (3700-(175}{3700}} g(175) = -113.1

como la referencia es

g(0) = -15 g(350) = 5.199

hay cambio de signo entre [175,350]

error estimado = 350-175 = 175

iteración  3

c= \frac{a+b}{2} = \frac{175-350}{2} = 262.5 p(262.5) = 85.0+ 0.05941 (262.5) - 0.0001015 (262.5)^2 +3.842x10^{-8}(262.5)^3 f(262.5) = 100 - \sqrt{\frac{1 (0.3278 (262.5) (3700-(262.5}{3700}} g(262.5) = 3.23

Instrucciones en Python

resultado usando el algoritmo:

[ i, a, c, b, f(a), f(c), f(b), tramo]
1 0.000 350.000 700.000 -15.000 5.199 3.670 700.000 
2 0.000 175.000 350.000 -15.000 -0.113 5.199 350.000 
3 175.000 262.500 350.000 -0.113 3.237 5.199 175.000 
4 175.000 218.750 262.500 -0.113 1.755 3.237 87.500 
5 175.000 196.875 218.750 -0.113 0.872 1.755 43.750 
6 175.000 185.938 196.875 -0.113 0.393 0.872 21.875 
7 175.000 180.469 185.938 -0.113 0.143 0.393 10.938 
8 175.000 177.734 180.469 -0.113 0.016 0.143 5.469 
9 175.000 176.367 177.734 -0.113 -0.048 0.016 2.734 
10 176.367 177.051 177.734 -0.048 -0.016 0.016 1.367 
11 177.051 177.393 177.734 -0.016 -0.000 0.016 0.684 
12 177.393 177.563 177.734 -0.000 0.008 0.016 0.342 
13 177.393 177.478 177.563 -0.000 0.004 0.008 0.171 
14 177.393 177.435 177.478 -0.000 0.002 0.004 0.085 
15 177.393 177.414 177.435 -0.000 0.001 0.002 0.043 
16 177.393 177.403 177.414 -0.000 0.000 0.001 0.021 
17 177.393 177.398 177.403 -0.000 0.000 0.000 0.011 
18 177.393 177.395 177.398 -0.000 -0.000 0.000 0.005 
19 177.395 177.397 177.398 -0.000 0.000 0.000 0.003 
20 177.395 177.396 177.397 -0.000 0.000 0.000 0.001 
21 177.395 177.396 177.396 -0.000 0.000 0.000 0.001 
raiz:  177.39558219909668
       raiz en:  177.39558219909668
error en tramo:  0.000667572021484375
# Algoritmo de Bisección
# [a,b] se escogen de la gráfica de la función
# error = tolera

import numpy as np
import matplotlib.pyplot as plt

# INGRESO
p = lambda d: 85.0 + 0.05941*d - 0.0001015*d**2 +3.842e-8*d**3
f = lambda d: 100 - np.sqrt(1*(0.3278)*d*(3700-d)/3700) 
fx = lambda d: p(d)- f(d)
a_todo = 0
b_todo = 1300
a = 0
b = 700
tolera = 0.001

muestras = 41

# PROCEDIMIENTO
# 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
np.set_printoptions(precision = 4)
print('[ i, a, c, b, f(a), f(c), f(b), tramo]')
# 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)
print('       raiz en: ', c)
print('error en tramo: ', tramo)


# GRAFICA
xi = np.linspace(a_todo,b_todo,muestras)
yi = fx(xi)
pi = p(xi)
ri = f(xi)
plt.plot(xi,yi,label='g(d)=p(d)-f(d)')
plt.plot(xi,pi,label='p(d)')
plt.plot(xi,ri,label='f(d)_Fresnel')
plt.axhline(0, color='grey')
plt.xlabel('d')
plt.legend()
plt.title('obstruccion Fresnel')
plt.show()

c. Desarrolle al menos 3 iteraciones con el método del Punto fijo para encontrar el segundo punto (derecha)

La siguiente raíz se encuentra con algoritmo presentado o también con el algoritmo del siguiente literal usando la misma expresion g(d).

a = 700 b= 1300

 raiz en: 867.1000480651855
error en tramo: 0.00057220458984375

s1Eva_2021PAOII_T1 Interpolación para perfil de terreno

Ejercicio: 1Eva_2021PAOII_T1 Interpolación para perfil de terreno

a. Plantee y desarrolle un polinomio P3(d1) de grado 3, que describa el perfil del terreno en el intervalo [0,1300] de distancias a la primera antena d1.

Para plantear el ejercicio, se observan los punto en el intervalo [0,1300], asi como los tamaños de paso.

Como el polinomio es de grado 3, solo se requien 4 tramos o 4 muestras, el intervalo presenta 5. Para incluir el punto de inicio, fin y máximo, las muestras seleccionadas son las que se muestran en la tabla siguiente:

Δd1 distancia d1 Perfil de Terreno DifDiv1 DifDiv2 DifDiv3
350 0 85 DifDiv11 = 0.02857 DifDiv12 = -6.122E-05 DifDiv13 =  -1.123E-05
350 350 95 DifDiv21 = -0.01428 DifDiv22 = -1.127E-05
600 700 90 DifDiv31 =
-0.025
1300 75

Los tamaños de paso, Δd1, son diferentes, por lo que se pueden usar el método de diferencias divididas de Newton o el Método de Lagrange.

Usando deiferencias divididas de Newton, se completa la tabla con las expresiones:

DifDiv_{11} = \frac{95-85}{350-0} = 0.02857143 DifDiv_{21} = \frac{90-95}{700-350} = -0.01428 DifDiv_{31} = \frac{75-90}{1300-700} =-0.025 DifDiv_{12} = \frac{-0.01428-0.02857}{700-0} = -6.122E-05 DifDiv_{22} = \frac{-0.025-(-0.01428)}{1300-350} = -1.127E-05 DifDiv_{13} = \frac{(-1.127E-05) -(6.122E-05)}{1300-0} = 1.123E-05

con lo que se puede constuir el polinomio

p_{3}(d) = 85 + 0.02857(d-0)+(-6.122E5)(d-0)(d-350)+ -(1.123E-5)(d-0)(d-350)(d-700)

simplificando el polinomio se obtiene:

p_{3}(d) = 85.0+ 0.05941 d - 0.0001015 d^2 +(3.842E-8)d^3

b. Calcule el error sobre el o los datos que no se usaron en el intervalo

El valor no usado que estaba en el intervalo es d=1000, que se evalúa en p3(d) y se compara con el valor muestreado para ver el error del modelo

p_{3}(1000) = 85.0+ 0.05941 (1000) - 0.0001015 (1000)^2 +(3.842E-8)(1000)^3 p_{3}(1000) = 81.26 Error = |81.26 - 80| = 1.26

encontrando que el error es de 1.26m de altitud del terreno a los 1000m de distancia desde la referencia a la izquierda.

c. Desarrolle y justifique una propuesta para disminuir los errores encontrados en el literal anterior, sobre el mismo intervalo, es decir obtiene un nuevo polinomio (use algoritmo).

En caso de requerir una mayor precisión, se puede aumentar el grado del polinomio al incluir el punto d=1000 no usado en el paso anterior.

Usando el algoritmo con Python se obtiene:

Tabla Diferencia Dividida
[['i   ', 'xi  ', 'fi  ', 'F[1]', 'F[2]', 'F[3]', 'F[4]', 'F[5]']]
[[ 0.0000e+00  0.0000e+00  8.5000e+01  2.8571e-02 -6.1224e-05  3.1920e-08  2.1666e-11  0.0000e+00]
 [ 1.0000e+00  3.5000e+02  9.5000e+01 -1.4286e-02 -2.9304e-05  6.0086e-08  0.0000e+00  0.0000e+00]
 [ 2.0000e+00  7.0000e+02  9.0000e+01 -3.3333e-02  2.7778e-05  0.0000e+00  0.0000e+00  0.0000e+00]
 [ 3.0000e+00  1.0000e+03  8.0000e+01 -1.6667e-02  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00]
 [ 4.0000e+00  1.3000e+03  7.5000e+01  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00  0.0000e+00]]
dDividida: 
[ 2.8571e-02 -6.1224e-05  3.1920e-08  2.1666e-11  0.0000e+00]
polinomio: 
2.16658863275405e-11*x*(x - 1000.0)*(x - 700.0)*(x - 350.0) + 3.19204604918891e-8*x*(x - 700.0)*(x - 350.0) - 6.12244897959184e-5*x*(x - 350.0) + 0.0285714285714286*x + 85.0
polinomio simplificado: 
2.16658863275405e-11*x**4 - 1.24946064795689e-8*x**3 - 6.6683650518237e-5*x**2 + 0.0525123706702654*x + 85.0

junto a la gráfica:

d. Escriba sus conclusiones y recomendaciones sobre los resultados obtenidos entre los dos polinomios.

Dado que el error se considera mínimo en el intervalo, y el coeficiente del polinomio para x4 es del orden de 10-11 se recomendaría usar el polinomio de grado 3.


Instrucciones en Python

# 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([0.0,350,700,1000,1300])
fi = np.array([85.0,95,90,80,75])


# 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()

s1Eva_2021PAOI_T2 Atención hospitalaria con medicamentos limitados

Ejercicio: 1Eva_2021PAOI_T2 Atención hospitalaria con medicamentos limitados

literal a

Se plantea el sistema de ecuaciones de acuerdo a la cantidad de medicamentos que se administra a cada tipo de paciente, siendo las variables X=[a,b,c,d] de acuerdo a la tabla:

Niños Adolescentes Adultos Adultos Mayores Medicamentos /semana
Medicamento_A 0.3 0.4 1.1 4.7 3500
Medicamento_B 1 3.9 0.15 0.25 3450
Medicamento_C 0 2.1 5.6 1.0 6500
0.3 a + 0.4 b + 1.1 c + 4.7 d = 3500 a + 3.9 b + 0.15 c + 0.25 d = 3450 0 a +2.1 b + 5.6 c + 1.0 d = 6500

Mostrando que el número de incógnitas no es igual al número de ecuaciones, por lo que sería necesario de disponer de otra ecuación, o  usar una variable libre.


literal b

Siguiendo las indicaciones sobre la variable libre que sea para el grupo de pacientes niños, te tiene que

0.4 b + 1.1 c + 4.7 d = 3500 - 0.3 K 3.9 b + 0.15 c + 0.25 d = 3450 - K 2.1 b + 5.6 c + 1.0 d = 6500 - 0 K

y haciendo K = 100, se convierte en :

0.4 b + 1.1 c + 4.7 d = 3500 - 0.3(100) = 3470 3.9 b + 0.15 c + 0.25 d = 3450 - 100 = 3350 2.1 b + 5.6 c + 1.0 d = 6500 - 0 = 6500

literal c

Antes de desarrollar el ejercicio para el algoritmo se requiere convertir el sistema de ecuaciones a su forma matricial. Por lo que se plantea la matriz aumentada:

\begin{pmatrix} 0.4 & 1.1 & 4.7 & \Big| & 3470 \\ 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500\end{pmatrix}

Se realiza el pivoteo parcial por filas, que al aplicar en la primera columa, se intercambia la primera y segunda fila, de tal forma que el mayor valor quede en la primera casilla de la diagonal.

\begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 0.4 & 1.1 & 4.7 & \Big| & 3470 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500\end{pmatrix}

se continua el mismo proceso para la segunda casilla de la diagonal hacia abajo, se aplica también pivoteo parcial,

\begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3350 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500 \\ 0.4 & 1.1 & 4.7 & \Big| & 3470 \end{pmatrix}

Con lo que el sistema queda listo para resolver por cualquier método.

Si seleccionamos Gauss-Seidel, el vector inicial de acuerdo al enunciado será X0=[100,100,100]

y las ecuaciones a resolver son:

b = \frac{3350 - 0.15 c - 0.25 d}{3.9} c = \frac{6500 - 2.1 b - 1.0 d}{5.6} d = \frac{3470- 0.4 b - 1.1 c}{4.7}

iteración 1

X0=[100,100,100]

b = \frac{3350 - 0.15(100) - 0.25 (100)}{3.9} = 848.71 c = \frac{6500 - 2.1 (848.71 ) - 1.0 (100)}{5.6} = 824.58 d = \frac{3470- 0.4 (848.71 ) - 1.1 (824.58)}{4.7} = 473.07

X1=[848.71, 824.58, 473.07]

diferencia = [848.71, 824.58, 473.07] – [100,100,100]

diferencia = [748.71, 724.58, 373.07]

error = max|diferencia| = 748.71

resultado del error es mayor que tolerancia de 0.01, se continua con la siguiente iteración.

iteración 2

X1=[848.71, 824.58, 473.07]

b = \frac{3350 - 0.15 (824.58) - 0.25 (473.07)}{3.9} = 796.93 c = \frac{6500 - 2.1 (796.93) - 1.0 (473.07)}{5.6} = 777.38 d = \frac{3470- 0.4 (796.93) - 1.1 (777.38)}{4.7} = 488.53

X2 = [796.93, 777.38, 488.53]

diferencia = [796.93, 777.38, 488.53]  – [848.71, 824.58, 473.07]

diferencia = [-51.78 ,  -47.2, 15.52]

error = max|diferencia| = 51.78

iteración 3

X2 = [796.93, 777.38, 488.53]

b = \frac{3350 - 0.15 (777.38) - 0.25 (488.53)}{3.9} = 797.75 c = \frac{6500 - 2.1 (797.75) - 1.0 (488.53)}{5.6} = 774.31 d = \frac{3470- 0.4 (797.75) - 1.1 (774.31)}{4.7} = 489.18

x3 = [ 797.75, 774.31, 489.18]

diferencias = [ 797.75, 774.31, 489.18] – [796.93, 777.38, 488.53]

diferencias = [0.82, -3.07, 0.65]

error = max|diferencia| = 3.07

Observación, el error disminuye en cada iteración, por lo que el sistema converge.

solución con el algoritmo, solo se toma la parte entera de la respuesta, pues los pacientes son números enteros.

respuesta X: [797.83, 774.16, 489.20]

con lo que el primer resultado es:

 X = [797, 774, 489]

literal d

Si la cantidad de pacientes presentados en una seman es la que se indica en el enunciado [350,1400,1500,1040], se determina que el sistema hospitalario estatal no podrá atender a todos los pacientes al compara con la capacidad encontrada en el literal c.

Hay un exceso de 2129 pacientes encontrados por grupo:

exceso = [350, 1400, 1500, 1040] – [100, 797, 774, 489] = [250, 603, 726, 551]

con lo que se recomienda una medida de confinamiento para disminuir los contagios y poder atender a los pacientes que se presenten.

literal e

Siendo K = 0, al estar vacunados todos los niños, y suponiendo que la vacuna es una cura, se tiene:

0.4 b + 1.1 c + 4.7 d = 3500 - 0.3 K = 3500 3.9 b + 0.15 c + 0.25 d = 3450 - K = 3450 2.1 b + 5.6 c + 1.0 d = 6500 - 0 k = 6500

pivoteado parcialmente por filas se tiene:

\begin{pmatrix} 3.9 & 0.15 &0.25 & \Big| & 3450 \\ 2.1 & 5.6 & 1.0 &\Big| & 6500 \\ 0.4 & 1.1 & 4.7 & \Big| & 3500 \end{pmatrix}

Resolviendo por un método directo:

se obtiene:

respuesta X: [823.46, 763.35, 495.94]

que al compararse con la capacidad anterior en números enteros se encuentra una diferencia de un incremento de 21 pacientes neto ante la condición de usar todos los medicamentos disponibles.

diferencia = [823, 763, 495] – [797, 774, 489] = [26, -11, 6]


 Instrucciones en Python

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

import numpy as np

# INGRESO
A = np.array([[0.4, 1.10, 4.7 ],
              [3.9, 0.15, 0.25],
              [2.1, 5.60, 1.0  ]])
k = 100
B = np.array([3500.-0.3*k,3450-1*k,6500-0*k])

X0  = np.array([100.0,100,100])

tolera = 0.001
iteramax = 100

# PROCEDIMIENTO
# Matriz aumentada
Bcolumna = np.transpose([B])
AB  = np.concatenate((A,Bcolumna),axis=1)
AB0 = np.copy(AB)

# Pivoteo parcial por filas
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

# Para cada fila en AB
for i in range(0,n-1,1):
    # columna desde diagonal i en adelante
    columna = abs(AB[i:,i])
    dondemax = np.argmax(columna)
    
    # dondemax no está en diagonal
    if (dondemax !=0):
        # intercambia filas
        temporal = np.copy(AB[i,:])
        AB[i,:]  = AB[dondemax+i,:]
        AB[dondemax+i,:] = temporal

A = np.copy(AB[:,:n])
B = np.copy(AB[:,n])

# 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
    print(X)
    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)

s1Eva_2021PAOI_T1 Función recursiva y raíces de ecuaciones

Ejercicio: 1Eva_2021PAOI_T1 Función recursiva y raíces de ecuaciones

Literal a

Evaluando las sucesión de la forma recursiva:

xi
[-0.45       -0.4383     -0.4458     -0.441      -0.4441 
 -0.4421     -0.4434     -0.4425     -0.4431     -0.4427
 -0.4429     -0.4428     -0.4429     -0.4428     -0.4429]
errores
[ 1.1745e-02 -7.5489e-03  4.8454e-03 -3.1127e-03  1.9986e-03
 -1.2837e-03  8.2430e-04 -5.2939e-04  3.3996e-04 -2.1833e-04
  1.4021e-04 -9.0044e-05  5.7826e-05 -3.7136e-05  0.0000e+00]

literal b

Se puede afirmar que converge, observe la diferencia entre cada dos valores consecutivos de la sucesión … (continuar de ser necesario)

literal c

Para el algoritmo se requiere la función f(x) y su derivada f'(x)

f(x) = x +ln(x+2) f'(x) = 1 + \frac{1}{x+2}

x0 = -0.45

itera = 1

x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)} x_{1} = x_0 -\frac{f(x_0)}{f'(x_0)} = -0.45 -\frac{-0.45+ln(-0.45+2)}{1 + \frac{1}{-0.45+2}}

x1 = -0.44286

error = |x1-x0| = |-0.44286 -(-0.45)| = 0.007139

itera = 2

x_{2} = -0.4428 -\frac{-0.4428 + ln(-0.45+2)}{1 + \frac{1}{-0.4428+2}}

x1 = -0.44286

error = |x1-x0| = |-0.44285 -(-4.4286)| = 6.4394e-06

con lo que se cumple el valor de tolerancia y no se requiere otra iteración

la raiz se encuentra en x = -0.44286

Solución con algoritmo

xi
[-0.45       -0.4383     -0.4458     -0.441      -0.4441 
 -0.4421     -0.4434     -0.4425     -0.4431     -0.4427
 -0.4429     -0.4428     -0.4429     -0.4428     -0.4429]
errores
[ 1.1745e-02 -7.5489e-03  4.8454e-03 -3.1127e-03  1.9986e-03
 -1.2837e-03  8.2430e-04 -5.2939e-04  3.3996e-04 -2.1833e-04
  1.4021e-04 -9.0044e-05  5.7826e-05 -3.7136e-05  0.0000e+00]
['xi', 'xnuevo', 'tramo']
[[-4.5000e-01 -4.4286e-01  7.1392e-03]
 [-4.4286e-01 -4.4285e-01  6.4394e-06]]
raiz en:  -0.44285440100759543
con error de:  6.439362322474551e-06
>>> 

Instrucciones en Python

import numpy as np
import matplotlib.pyplot as plt

# literal a sucesión en forma recursiva
def secuenciaL(n,x0):
    if n == 0:
        xn = x0
    if n>0:
        xn = np.log(1/(2+secuenciaL(n-1,x0)))
    return(xn)
x0 = -0.45
n = 15
xi = np.zeros(n,dtype=float)
xi[0] = x0
errado = np.zeros(n,dtype=float)
for i in range(1,n,1):
    xi[i] = secuenciaL(i,x0)
    errado[i-1] = xi[i] - xi[i-1]
   
np.set_printoptions(precision=4)
print('xi: ')
print(xi)
print('errado: ')
print(errado)

#Grafica literal a y b
plt.plot(xi)
plt.plot(xi,'o')
plt.xlabel('n')
plt.ylabel('xn')
plt.show()

# Método de Newton-Raphson
# Ejemplo 1 (Burden ejemplo 1 p.51/pdf.61)

import numpy as np

# INGRESO
fx  = lambda x: x+np.log(x+2)
dfx = lambda x: 1+ 1/(x+2)

x0 = -0.45
tolera = 0.0001

a = -0.5
b = -0.2
muestras = 21
# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo  = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

# convierte la lista a un arreglo.
tabla = np.array(tabla)
n = len(tabla)

# para la gráfica
xj = np.linspace(a,b,muestras)
fj = fx(xj)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)

plt.plot(xj,fj)
plt.axhline(0, color='grey')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()

litera d: tarea

s1Eva_2021PAOI_T3 Interpolar, modelo de contagios 2020

Ejercicio: 1Eva_2021PAOI_T3 Interpolar, modelo de contagios 2020

literal a

Los datos de los pacientes casos graves entre las semanas 11 a la 20, que son el intervalo donde será válido el polinomio de interpolación son:

semana 11 12 13 14 15 16 17 18 19 20
casos graves 1503 3728 7154 6344 4417 3439 2791 2576 2290 2123

de los cuales solo se usarán los indicados en el literal a : 11,13,16,18,20.

xi0 = [    9,   10,   11,   12,   13,   14,
          15,   16,   17,   18,   19,   20,
          21,   22,   23,   24,   25,   26 ])
fi0 = [ 1435, 1645, 1503, 3728, 7154, 6344,
        4417, 3439, 2791, 2576, 2290, 2123,
        2023, 2067, 2163, 2120, 2125, 2224 ])
xi = [  11,   13,   16,   18,   20]
fi = [1503, 7154, 3439, 2576, 2123]

Se observa que los datos estan ordenados en forma ascendente respecto a la variable independiente, tambien se determina que no se encuentran equidistantes entre si (13-11=2, 16-13=3). Por lo que se descarta usar el método de diferencias finitas avanzadas.

Los métodos que se podrían usar con puntos no equidistantes en el eje semanas serían el método de diferencias divididas de Newton o el  método de Lagrange.

Seleccionando por ejemplo, Diferencias divididas de Newton, donde primero se realiza la tabla:

xi fi f[x1,x0] f[x2,x1,x0] f[x3,x2,x1,x0] f[x4,x3,x2,x1,x0]
11 1503 =(7154-1503) /(13-11) = 2825.5 =(-1238.33-2835.5) /(16-11) = -812.76 =(161.36-(-812.76)) /(18-11) = 139.16 =(-15.73-139.16) /(20-11) = -17.21
13 7154 (3439-7154) /(16-13) = -1238.33 (-431.5-(-1238.33)) /(18-13) = 161.36 (51.25-161.36) /(20-13)= -15.73 —-
16 3439 (2576-3439) /(18-16) = -431.5 (-226.5-(-431.5)) /(20-16) = 51.25 —-
18 2576 (2123-2576) /(20-18) = -226.5 —-
20 2123 —-

con lo que se puede contruir el polinomio usando las diferencias divididas para el intervalo dado:

[2825.5    -812.76  139.16  -17.21]
p_4(x) = 1503 + 2825.5(x-11) - 812.76(x - 13)(x - 11) + 139.16(x - 16)(x - 13)(x - 11) - 17.21(x - 18)(x - 16) (x - 13) (x - 11)

Simplificando el algoritmo se tiene:

p_4(x) = - 1172995.28 + 298304.50 x - 27840.50x^2 + 1137.36x^3 - 17.21 x^4


literal b

El cálculo de los errores se puede realizar usando el polinomio de grado 4 encontrado, notando que los errores deberían ser cero para los puntos usados para el modelo del polinomio.

xi fi p4(x) |error|
11 1503 1503 0
12 3728 6110.96 2382.96
13 7154 7154 0
14 6344 6293.18 50.81
15 4417 4776.52 359.52
16 3439 3439 0
17 2791 2702.53 88.46
18 2576 2576 0
19 2290 2655.22 365.22
20 2123 2123 0

literal c

Podría aplicarse uno de varios criterios, lo importante por lo limitado del tiempo en la evaluación son las conclusiones y recomendaciones expresadas en el literal e, basadas en lo realizado en los literales c y d. Teniendo como opciones:

– cambiar uno de los puntos selecionados, mateniendo así el grado del polinomio
– aumentar el número de puntos usados para armar el polinomio con grado mayor
– dividir el intervalo en uno o mas segmentos, con el correspondiente número de polinomios.

Se desarrolla la opción de cambiar uno de los puntos seleccionados, usando para esta ocasión como repaso la interpolación de Lagrange. Para los puntos se usa el punto con mayor error de la tabla del literal anterior y se elimina el punto penúltimo, es decir se usa la semana 12 en lugar de la semana 18 de la siguiente forma:

xi = [  11,   12,   13,   16,   20]
fi = [1503, 3728, 7154, 3439, 2123]
p_4(x) = 1503 \frac{(x-12)(x-13)(x-16)(x-20)}{(11-12)(11-13)(11-16)(11-20)} + 3728\frac{(x-11)(x-13)(x-16)(x-20)}{(12-11)(12-13)(12-16)(12-20)} + 7154\frac{(x-11)(x-12)(x-16)(x-20)}{(13-11)(13-12)(13-16)(13-20)} + 3439\frac{(x-11)(x-12)(x-13)(x-20)}{(16-11)(16-12)(16-13)(16-20)} + 2123\frac{(x-11)(x-12)(x-13)(x-16)}{(20-11)(20-12)(20-13)(20-16)}

Simplificando el polinomio:

p_4(x) = \frac{46927445}{21} - \frac{1655552687}{2520} x + \frac{715457663}{10080}x^2 - \frac{8393347}{2520} x^3 + \frac{577153}{10080} x^4

literal d

El cálculo de los errores se puede realizar usando el polinomio de grado 4 encontrado, notando que los errores deberían ser cero para los puntos usados para el modelo del polinomio.

xi fi p4(x) error
11 1503 1503 0
12 3728 3728 0
13 7154 7154 0
14 6344 8974.01 2630.01
15 4417 7755.22 3338.22
16 3439 3439 0
17 2791 -2659.13 -5450.13
18 2576 -7849.45 -10425.45
19 2290 -8068.1 -10358.1
20 2123 2123 0

literal e

El cambio aplicado a los puntos usados en el modelo del polinomio disminuyó el error entre las semanas 11 a 13. Sin embargo la magnitud del error aumentó  para las semanas posteriores a la 13, es decir aumentó la distorsión de la estimación y se recomienda realizar otras pruebas para mejorar el modelo aplicando los otros criterios para determinar el que tenga mejor desempeño respecto a la medida de error.



Intrucciones en Python

Literal a y b. Desarrollado a partir del algoritmo desarrollado en clases:

# 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
xi0 = np.array([    9,   10,   11,   12,   13,   14,
                   15,   16,   17,   18,   19,   20,
                   21,   22,   23,   24,   25,   26 ])
fi0 = np.array([ 1435, 1645, 1503, 3728, 7154, 6344,
                 4417, 3439, 2791, 2576, 2290, 2123,
                 2023, 2067, 2163, 2120, 2125, 2224 ])

xi1 = np.array([   11,   12,   13,   14,   15,   16,
                   17,   18,   19,   20 ])
fi1 = np.array([ 1503, 3728, 7154, 6344, 4417, 3439,
                 2791, 2576, 2290, 2123 ])

xi = np.array([  11,   13,   16,   18,   20])
fi = np.array([1503, 7154, 3439, 2576, 2123])

# 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)

# calcula errores en intervalo usado
pfi1 = px(xi1)
errado1 = np.abs(fi1-pfi1)

# 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)
print('errores en intervalo:')
print(xi1)
print(errado1)

# Gráfica
plt.plot(xi0,fi0,'o', label = 'Puntos')
plt.plot(xi,fi,'ro', label = 'Puntos')
for i in range(0,n,1):
    etiqueta = '('+str(xi[i])+','+str(fi[i])+')'
    plt.annotate(etiqueta,(xi[i],fi[i]))
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencias Divididas - Newton')
plt.grid()
plt.show()

Literal c y d. Se puede continuar con el algoritmo anterior. Como repaso se adjunta un método diferente al anterior.

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

# INGRESO , Datos de prueba
xi0 = np.array([    9,   10,   11,   12,   13,   14,
                   15,   16,   17,   18,   19,   20,
                   21,   22,   23,   24,   25,   26 ])
fi0 = np.array([ 1435, 1645, 1503, 3728, 7154, 6344,
                 4417, 3439, 2791, 2576, 2290, 2123,
                 2023, 2067, 2163, 2120, 2125, 2224 ])

xi2 = np.array([   11,   12,   13,   14,   15,   16,
                   17,   18,   19,   20 ])
fi2 = np.array([ 1503, 3728, 7154, 6344, 4417, 3439,
                 2791, 2576, 2290, 2123 ])

xi = np.array([  11,   12,   13,   16,   20])
fi = np.array([1503, 3728, 7154, 3439, 2123])

# 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)

# calcula errores en intervalo usado
pfi2 = px(xi2)
errado2 = np.abs(fi2-pfi2)

# 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)
print('errores en intervalo:')
print(xi2)
print(errado2)

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