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

s3Eva_2020PAOI_T2 Modelo epidemiológico no letal

Ejercicio: 3Eva_2020PAOI_T2 Modelo epidemiológico no letal

El ejercicio representa un sistema de ecuaciones diferenciales ordinarias, que serán resueltas usando Runge-Kutta de 2do Orden.

De compararse con la curva de contagios de Covid-19 se tienen diferencias en la población recuperada, pues el modelo se considera no letal por lo que no se contabiliza el número de fallecidos.

El módelo es el más básico y permite cambiar por ejemplo la tasa de infección, y se ve los cambios en la curva de infectados. Se puede observar lo que se indicaba como objetivo de «aplanar la curva» al disminuir la población expuesta mediante cambiar la tasa de infección al exponer más o menos población al contagio por iteacción entre «suceptibles» e «infectados.


Desarrollo analítico

Las fórmulas para el algoritmo se identifican como:
binfecta = 1.4
grecupera = 1/4

f(t,S,I,R) = -1.4(S)(I) g(t,S,I,R) = 1.4(S)(I) - \frac{1}{4}(I) w(t,S,I,R) = \frac{1}{4}(I)

que luego se usan en cada iteración que se registra en la tabla empezando con las condiciones iniciales

itera Si Ii Ri
0 1 0.001 0
1 0.9977 0.002809 0.0003937
2 0.9916 0.007862 0.0014987
3 tarea

itera = 1
K_{1S} = h. f(t_i,S_i,I_i,R_i) =

= 1(-1.4)( 1)(0.001) = -0.0014 K_{1I} = h. g(t_i,S_i,I_i,R_i) = = 1\Big((1.4)(1)(0.001) - \frac{1}{4}0.001 \Big) = 0.00115 K_{1R} = h. w(t_i,S_i,I_i,R_i) = = 1 \Big( \frac{1}{4} 0.001 \Big) = 0.00025 K_{2S} = h. f(t_i+h, S_i + K_{1S}, I_i+K_{1I}, R_i +K_{1R}) = = 1 \Big( (-1.4)(1-0.0014)(0.001+0.00115)\Big) = = -0.003005786 K_{2I} = h. g(t_i+h, S_i + K_{1S}, I_i+K_{1I}, R_i + K_{1R}) = = 1\Big((1.4)(1-0.0014)(0.001+0.00115) -\frac{1}{4} (0.001+0.00115)\Big) = = 0.002468286 K_{2R} = h w(t_i+h, S_i + K_{1S}, I_i+K_{1I}, R_i +K_{1R}) = 1\Big( \frac{1}{4}(0.001+0.00115) \Big) = 0.0005375 S_i = S_i + \frac{K_{1S}+K_{2S}}{2} = 1 + \frac{-0.0014 -0.003005786}{2} = 0.997797107 I_i = I_i + \frac{K_{1I}+K_{2I}}{2} = 0.001 + \frac{0.00115+0.002468286}{2} = 0.002809143 R_i = R_i + \frac{K_{1R}+K_{2R}}{2} = 0 + \frac{0.00025 + 0.0005375}{2} = 0.00039375 t_i = t_i + h = 1 + 1 = 2

 

itera = 2

K_{1S} = 1(-1.4)(0.9977)(0.002809) = -0.003924 K_{1I} = 1 \Big(1.4(0.9977)(0.002809) - \frac{1}{4}0.002809 \Big) = 0.003221 K_{1R} = 1(\frac{1}{4} 0.002809) = 0.0007022 K_{2S} = 1 \Big( -1.4*(0.9977-0.003924)(0.002809+0.003221) \Big) = = -0.008391 K_{2I} = 1 \Big(1.4(0.9977-0.003924)(0.002809+0.003221) - \frac{1}{4}(0.002809+0.003221) \Big) = 0.006883 K_{2R} = 1 \Big( \frac{1}{4}(0.002809+0.0032218) \Big) = 0.001507

 

S_i = 0.9977 + \frac{-0.003924 -0.008391}{2} = 0.9916

 

I_i = 0.002809 + \frac{0.003221+0.006883}{2} = 0.007862

 

R_i = 0.0003937 + \frac{0.0007022 + 0.001507}{2} = 0.001498

 

t_i = t_i + h = 2 + 1 = 3

 

itera 3 – TAREA


Instrucciones en Python

# 3Eva_2020PAOI_T2 Modelo epidemiológico no letal
# Sistemas EDO con Runge-Kutta de 2do Orden
import numpy as np

def rungekutta2_fgw(f,g,w,t0,x0,y0,z0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,4),dtype=float)
    tabla[0] = [t0,x0,y0,z0]
    ti = t0
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi,zi)
        K1y = h * g(ti,xi,yi,zi)
        K1z = h * w(ti,xi,yi,zi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y, zi +K1z)
        K2y = h * g(ti+h, xi + K1x, yi+K1y, zi +K1z)
        K2z = h * w(ti+h, xi + K1x, yi+K1y, zi +K1z)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        zi = zi + (1/2)*(K1z+K2z)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi,zi]
    tabla = np.array(tabla)
    return(tabla)

# Programa
# Parámetros de las ecuaciones

binfecta = 1.4
grecupera = 1/4
# Ecuaciones
f = lambda t,S,I,R : -binfecta*S*I
g = lambda t,S,I,R : binfecta*S*I - grecupera*I
w = lambda t,S,I,R : grecupera*I
# Condiciones iniciales
t0 = 0
S0 = 1.0
I0 = 0.001
R0 = 0.0

# parámetros del algoritmo
h = 1.0
muestras = 51

# PROCEDIMIENTO
tabla = rungekutta2_fgw(f,g,w,t0,S0,I0,R0,h,muestras)
ti = tabla[:,0]
Si = tabla[:,1]
Ii = tabla[:,2]
Ri = tabla[:,3]
# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, Si, Ii, Ri]')
print(tabla)

# Grafica tiempos vs población
import matplotlib.pyplot as plt
plt.plot(ti,Si, label='S')
plt.plot(ti,Ii, label='I')
plt.plot(ti,Ri, label='R')
plt.title('Modelo SIR')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

s1Eva_IIT2009_T1 Movimiento de partícula en plano

Ejercicio: 1Eva_IIT2009_T1 Movimiento de partícula en plano

a. Planteamiento del problema

Las ecuaciones expresan las trayectorias de dos partículas,

x(t) = 3 \sin ^{3}(t)-1 y(t) = 4 \sin (t)\cos (t)

que para que se encuentren o choquen, sus coordenadas deberían ser iguales.

x(t) = y(t) 3 \sin ^{3}(t)-1 = 4 \sin (t)\cos (t)

Se reordena la igualdad, de tal manera que uno de los lados sea cero, de la forma f(t) = 0 que es usada en el algoritmo de búsqueda de raíces.

3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t) = 0 f(t) = 3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t)

b. Intervalo de búsqueda de raíz

Como la variable independiente es tiempo, el evento a buscar se suponte sucede en tiempos positivos t>=0, por lo que el valor inicial a la izquierda del intervalo será a=0

Para el caso de b, a la derecha, se usa lo indicado en el enunciado para la pregunta dell literal b), donde se indica t ∈ [0, π/2], por lo que b = π/2

[0, π/2]

verificando que exista cambio de signo entre f(a) y f(b)

f(0) = 3 \sin ^{3}(0)-1 - 4 \sin (0)\cos (0) = -1 f(\pi /2) = 3 \sin ^{3}(\pi /2)-1 - 4 \sin (\pi /2)\cos (\pi /2) = 3 (1)^3-1 - 4 (1)(0) = 2

con lo que se comprueba que al existir cambio de signo, debe existir una raíz en el intervalo.


c. Método de Falsa Posición


Desarrollo analítico con lápiz y papel

se trata de mostrar los pasos en al menos tres iteraciones del método, usando las siguientes expresiones:

f(t) = 3 \sin ^{3}(t)-1 - 4 \sin (t)\cos (t) c = b - f(b) \frac{a-b}{f(a)-f(b)}

[0, π/2]

iteración 1

a = 0 , b = \pi/2

tomando los datos al validar los puntos extremos

f(0) = -1 f(\pi /2) = 2 c = \pi/2 - 2 \frac{0-\pi/2}{-1-2} = \pi/6 f(\pi/6) = 3 \sin ^{3}(\pi/6)-1 - 4 \sin (\pi/6)\cos (\pi/6) = -2.3570

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = |\pi /6 - 0| = \pi/6 a = c = \pi/6 , b = \pi/2

iteración 2

a = \pi/6 , b = \pi/2 f(\pi/6) = -2.3570 f(\pi /2) = 2 c = \pi/2 - 2 \frac{\pi/6-\pi/2}{-1-2} = 1.0901 f(1.0901) = 3 \sin ^{3}(1.0901)-1 - 4 \sin (1.0901)\cos (1.0901) = -0.5486

el signo de f(c) es el mismo que f(a), se ajusta el lado izquierdo

tramo = |c-a| = | 1.0901-\pi/6| = 1.0722 a = c = 1.0901 , b = \pi/2

iteración 3

a = 1.0901 , b = \pi/2 f(1.0901) = -0.5486 f(\pi /2) = 2 c = \pi/2 - 2 \frac{-0.5486-\pi/2}{-0.5486-2} = 1.19358 f(1.19358) = 3 \sin ^{3}(1.19358)-1 - 4 \sin (1.19358)\cos (1.19358) = 0.0409

el signo de f(c) es el mismo que f(b), se ajusta el lado derecho

tramo = |b-c| = | \pi/2- 1.19358| = 0.3772 a = 1.0901 , b = 1.19358


Algoritmo con Python

Los parámetros aplicados en el algoritmo son los desarrollados en el planteamiento del problema e intervalo de búsqueda, con lo que se obtiene los siguientes resultados:

 raiz: 1.1864949811467547
error: 9.919955449055884e-05

las instrucciones en python son:

# 1Eva_IIT2009_T1 Movimiento de partícula en plano
import numpy as np
import matplotlib.pyplot as plt

#INGRESO
xt = lambda t: 3*(np.sin(t)**3)-1
yt = lambda t: 4*np.sin(t)*np.cos(t)
fx = lambda t: 3*(np.sin(t)**3)-1 - 4*np.sin(t)*np.cos(t) 

a = 0
b = np.pi/2
tolera = 0.001

# intervalo para gráfica
La = a
Lb = b
muestras = 21

# PROCEDIMIENTO
# Posicion Falsa
tramo = abs(b-a)
fa = fx(a)
fb = fx(b)
while not(tramo<=tolera):
    c = b - fb*(a-b)/(fa-fb)
    fc = fx(c)
    cambio = np.sign(fa)*np.sign(fc)
    if cambio>0:
        # actualiza en izquierda
        tramo = abs(c-a)
        a = c
        b = b
        fa = fc
    else:
        # actualiza en derecha
        tramo = abs(b-c)
        a = a
        b = c
        fb = fc

# para grafica
ti = np.linspace(La,Lb,muestras)
xi = xt(ti)
yi = yt(ti)
fi = fx(ti)

# SALIDA
print(' raiz:', c)
print('error:', tramo)
plt.plot(ti,xi, label='x(t)')
plt.plot(ti,yi, label='y(t)')
plt.plot(ti,fi, label='f(t)')
plt.plot(c,fx(c),'ro')
plt.axhline(0, color='green')
plt.axvline(c, color='magenta')
plt.legend()
plt.xlabel('t')
plt.title('Método de Falsa Posición')
plt.show()

s1Eva_IIT2010_T1 Aproximar con polinomio

Ejercicio: 1Eva_IIT2010_T1 Aproximar con polinomio

Desarrollo Analítico

Ejemplo para Lápiz  y Papel


Tarea 01 Semana 01 Fecha: año/mes/dia
Apellidos Nombres
Referencia: 1Eva_IIT2010_T1 Aproximar con polinomio

Opción 1

Para el ejemplo, supondremos que x0=0

El polinomio de Taylor requerido es de grado 2

P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 +

función f(x) y sus derivadas:

f(x) = e^x \cos (x) +1
f'(x) = e^x \cos (x) - e^x \sin(x) f'(x) = e^x (\cos (x) - \sin(x))
f''(x) = e^x( \cos (x) - \sin(x))+ + e^x (-\sin(x) - \cos(x)) f''(x) = -2 e^x \sin(x))

Punto x0 = 0 (ejemplo), dentro del intervalo.

Observación: escriba las expresiones, reemplazando los valores, asi si en la lección o examen no tuvo tiempo para usar la calculadora, se puede evaluar si realizaba las operaciones con el punto de referencia y expresiones correctas.

f(0) = e^0 \cos (0) +1 = 2 f'(0) = e^0(\cos (0) - \sin(0)) = 1 f''(0) = -2 e^0 \sin(0)) = 0

Sustitución en fórmula de polinomio:

p_2(x) = f(x_0) + \frac{f'(x_0)}{1!}(x-x_0) + \frac{f''(x_0)}{2!}(x-x_0)^2 P_{2}(x) = 2+\frac{1}{1} (x-0) + + \frac{0}{2}(x-0)^2 + P_{2}(x) = 2+ x

Tarea: realizar el ejercicio para x0 = π/2, verificando que pase por los puntos requeridos. Respuesta usando el algoritmo de Taylor:

p(x) = -4.81047738096535*x**2 + 10.3020830193353*x – 3.31309698247881


Opción 2

El polinomio requerido tiene la forma:

p(x) = a + bx + cx^2

por lo que conocido los pares ordenados por donde debe pasar se puede plantear las ecuaciones y encontrar a,b,c.

f(0) = e^0 \cos (0) +1 = 2 f(\pi/2) = e^{\pi/2} \cos (\pi /2) +1 = 1(0)+1 =1 f(\pi) = e^{\pi} \cos (\pi) +1 = e^{\pi} +1

se encuentra que a = 2 cuando x = 0 y que reemplazando los valores de x =π/2 y x=π se tiene:

2 + (π/2) b + (π/2)2 c = 1
2 +    π  b +   (π)2 c = eπ +1

que se convierte en:

(π/2) b + (π/2)2 c = -1
   π  b +   (π)2 c = -(eπ +1)

al multiplicar la primera ecuacion por 2 y restando de la segunda

- π2/2 c = eπ -1
c = (-2/π2)(eπ -1)

y sustituir c en la segunda ecuación:

π b + (π)2 (-2/π2)(eπ -1) = -(eπ +1)
π b = -(eπ +1) + 2(eπ -1) = -eπ -1 + 2eπ -2
b = (eπ -3)/π

El polinomio resultante es:

p(x) = 2 + \frac{e^{\pi}-3}{\pi}x + \frac{-1(e^{\pi}-1)}{\pi ^2}x^2

Probando respuesta con los valores en la función y polinomio usando python, se encuentra que el polinomio pasa por los puntos. Al observar la gráfica observa que se cumple lo requerido pero visualiza el error de aproximación usando el método de la opción 2.


Desarrollo con Python

Algoritmo desarrollado en clase, usado como taller, modificado para el problema planteado.

Observación: Se reordena el algoritmo para mantener ordenados y separados los bloques de ingreso, procedimiento y salida. Así los bloques pueden ser convertidos fácilmente a funciones algoritmicas def-return.

Observe que la variable n se interprete correctamente como «términos» o «grados» del polinomio de Taylor.

# Aproximación Polinomio de Taylor alrededor de x0
# f(x) en forma simbólica con sympy
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO --------------------
x = sym.Symbol('x')
fx = sym.exp(x)*sym.cos(x) + 1

x0 = 0
n  = 3 # grado de polinomio

# Intervalo para Gráfica
a = 0
b = np.pi
muestras = 21

# PROCEDIMIENTO  -------------
# construye polinomio Taylor
k = 0 # contador de términos
polinomio = 0
while (k <= n):
    derivada   = fx.diff(x,k)
    derivadax0 = derivada.subs(x,x0)
    divisor   = np.math.factorial(k)
    terminok  = (derivadax0/divisor)*(x-x0)**k
    polinomio = polinomio + terminok
    k = k + 1

# forma lambda para evaluación numérica
fxn = sym.lambdify(x,fx,'numpy')
pxn = sym.lambdify(x,polinomio,'numpy')

# evaluar en intervalo para gráfica
xi = np.linspace(a,b,muestras)
fxi = fxn(xi)
pxi = pxn(xi)

# SALIDA  --------------------
print('polinomio p(x)=')
print(polinomio)
print()
sym.pprint(polinomio)

# Gráfica
plt.plot(xi,fxi,label='f(x)')
plt.plot(xi,pxi,label='p(x)')
# franja de error
plt.fill_between(xi,pxi,fxi,color='yellow')
plt.xlabel('xi')
plt.axvline(x0,color='green', label='x0')
plt.axhline(0,color='grey')
plt.title('Polinomio Taylor: f(x) vs p(x)')
plt.legend()
plt.show()

Resultado del algoritmo

Revisar si el polinomio es concordante con lo realizado a lápiz y papel, de no ser así revisar el algoritmo o los pasos realizados en papel, deben ser iguales.
Comprobando que el algoritmo esté correcto y pueda ser usado en otros ejercicios.

 RESTART: D:\MATG1052Ejemplos\Taylot01Tarea01.py 
polinomio p(x)=
-x**3/3 + x + 2

   3        
  x         
- -- + x + 2
  3         

Resultados gráficos para x0=0

Continuar con el ejercicio con x0 = π y luego con el siguiente punto x0 = π/2.

Comparar resultados y presentar: Observaciones  y recomendaciones semejantes a las indicadas durante el desarrollo de la clase.

s3Eva_IIT2019_T4 completar polinomio de interpolación

Ejercicio: 3Eva_IIT2019_T4 completar polinomio de interpolación

Para visualizar la solución, se plantea graficar los polinomios que están completos S0(x) y S2(x).

S_0(x) = 1 + 1.1186x + 0.6938 x^3

 0.0 ≤ x ≤ 0.4

S_1(x) = 1.4918 + 1.4516(x-0.4) + c(x-0.4)^2 +d(x-0.4)^3

0.4 ≤ x ≤ 0.6

S_2(x) = 1.8221 + 1.8848(x-0.6) + +1.3336(x-0.6)^2 - 1.1113(x-0.6)^3

0.6 ≤ x ≤ 1.0

Como en trazadores cúbicos, lo que se usa es un polinomio por cada tramo muestreado para una curva contínua, etc. se tiene que los polinomios deben tener valores iguales en los puntos el eje x = 0.4 y 0.6

Por lo que se evalua con  los polinomios completos:

S_0(0.4) = 1 + 1.1186(0.4) + 0.6938 (0.4)^3 = 1.4918432 S_2(0.6) = 1.8221 + 1.8848(0.6-0.6) + +1.3336(0.6-0.6)^2-1.1113(0.6-0.6)^3=1.8221

Opción 1

Valores que se usan en los extremos del polinomio S1(x) para crear un sistema de dos ecuaciones y determinar los valores de c y d, completando el polinomio.

S_1(0.4) = 1.4918 + 1.4516(0.4-0.4) + c(0.4-0.4)^2 +d(0.4-0.4)^3

como los términos de c y d se hacen cero, hace falta una ecuación.

S_1(0.6) = 1.4918 + 1.4516(0.6-0.4) + + c(0.6-0.4)^2 +d(0.6-0.4)^3 = 1.8221 1.8221 = 1.4918 + 1.4516(0.2) + c(0.2)^2 +d(0.2)^3 0,03998 = 0.04c +0,008d

la otra ecuación se podría obtener usando la propiedad que las primeras derivadas de los polinomios deben ser iguales en los puntos x=0,4 y x= 0.6

S’0(0.2) =S’1(0.2)

Tarea: Desarrollar la siguiente ecuación y resolver


Opción 2

Si no recuerda la propiedad anterior, puede optar por usar otros conceptos para aproximar el resultado.

Si para el tramo en que se busca el polinomio se puede retroceder un tamaño de paso x = 0.2 y evualuar usando S0(0.2), se obtiene otropunto de referencia para crear un polinomio que pase por los mismos puntos.

S_0(0.2) = 1 + 1.1186*(0.2) + 0.6938 (0.2)^3

Se aplica lo mismo para un tamaño de paso más adelante de x = 0.6 es x = 0.8m se evalua S2(0.8) y se tienen suficientes puntos para usar cualquier método de interpolación y determinar el polinomio para el tramo faltante.

S_2(0.8) = 1.8221 + 1.8848(0.8-0.6) + +1.3336(0.8-0.6)^2 - 1.1113(0.8-0.6)^3
xi     = [0.        0.2        0.4      ]
S0(xi) = [1.        1.2292704  1.4918432]
xi     = [0.6       0.8        1.       ]
S2(xi) = [1.8221    2.2435136  2.7182728]

que permite hacer una tabla de puntos, y usando por ejemplo el método de interpolación de Lagrange  con x entre [0.2, 0.8] se obtiene otra forma del polinomio buscado:

p(x)=1.2292704 \frac{(x-0.4)(x-0.6)(x-0.8)}{(0.2-0.4)(0.2-0.6)(0.2-0.8)} + +1.4918432\frac{(x-0.2)(x-0.6)(x-0.8)}{0.4-0.2)(0.4-0.6)(0.4-0.8)}+ +1.8221\frac{(x-0.2)(x-0.4)(x-0.8)}{(0.6-0.2)(0.6-0.4)(0.6-0.8)} + + 2.2435136\frac{(x-0.2)(x-0.4)(x-0.6)}{(0.8-0.2)(0.8-0.4)(0.8-0.6)}

Tarea: continuar con el desarrollo


El literal b, requiere usar un metodo de búsqueda de raíces, para el cual se puede usar incluso bisección.

Tarea: continuar con el desarrollo

s3Eva_IIT2019_T3 Preparación de terreno en refineria

Ejercicio: 3Eva_IIT2019_T3 Preparación de terreno en refineria

Se requiere usar el nivel inicial en la matriz, para restar del nivel requerido que es constante 220

Nivel inicio (m) 0 50 100 150 200
0 241 239 238 236 234
25 241 239 237 235 233
50 241 239 236 234 231
75 242 239 236 232 229
100 243 239 235 231 227

lo que genera la matriz de diferencias. El valor es positivo indica remoción, el valor negativo indica por rellenar.

Diferencia (m) 0 50 100 150 200
0 21 19 18 16 14
25 21 19 17 15 13
50 21 19 16 14 11
75 22 19 16 12 9
100 23 19 15 11 7

El volumen se puede calcular por un método en cada fila, y luego los resultados por columnas por otro método o el mismo.
Por ejemplo Simpson de 1/3

I= \frac{hx}{3}(f(x_0) +4f(x_1)+f(x_2))

con lo que se obtiene:

I_{fila}(0) = \frac{50}{3}(21 +4(19)+18) +\frac{50}{3}(18 +4(16)+14) =24750 I_{fila}(25) = = \frac{50}{3}(21 +4(19)+17) + \frac{50}{3}(17 +4(15)+13) = 23383,33 I_{fila}(50) = \frac{50}{3}(21 +4(19)+16) + \frac{50}{3}(16 +4(14)+11) = 22000 I_{fila}(75) = \frac{50}{3}(22 +4(19)+16) + \frac{50}{3}(16 +4(12)+5) =21850 I_{fila}(100) = \frac{50}{3}(23 +4(19)+15) + \frac{50}{3}(15 +4(11)+7) = 20483,33

y usando el otro eje, se completa el volumen usando dos veces simpson:

Volumen = \frac{h_y}{3}(f(x_0) +4f(x_1)+f(x_2)) Remover = \frac{25}{3}(24750 +4(23383,33)+22000) + + \frac{25}{3}(22000 +4(21850)+20483,33)=2251388,89

El signo lo trae desde la diferencia, y muestra el sentido del desnivel.

Se adjunta la gráfica de superficie en azul como referencia del signo,  respecto al nivel requerido en color verde.

Error de truncamiento

la cota del error de truncamiento se estima como O(h5)

error_{trunca} = -\frac{h^5}{90} f^{(4)}(z)

para un valor de z entre [a,b]

para cuantificar el valor, se puede usar la diferencia finita Δ4f, pues con la derivada sería muy laborioso.

s3Eva_IIT2019_T2 Diferenciación, valor en frontera

Ejercicio: 3Eva_IIT2019_T2 Diferenciación, valor en frontera

La ecuación de problema de valor de frontera, con h = 1/4 = 0.25:

y'' = -(x+1)y' + 2y + (1-x^2) e^{-x}

0 ≤ x ≤ 1
y(0) = -1
y(1) = 0

Se interpreta en la tabla como los valores que faltan por encontrar:

i 0 1 2 3 4
xi 0 1/4 1/2 3/4 1
yi -1 0

Por ejemplo, se usan las derivadas en diferencias finitas divididas centradas para segunda derivada y hacia adelante para primera derivada.

Semejante al procedimiento aplicado para métodos con EDO.

f''(x_i) = \frac{f(x_{i+1})-2f(x_{i})+f(x_{i-1})}{h^2} + O(h^2) f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} + O(h)

Se formula entonces la ecuación en forma discreta, usando solo los índices para los puntos yi:

\frac{y[i+1]-2y[i]+y[i-1]}{h^2} = -(x_i +1) \frac{y[i+1]-y[i]}{h} + 2y[i] + (1-x_i ^2) e^{-x_i}

se multiplica todo por h2

y[i+1]-2y[i]+y[i-1] = -h(x_i +1)(y[i+1]-y[i]) + + 2 h^2 y[i] + h^2 (1-x_i ^2) e^{-x_i}
y[i+1]-2y[i]+y[i-1] = -h(x_i +1)y[i+1] + + h(x_i +1)y[i] + 2 h^2 y[i] + h^2 (1-x_i ^2) e^{-x_i}
y[i+1](1 +h(x_i +1)) + y[i](-2 - h(x_i +1)-2 h^2)+y[i-1] = = h^2 (1-x_i ^2) e^{-x_i}

Ecuación que se aplica en cada uno de los puntos desconocidos con i =1,2,3


i = 1

y[2](1 +h(x_1 +1)) + y[1](-2 - h(x_1 +1)-2 h^2)+y[0] = = h^2 (1-x_1 ^2) e^{-x_1} y[2]\Big(1 +\frac{1}{4}\Big(\frac{1}{4} +1\Big)\Big) + y[1]\Big(-2 - \frac{1}{4}\Big(\frac{1}{4} +1\Big)-2 \Big(\frac{1}{4}\Big)^2\Big) -1 = = \Big(\frac{1}{4}\Big)^2 \Big(1-\Big(\frac{1}{4}\Big) ^2\Big) e^{-1/4} y[2]\Big(1 +\frac{5}{16} \Big) + y[1]\Big(-2 - \frac{5}{16}- \frac{2}{16}\Big) = = 1+ \frac{1}{16} \Big(1-\frac{1}{16}\Big) e^{-1/4} \frac{21}{16} y[2]- \frac{39}{16}y[1] = 1+ \frac{15}{16^2} e^{-1/4} - \frac{39}{16}y[1] + \frac{21}{16} y[2] = 1+ \frac{15}{16^2} e^{-1/4}

multiplicando ambos lados de la ecuacion por 16 y reordenando

- 39 y[1] + 21 y[2] = 16+ \frac{15}{16} e^{-1/4}

i = 2

y[3](1 +h(x_2 +1)) + y[2](-2 - h(x_2 +1)-2 h^2)+y[1] = = h^2 (1-x_2 ^2) e^{-x_2} y[3]\Big(1 +\frac{1}{4}\Big(\frac{1}{2} +1\Big)\Big) + y[2]\Big(-2 - \frac{1}{4}\Big(\frac{1}{2} +1\Big)-2 \Big(\frac{1}{4}\Big)^2\Big)+y[1] = = \Big(\frac{1}{4}\Big)^2 \Big(1-\Big(\frac{1}{2}\Big) ^2\Big) e^{-\frac{1}{2}} y[3]\Big(1 +\frac{3}{8}\Big) + y[2]\Big(-2 - \frac{1}{4}\Big(\frac{3}{2}\Big)- \frac{1}{8}\Big)+y[1] = = \Big(\frac{1}{16}\Big) \Big(1-\frac{1}{4}\Big) e^{-\frac{1}{2}} \frac{11}{8} y[3] - \frac{21}{8} y[2]+y[1] = \frac{1}{16} \Big(\frac{3}{4}\Big) e^{-\frac{1}{2}}

multiplicando ambos lados por 8 y reordenando,

8y[1] - 21 y[2] + 11 y[3] = \frac{3}{8} e^{-\frac{1}{2}}

i = 3

y[4](1 +h(x_3 +1)) + y[3](-2 - h(x_3 +1)-2 h^2)+y[2] = = h^2 (1-x_3 ^2) e^{-x_3} (0) \Big(1 +h\Big(x_3 +1\Big)\Big) + y[3]\Big(-2 - \frac{1}{4}\Big(\frac{3}{4} +1\Big)-2 \frac{1}{16}\Big)+y[2] = = \frac{1}{16}\Big (1-\Big(\frac{3}{4}\Big) ^2 \Big) e^{-3/4} y[3]\Big(-2 - \frac{7}{16} - \frac{2}{16}\Big)+y[2] = \frac{1}{16}\Big (1-\frac{9}{16}\Big) e^{-3/4}

multiplicando todo por 16 y reordenando:

16 y[2] - 41 y[3] = \frac{7}{16} e^{-3/4}

Con lo que se puede crear la forma matricial del sistema de ecuaciones Ax=B

\begin{bmatrix} -39 && 21 && 0\\8 && 21 && 11 \\ 0 && 16 && 41 \end{bmatrix}\begin{bmatrix} y[1]\\ y[2] \\y[3] \end{bmatrix} =\begin{bmatrix} 16+ \frac{15}{16} e^{-1/4} \\ \frac{3}{8} e^{-\frac{1}{2}} \\ \frac{7}{16} e^{-3/4} \end{bmatrix}

con lo que resolviendo la matriz se obtienen los valroes para y[1], y[2], y[3]

solución de matriz: 
[-0.59029143 -0.29958287 -0.12195088]

con lo que se completan los puntos de la tabla,

Solución de ecuación
x[i] =  [ 0.   0.25        0.5         0.75        1.  ]
y[i] =  [-1.  -0.59029143 -0.29958287 -0.12195088  0.  ]

con la siguiente gráfica:


Algoritmo para solución de matriz con Python

tarea: modificar para cambiar el valor del tamaño de paso.

# Problema de Frontera
import numpy as np
import matplotlib.pyplot as plt

# INGRESO

h = 1/4
y0 = -1
y1 = 0

xi = np.arange(0,1+h,h)
n = len(xi)
yi = np.zeros(n,dtype=float)
yi[0] = y0
yi[n-1] = y1
    
A = np.array([[-39.0, 21,  0],
              [  8.0, -21, 11],
              [  0.0, 16, -41]])
B = np.array([16+(15/16)*np.exp(-1/4),
              (3/8)*np.exp(-1/2),
              (7/16)*np.exp(-3/4)])

# PROCEDIMIENTO
x = np.linalg.solve(A,B)

for i in range(1,n-1,1):
    yi[i] = x[i-1]

# SALIDA
print('solución de matriz: ')
print(x)
print('Solución de ecuación')
print('x[i] = ',xi)
print('y[i] = ',yi)

# Grafica
plt.plot(xi,yi,'ro')
plt.plot(xi,yi)
plt.show()

s3Eva_IIT2019_T1 Lanzamiento de Cohete

Ejercicio: 3Eva_IIT2019_T1 Lanzamiento de Cohete

A partir de la tabla del enunciado  se realiza la tabla de diferencias finitas.

i ti fi Δfi Δ2fi Δ3fi Δ4fi Δ5fi
1 0 0 32 -6 0 0 0
2 25 32 26 -6 0 0
3 50 58 20 -6 0
4 75 78 14 -6
5 100 92 8
6 125 100

Observando que a partir de la tercera diferencia finita  los valores son cero, por lo que se usa la fórmula general de diferencias finitas divididas hasta el polinomio de grado 2.

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

al sustituir los valores conocidos, se convierte en,

p_2 (t) =0 + \frac{32}{25} (t -0) + + \frac{-6}{2(25)^2} (t -0)(t - 25) =\frac{32}{25}t + \frac{-3}{(25)^2} (t^2 - 25t) =\frac{32}{25}t + \frac{-3}{(25)^2} t^2 - \frac{-3}{(25)^2}25t =\frac{7}{5}t - \frac{3}{625} t^2 y(t) =p_2 (t) =1.4 t - 0.0048 t^2

Con lo que se puede obtener la velocidad:

y'(t) = 1.4 - 0.0096 t

y luego la aceleración:

y''(t) = - 0.0096

Si el error es el próximo término del polinomio Δ3fi  entonces se estima en cero.

Tarea:  Evaluar la velocidad y aceleración para cada punto de la tabla

La gráfica del polinomio encontrado es:

Algoritmo en Python

El algoritmo realizado en Python entrega los siguientes resultados:

[[  i,  ti,  fi, df1, df2, df3, df4, df5,  df6]]
[[  1.   0.   0.  32.  -6.   0.   0.   0.   0.]
 [  2.  25.  32.  26.  -6.   0.   0.   0.   0.]
 [  3.  50.  58.  20.  -6.   0.   0.   0.   0.]
 [  4.  75.  78.  14.  -6.   0.   0.   0.   0.]
 [  5. 100.  92.   8.   0.   0.   0.   0.   0.]
 [  6. 125. 100.   0.   0.   0.   0.   0.   0.]]
polinomio:
-0.0048*t**2 + 1.4*t

las instrucciones en Python son:

# 3Eva_IIT2019_T1 Lanzamiento de Cohete
# Tarea: Verificar tamaño de vectores
#        considerar puntos no equidistantes en eje t
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO , Datos de prueba
ti = np.array([0.0, 25, 50, 75, 100, 125])
fi = np.array([0.0, 32, 58, 78, 92, 100])

# PROCEDIMIENTO
# Tabla de diferencias finitas
titulo = ['i','ti','fi']
n = len(ti)

# cambia a forma de columnas
i = np.arange(1,n+1,1)
i = np.transpose([i])
ti = np.transpose([ti])
fi = np.transpose([fi])

# Añade matriz de diferencias
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((i,ti,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 t equidistantes
dfinita = tabla[:,3:]
n = len(dfinita)
t = sym.Symbol('t')
h = ti[1,0]-ti[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*(t-ti[f])
    polinomio = polinomio + termino*factor

# simplifica polinomio, multiplica los (t-ti)
polinomio = polinomio.expand()

# para evaluacion numérica
pt = sym.lambdify(t,polinomio)

# Puntos para la gráfica
a = np.min(ti)
b = np.max(ti)
muestras = 101
ti_p = np.linspace(a,b,muestras)
fi_p = pt(ti_p)

# SALIDA
print([titulo])
print(tabla)
print('polinomio:')
print(polinomio)

# Gráfica
plt.title('Interpolación polinómica')
plt.plot(ti,fi,'o', label = 'Puntos')
plt.plot(ti_p,fi_p, label = 'Polinomio')
plt.legend()
plt.show()