Publicado en

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
Publicado en

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 3 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 diferencias 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 construir 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()
Publicado en

s3Eva_2021PAOI_T3 Respuesta a entrada cero en un sistema LTIC

Ejercicio: 3Eva_2021PAOI_T3 Respuesta a entrada cero en un sistema LTIC

la ecuación a resolver es:

\frac{\delta^2 y(t)}{\delta t^2}+3 \frac{\delta y(t)}{ \delta t}+2 y(t) =0

con valores iniciales: y0(t)=0 , y’0(t) =-5

se puede escribir como:

y"+3 y'+2y = 0 y" = -3y'-2y

sutituyendo las expresiones de las derivadas como las funciones f(x) por z y g(x) por z’:

y’ = z = f(x)

y» = z’= -3z-2y = g(x)

Los valores iniciales de t0=0, y0=0, z0=-5 se usan en el algoritmo.

En este caso también se requiere conocer un intervalo de tiempo a observar [0,tn=6] y definir el tamaño de paso o resolución del resultado

h = \delta t = \frac{b-a}{n} = \frac{6-0}{60} = 0.1

t0 = 0, y0 = 0,  y’0 = z0 = -5

iteración 1

K1y = h * f(ti,yi,zi) = 0.1 (-5) = -0.5

K1z = h * g(ti,yi,zi) ) = 0.1*(-3(-5)-2(0)) = 1.5

K2y = h * f(ti+h, yi + K1y, zi + K1z) = 0.1(-5+1.5) = -0.35

K2z = h * g(ti+h, yi + K1y, zi + K1z) = 0.1 ( -3(-5+1.5) – 2(0-0.5)) = 1.15

yi = yi + (K1y+K2y)/2 =0+(-0.5-0.35)/2=-0.425

zi = zi + (K1z+K2z)/2 = -5+(1.5+1.15)/2 = -3.675

ti = ti + h = 0+0.1 = 0.1

iteración 2

K1y = 0.1 (-3.675) = -0.3675

K1z = 0.1*(-3(-3.675)-2(-0.425)) = 1.1875

K2y = 0.1(-3.675+ 1.1875) = -0.24875

K2z = 0.1 ( -3(-3.675+ 1.1875) – 2(-0.425-0.3675)) = 0.90475

yi = -0.425+(-0.3675-0.24875)/2=-0.7331

zi = -3.675+( 1.1875+0.90475)/2 = -2.6288

ti =0.1+0.1 = 0.2

iteración 3

continuar como ejercicio

El algoritmo permite obtener la gráfica y la tabla de datos.

los valores de las iteraciones son:

t, y, z
[[ 0.        0.       -5.      ]
 [ 0.1      -0.425    -3.675   ]
 [ 0.2      -0.733125 -2.628875]
 [ 0.3      -0.949248 -1.807592]
 [ 0.4      -1.093401 -1.167208]
 [ 0.5      -1.18168  -0.67202 ]
 [ 0.6      -1.226984 -0.293049]
 [ 0.7      -1.239624 -0.006804]
 [ 0.8      -1.227806  0.205735]
 [ 0.9      -1.19804   0.359943]
 [ 1.       -1.155465  0.468225]
 [ 1.1      -1.104111  0.540574]
 [ 1.2      -1.047121  0.585021]
 [ 1.3      -0.986923  0.608001]
 [ 1.4      -0.925374  0.614658]
 [ 1.5      -0.863874  0.609087]
 [ 1.6      -0.803463  0.594537]
 [ 1.7      -0.744893  0.573574]
 [ 1.8      -0.68869   0.548208]
 [ 1.9      -0.635205  0.520011]
 [ 2.       -0.584652  0.490193]
 [ 2.1      -0.53714   0.459683]
 [ 2.2      -0.492695  0.42918 ]
 [ 2.3      -0.451288  0.399206]
 [ 2.4      -0.412843  0.370135]
 [ 2.5      -0.377253  0.342233]
 [ 2.6      -0.34439   0.315674]
 [ 2.7      -0.314114  0.290567]
 [ 2.8      -0.286275  0.266966]
 [ 2.9      -0.26072   0.244887]
 [ 3.       -0.237297  0.224314]
 [ 3.1      -0.215858  0.205211]
 [ 3.2      -0.196256  0.187526]
 [ 3.3      -0.178354  0.171195]
 [ 3.4      -0.162019  0.156149]
 [ 3.5      -0.147126  0.142312]
 [ 3.6      -0.133558  0.129611]
 [ 3.7      -0.121206  0.117969]
 [ 3.8      -0.109966  0.107312]
 [ 3.9      -0.099745  0.097569]
 [ 4.       -0.090454  0.08867 ]
 [ 4.1      -0.082013  0.080549]
 [ 4.2      -0.074346  0.073146]
 [ 4.3      -0.067385  0.066401]
 [ 4.4      -0.061067  0.06026 ]
 [ 4.5      -0.055334  0.054673]
 [ 4.6      -0.050134  0.049591]
 [ 4.7      -0.045417  0.044972]
 [ 4.8      -0.04114   0.040776]
 [ 4.9      -0.037263  0.036964]
 [ 5.       -0.033748  0.033503]
 [ 5.1      -0.030563  0.030362]
 [ 5.2      -0.027677  0.027512]
 [ 5.3      -0.025062  0.024926]
 [ 5.4      -0.022692  0.022581]
 [ 5.5      -0.020546  0.020455]
 [ 5.6      -0.018602  0.018527]
 [ 5.7      -0.016841  0.01678 ]
 [ 5.8      -0.015246  0.015196]
 [ 5.9      -0.013802  0.013761]
 [ 6.       -0.012494  0.012461]]

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y

t0 = 0
y0 = 0
z0 = -5

h = 0.1
tn = 6
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print('t, y, z')
print(tabla)

# GRAFICA
plt.plot(ti,yi, label='y(t)')

plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()
Publicado en

s3Eva_2021PAOI_T2 Tensiones mínimas en cables por carga variable

Ejercicio: 3Eva_2021PAOI_T2 Tensiones mínimas en cables por carga variable

El ejercicio usa el resultado del tema anterior, planteando una función de Python como la solución para valores dados. Se requiere una función, para disponer de los valores solución en cada llamada para el intervalo de análisis.

Por lo que básicamente lo que se pide es usar algún algoritmo de búsqueda de raíces. Para simplicidad en la explicación se toma el algoritmo de la bisección.

Los resultados se grafican como theta vs Tensiones, y el punto a buscar es cuando las tensiones en los cables son de igual magnitud, es decir:

TCA = TCB

Resultando en :

Resultado: [TCA, TCB], diferencia
[array([3.46965006e-14, 4.00000000e+02]), -399.99999999999994]
tetha, TCA, TCB
[[-2.61799388e-01  3.46965006e-14  4.00000000e+02]
 [-1.74532925e-01  3.70996817e+01  3.85789041e+02]
 [-8.72664626e-02  7.39170124e+01  3.68641994e+02]
 [ 8.32667268e-17  1.10171790e+02  3.48689359e+02]
 [ 8.72664626e-02  1.45588094e+02  3.26082988e+02]
 [ 1.74532925e-01  1.79896384e+02  3.00994928e+02]
 [ 2.61799388e-01  2.12835554e+02  2.73616115e+02]
 [ 3.49065850e-01  2.44154918e+02  2.44154918e+02]
 [ 4.36332313e-01  2.73616115e+02  2.12835554e+02]
 [ 5.23598776e-01  3.00994928e+02  1.79896384e+02]
 [ 6.10865238e-01  3.26082988e+02  1.45588094e+02]
 [ 6.98131701e-01  3.48689359e+02  1.10171790e+02]
 [ 7.85398163e-01  3.68641994e+02  7.39170124e+01]
 [ 8.72664626e-01  3.85789041e+02  3.70996817e+01]]
       raiz en:  0.34898062924398343 radianes
       raiz en:  19.995117187500004 grados
error en tramo:  8.522115488257542e-05
>>>

Instrucciones en Python

se añaden las instrucciones de la bisección al algoritmo del tema anterior, para encontrar el punto de intersección,

import numpy as np
import matplotlib.pyplot as plt

# Tema 1
def funcion(P,theta,alfa,beta):
    # ecuaciones
    A = np.array([[np.cos(alfa), -np.cos(beta)],
                  [np.sin(alfa),  np.sin(beta)]])
    B = np.array([P*np.sin(theta), P*np.cos(theta)])

    # usar un algoritmo directo
    X = np.linalg.solve(A,B)
    
    diferencia = X[0]-X[1]
    return([X,diferencia])    

# INGRESO
alfa = np.radians(35)
beta = np.radians(75)
P = 400

# PROCEDIMIENTO
theta = beta-np.radians(90)
resultado = funcion(P,theta,alfa, beta)

# SALIDA
print("Resultado: [TCA, TCB], diferencia")
print(resultado)

# Tema 1b --------------
# PROCEDIMIENTO
dtheta = np.radians(5)
theta1 = beta-np.radians(90)
theta2 = np.radians(90)-alfa

tabla = []
theta = theta1
while not(theta>=theta2):
    X = funcion(P,theta,alfa,beta)[0] # usa vector X
    tabla.append([theta,X[0],X[1]])
    theta = theta + dtheta
    
tabla = np.array(tabla)
thetai = np.degrees(tabla[:,0])
Tca = tabla[:,1]
Tcb = tabla[:,2]

# SALIDA
print('tetha, TCA, TCB')
print(tabla)

# Grafica
plt.plot(thetai,Tca, label='Tca')
plt.plot(thetai,Tcb, label='Tcb')
# plt.axvline(np.degrees(c))
plt.legend()
plt.xlabel('theta')
plt.ylabel('Tensión')
plt.show()


# Tema 2 -------------------------
# busca intersección con Bisección
diferencia = Tca-Tcb
donde_min  = np.argmin(np.abs(diferencia))
a = tabla[donde_min-1,0]
b = tabla[donde_min+1,0]
tolera = 0.0001

tramo = b-a
while not(tramo<tolera):
    c = (a+b)/2
    fa = funcion(P,a,alfa,beta)[1] # usa delta
    fb = funcion(P,b,alfa,beta)[1]
    fc = funcion(P,c,alfa,beta)[1]
    cambia = np.sign(fa)*np.sign(fc)
    if cambia < 0: 
        a = a
        b = c
    if cambia > 0:
        a = c
        b = b
    tramo = b-a

# SALIDA
print('       raiz en: ', c,"radianes")
print('       raiz en: ', np.degrees(c),"grados")
print('error en tramo: ', tramo)

# Grafica
plt.plot(thetai,Tca, label='Tca')
plt.plot(thetai,Tcb, label='Tcb')
plt.axvline(np.degrees(c))
plt.legend()
plt.xlabel('theta')
plt.ylabel('Tensión')
plt.show()
Publicado en

s3Eva_2021PAOI_T1 Tensiones en cables por carga variable

Ejercicio: 3Eva_2021PAOI_T1 Tensiones en cables por carga variable

Planteamiento del problema

Las ecuaciones de equilibrio del sistema corresponden a:

-T_{CA} \cos (\alpha) + T_{CB} \cos (\beta) + P \sin (\theta) = 0 T_{CA} \sin (\alpha) + T_{CB} \sin (\beta) - P \cos (\theta) = 0

se reordenan considerando que P y θ son valores constantes para cualquier caso

T_{CA} \cos (\alpha) - T_{CB} \cos (\beta) = P \sin (\theta) T_{CA} \sin (\alpha) + T_{CB} \sin (\beta) = P \cos (\theta)

se convierte a la forma matricial

\begin{bmatrix} \cos (\alpha) & -\cos (\beta) \\ \sin (\alpha) & \sin (\beta) \end{bmatrix} \begin{bmatrix} T_{CA} \\ T_{CB} \end{bmatrix} = \begin{bmatrix} P \sin (\theta) \\ P \cos (\theta) \end{bmatrix}

tomando valores por ejemplo:

α = 35°, β = 75°, P = 400 lb, Δθ = 5°

θ = 75-90 = -15

\begin{bmatrix} \cos (35) & -\cos (75) \\ \sin (35) & \sin (75) \end{bmatrix} \begin{bmatrix}T_{CA} \\ T_{CB} \end{bmatrix} = \begin{bmatrix} 400 \sin (-15) \\ 400 \cos (15) \end{bmatrix}

Desarrollo analítico

matriz aumentada

\begin{bmatrix} \cos (35) & - \cos (75) & 400 \sin (-15) \\ \sin (35) & \sin (75 ) & 400 \cos (15 ) \end{bmatrix} 
A = 
[[ 0.81915204 -0.25881905]
 [ 0.57357644  0.96592583]]
B = 
[-103.52761804  386.37033052]

AB = 
[[ 0.81915204 -0.25881905 -103.52761804]
 [ 0.57357644  0.96592583 386.37033052]]

Pivoteo parcial por filas

cos(-15°) tendría mayor magnitud que sin(-15°), por lo que la matriz A se encuentra pivoteada

Eliminación hacia adelante

pivote =  0.81915204/0.57357644
[[ 0.81915204 -0.25881905 -103.52761804]
 [ 0.0         1.63830408  655.32162903]]

Sustitución hacia atras

usando la última fila:

1.63830408 TCB = 655.32162903
TCB = 400

luego la primera fila:

0.81915204TCA -0.25881905TCB = -103.52761804

0.81915204TCA = 0.25881905TCB  -103.52761804

TCA = 2,392 x10-6 ≈ 0

Se deja como tarea realizar el cálculo para:  θ+Δθ

Instrucciones en Python

Resultado:

Resultado: [TCA, TCB], diferencia 
[array([3.46965006e-14, 4.00000000e+02]), -399.99999999999994]

usando el intervalo para θ1 y θ2:

con las instrucciones:

import numpy as np
import matplotlib.pyplot as plt

# Tema 1
def funcion(P,theta,alfa,beta):
    # ecuaciones
    A = np.array([[np.cos(alfa), -np.cos(beta)],
                  [np.sin(alfa),  np.sin(beta)]])
    B = np.array([P*np.sin(theta), P*np.cos(theta)])

    # usar un algoritmo directo
    X = np.linalg.solve(A,B)
    
    diferencia = X[0]-X[1]
    return([X,diferencia])    

# INGRESO
alfa = np.radians(35)
beta = np.radians(75)
P = 400

# PROCEDIMIENTO
theta = beta-np.radians(90)
resultado = funcion(P,theta,alfa, beta)

# SALIDA
print("Resultado: [TCA, TCB], diferencia")
print(resultado)

# Tema 1b --------------
# PROCEDIMIENTO
dtheta = np.radians(5)
theta1 = beta-np.radians(90)
theta2 = np.radians(90)-alfa

tabla = []
theta = theta1
while not(theta>=theta2):
    X = funcion(P,theta,alfa,beta)[0] # usa vector X
    tabla.append([theta,X[0],X[1]])
    theta = theta + dtheta
    
tabla = np.array(tabla)
thetai = np.degrees(tabla[:,0])
Tca = tabla[:,1]
Tcb = tabla[:,2]

# SALIDA
print('tetha, TCA, TCB')
print(tabla)

# Grafica
plt.plot(thetai,Tca, label='Tca')
plt.plot(thetai,Tcb, label='Tcb')
# plt.axvline(np.degrees(c))
plt.legend()
plt.xlabel('theta')
plt.ylabel('Tensión')
plt.show()
Publicado en

s2Eva_2021PAOI_T3 EDP Elíptica con valores en la frontera f(x) g(y)

Ejercicio: 2Eva_2021PAOI_T3 EDP Elíptica con valores en la frontera f(x) g(y)

Dada la EDP elíptica

\frac{\partial ^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2} = 0 0 \lt x \lt \frac{1}{2}, 0 \lt y\lt \frac{1}{2}

Se convierte a la versión discreta usando diferencias divididas centradas:

 

\frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} = 0

Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy2

\frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) + + \frac{\Delta y^2}{\Delta y^2}\Big(u[i,j-1]-2u[i,j]+u[i,j+1]\Big) = 0

los tamaños de paso en ambos ejes son de igual valor, se simplifica la ecuación

\lambda= \frac{\Delta y^2}{\Delta x^2} = 1
u[i-1,j]-2u[i,j]+u[i+1,j] + + u[i,j-1]-2u[i,j]+u[i,j+1] = 0
u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] = 0

que permite plantear las ecuaciones para cada punto en posición [i,j]

En cada iteración se requiere el uso de los valores en la frontera

u(x,0)=0, 0 \leq x \leq \frac{1}{2} u(0,y)=0 , 0\leq y \leq \frac{1}{2} u\Big(x,\frac{1}{2} \Big) = 200 x , 0 \leq x \leq \frac{1}{2} u\Big(\frac{1}{2} ,y \Big) = 200 y , 0 \leq y \leq \frac{1}{2}

Iteraciones

i=1, j=1

u[0,1]-4u[1,1]+u[2,1] + + u[1,0]+u[1,2] = 0

usando los valores en la frontera,

0-4u[1,1]+u[2,1] + 0+u[1,2] = 0 -4u[1,1]+u[2,1] + u[1,2] = 0

i=2, j=1

u[1,1]-4u[2,1]+u[3,1] + + u[2,0]+u[2,2] = 0

usando los valores en la frontera,

u[1,1]-4u[2,1]+200(1/6) + 0+u[2,2] = 0 u[1,1]-4u[2,1] + u[2,2] = -200(1/6)

i=1, j=2

u[0,2]-4u[1,2]+u[2,2] + + u[1,1]+u[1,3] = 0 0 - 4u[1,2]+u[2,2] + u[1,1]+200\frac{1}{6} = 0 - 4u[1,2] + u[2,2]+u[1,1] = -200\frac{1}{6}

i=2, j=2

u[1,2]-4u[2,2]+u[3,2] + + u[2,1]+u[2,3] = 0 u[1,2]-4u[2,2]+200\frac{2}{6} + u[2,1]+200\frac{2}{6} = 0 u[1,2]-4u[2,2]+ u[2,1] = -(2)200\frac{2}{6}

Sistema de ecuaciones a resolver:

\begin{bmatrix} -4 & 1 &1&0\\1 &-4&0&1\\1&0&-4&1 \\0&1&1&-4\end{bmatrix} \begin{bmatrix} u[1,1]\\u[2,1] \\u[1,2]\\u[2,2] \end{bmatrix} = \begin{bmatrix} 0\\-200(1/6)\\-200(1/6)\\-200(4/6) \end{bmatrix}

Resolviendo el sistema se tiene:

[11.11111111 22.22222222 22.22222222 44.44444444]

Instrucciones Python

import numpy as np

A = np.array([[-4, 1, 1, 0],
              [ 1,-4, 0, 1],
              [ 1, 0,-4, 1],
              [ 0, 1, 1,-4]])

B = np.array([0,-200*(1/6),-200*(1/6),-200*(4/6)])

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

print(x)
Publicado en

s2Eva_2021PAOI_T1 Masa transportada por tubo

Ejercicio: 2Eva_2021PAOI_T1 Masa transportada por tubo

Las expresiones siguientes se usan dentro de la expresión del integral

Q(t)=9+4 \cos ^2 (0.4t) c(t)=5e^{-0.5t}+2 e^{-0.15 t} M = \int_{t_1}^{t_2} Q(t)c(t) dt

literal a

Usando los valores dados para el intervalo [2,8] con 6 tramos h = (8-2)/6 =1

Se usa los valores de cada ti en Se puede obtener una tabla de valores muestreados para integrar f(t) = Q(i)c(t)

[ti,	 Qi,	 Ci,	 fi]
[ 2.     10.9416  3.321  36.3374]
[ 3.      9.5252  2.3909 22.7739]
[ 4.      9.0034  1.7743 15.9747]
[ 5.      9.6927  1.3552 13.1352]
[ 6.     11.175   1.0621 11.8687]
[ 7.     12.5511  0.8509 10.6793]
[ 8.     12.9864  0.694   9.0121]

Para el integral se usan los valores por cada dos tramos

I\cong \frac{1}{3}[36.3374+4(22.7739) + 15.9747] + \frac{1}{3}[15.9747+4(13.1352) + 11.8687] + \frac{1}{3}[11.8687+4(10.6793) + 9.0121] I = 95.7965

literal b

L acota de error de truncamiento por cada fórmula usada, se estima como O(h5),

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

para un valor de z entre [a,b]

por lo que al usar 3 veces la formula de Simpson se podría estimar en:

error_{trunca} = 3(1^5/90) = 0.033

literal c

El resultado se puede mejorar de dos formas:

1. Dado que el número de tramos es múltiplo de 3, se puede cambiar la fórmula a Simpon de 3/8, que tendría una cota de error menor

2. Aumentar el número de tramos disminuyendo el valor de h para que el error disminuya. Por ejemplo si se reduce a 0.5, el error disminuye en el orden de 0.55

Podría recomendar la segunda opión, pues a pesar que se aumenta la cota de error por cada vez que se usa la fórmula, el error de cada una disminuye en ordenes de magnitud 0,03125

La gráfica del ejercicio es:

Instrucciones en Python

import numpy as np
import matplotlib.pyplot as plt

# INGRESO
Q = lambda t: 9+4*(np.cos(0.4*t)**2)
C = lambda t: 5*np.exp(-0.5*t)+2*np.exp(-0.15*t)

t1 = 2
t2 = 8
n  = 6

# PROCEDIMIENTO
muestras = n+1 
dt = (t2-t1)/n
ti = np.arange(t1,t2+dt,dt)
Qi = Q(ti)
Ci = C(ti)
fi = Qi*Ci

# integración con Simpson 1/3
h= dt
I13 = 0
for i in range(0,6,2):
    S13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
    I13 = I13 + S13

# SALIDA
np.set_printoptions(precision=4)
print("[ti,\t Qi,\t Ci,\t fi]")
for i in range(0,muestras,1):
    print(np.array([ti[i],Qi[i],Ci[i],fi[i]]))
print("Integral S13: ",I13)
# grafica
plt.plot(ti,Qi, label = "Q(t)")
plt.plot(ti,Ci, label = "c(t)")
plt.plot(ti,fi, label = "f(t)")
plt.plot(ti,Qi,'.b')
plt.plot(ti,Ci,'.r')
plt.plot(ti,fi,'.g')
plt.xlabel("t")
plt.ylabel("f(t)=Q(t)*c(t)")
plt.show()
Publicado en

s2Eva_2021PAOI_T2 EDO para cultivo de peces

Ejercicio: 2Eva_2021PAOI_T2 EDO para cultivo de peces

Siendo la captura una constante mas una función periódica,

h(t) = a + b \sin (2 \pi t)

La ecuación EDO del ejercicio, junto a las constantes a=0.9 y b=0.75, r=1

\frac{\delta y(t)}{\delta t} = r y(t)-h(t)

se convierte en:

\frac{\delta y(t)}{\delta t} = (1) y(t)- \Big( 0.9 + .75 \sin (2 \pi t)\Big) \frac{\delta y(t)}{\delta t} = y(t)- 0.9 - .75 \sin (2 \pi t)

Considerando que la población inicial de peces es 1 o 100%, y(0)=1

literal a

h=1/12
tamano = muestras + 1
estimado = np.zeros(shape=(tamano,2),dtype=float)
estimado[0] = [0,1]
ti = 0
yi = 1
for i in range(1,tamano,1):
    K1 = 1/12 * d1y(ti,yi)
    K2 = 1/12 * d1y(ti+1/24, yi + K1/2)
    K3 = 1/12 * d1y(ti+1/24, yi + K2/2)
    K4 = 1/12 * d1y(ti+1/12, yi + K3)

    yi = yi + (1/6)*(K1+2*K2+2*K3 +K4)
    ti = ti + 1/12
        
    estimado[i] = [ti,yi]

literal b

iteración i=0

t(0) = 0

y(0) = 1

K1 = \frac{1}{12} \Big(1- 0.9 - .75 \sin (2 \pi 0)\Big) = 0,008333 K2 = \frac{1}{12} \Big(1- 0.9 - .75 \sin \Big(2 \pi (0+\frac{1}{12})\Big)\Big) = -0.02222 y(1) = 0 + \frac{0.008333+(-0.02222)}{2} = 0.9930 t(1) = 0 + \frac{1}{12} = \frac{1}{12}

iteración i=1

t(1) = \frac{1}{12}

y(1) = 0.9930

K1 = \frac{1}{12} \Big(0.9930 - 0.9 - .75 \sin \Big( 2 \pi\frac{1}{12}\Big)\Big) = -0.02349 K2 = \frac{1}{12} \Big(0.9930 - 0.9 - .75 \sin \Big(2 \pi (\frac{1}{12}+\frac{1}{12})\Big)\Big) = -0.04832 y(1) = 0.9930 + \frac{-0.02349+(-0.04832)}{2} = 0.9571 t(1) = \frac{1}{12} + \frac{1}{12} = \frac{2}{12}

iteración i=2

t(2) = \frac{2}{12}

y(1) = 0.9571

K1 = \frac{1}{12} \Big(0.9571 - 0.9 - .75 \sin \Big( 2 \pi\frac{2}{12}\Big)\Big) = -0.04936 K2 = \frac{1}{12} \Big(0.9571 - 0.9 - .75 \sin \Big(2 \pi (\frac{2}{12}+\frac{1}{12})\Big)\Big) = -0.06185 y(1) = 0.9571 + \frac{-0.04936+(-0.06185)}{2} = 0.9015 t(3) = \frac{2}{12} + \frac{1}{12} = \frac{3}{12}

literal c

Resultado del algoritmo, muestra que la estragegia de cosecha, en el tiempo no es sostenible, dado que la población de peces en el tiempo decrece.

estimado[xi,yi,K1,K2]
[[ 0.0000e+00  1.0000e+00  8.3333e-03 -2.2222e-02]
 [ 8.3333e-02  9.9306e-01 -2.3495e-02 -4.8330e-02]
 [ 1.6667e-01  9.5714e-01 -4.9365e-02 -6.1852e-02]
 [ 2.5000e-01  9.0153e-01 -6.2372e-02 -5.9196e-02]
 [ 3.3333e-01  8.4075e-01 -5.9064e-02 -4.1109e-02]
 [ 4.1667e-01  7.9066e-01 -4.0361e-02 -1.2475e-02]
 [ 5.0000e-01  7.6425e-01 -1.1313e-02  1.8994e-02]
 [ 5.8333e-01  7.6809e-01  2.0257e-02  4.4822e-02]
 [ 6.6667e-01  8.0063e-01  4.5845e-02  5.8039e-02]
 [ 7.5000e-01  8.5257e-01  5.8547e-02  5.5053e-02]
 [ 8.3333e-01  9.0937e-01  5.4907e-02  3.6606e-02]
 [ 9.1667e-01  9.5513e-01  3.5844e-02  7.5807e-03]
 [ 1.0000e+00  9.7684e-01  6.4031e-03 -2.4313e-02]
 [ 1.0833e+00  9.6788e-01 -2.5593e-02 -5.0602e-02]
 [ 1.1667e+00  9.2978e-01 -5.1645e-02 -6.4322e-02]
 [ 1.2500e+00  8.7180e-01 -6.4850e-02 -6.1881e-02]
 [ 1.3333e+00  8.0844e-01 -6.1757e-02 -4.4027e-02]
 [ 1.4167e+00  7.5554e-01 -4.3288e-02 -1.5645e-02]
 [ 1.5000e+00  7.2608e-01 -1.4494e-02  1.5549e-02]
 [ 1.5833e+00  7.2661e-01  1.6800e-02  4.1077e-02]
 [ 1.6667e+00  7.5554e-01  4.2089e-02  5.3969e-02]
 [ 1.7500e+00  8.0357e-01  5.4464e-02  5.0630e-02]
 [ 1.8333e+00  8.5612e-01  5.0470e-02  3.1799e-02]
 [ 1.9167e+00  8.9725e-01  3.1021e-02  2.3563e-03]
 [ 2.0000e+00  9.1394e-01  1.1619e-03 -2.9991e-02]
 [ 2.0833e+00  8.9953e-01 -3.1289e-02 -5.6773e-02]
 [ 2.1667e+00  8.5550e-01 -5.7835e-02 -7.1028e-02]
 [ 2.2500e+00  7.9107e-01 -7.1578e-02 -6.9169e-02]
 [ 2.3333e+00  7.2069e-01 -6.9069e-02 -5.1948e-02]
 [ 2.4167e+00  6.6018e-01 -5.1235e-02 -2.4254e-02]
 [ 2.5000e+00  6.2244e-01 -2.3130e-02  6.1924e-03]
 [ 2.5833e+00  6.1397e-01  7.4142e-03  3.0909e-02]
 [ 2.6667e+00  6.3313e-01  3.1888e-02  4.2918e-02]
 [ 2.7500e+00  6.7053e-01  4.3378e-02  3.8619e-02]
 [ 2.8333e+00  7.1153e-01  3.8421e-02  1.8746e-02]
 [ 2.9167e+00  7.4012e-01  1.7926e-02 -1.1830e-02]
 [ 3.0000e+00  7.4317e-01  0.0000e+00  0.0000e+00]]

Instrucciones en Python

# EDO. Método de RungeKutta 2do Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano   = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (K1+K2)/2
        xi = xi + h
        estimado[i-1,2:]=[K1,K2]
        estimado[i] = [xi,yi,0,0]
    return(estimado)

# PROGRAMA PRUEBA
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1

# INGRESO
# d1y = y' = f, d2y = y'' = f'
a =0.9; b=0.75; r=1
d1y = lambda t,y: r*y-(a+b*np.sin(2*np.pi*t))
x0 = 0
y0 = 1
h  = 1/12
muestras = 12*3

# PROCEDIMIENTO
puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
xi = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# SALIDA
np.set_printoptions(precision=4)
print('estimado[xi,yi,K1,K2]')
print(puntosRK2)


# Gráfica
import matplotlib.pyplot as plt


plt.plot(xi[0],yiRK2[0],
         'o',color='r', label ='[x0,y0]')
plt.plot(xi[1:],yiRK2[1:],
         'o',color='m',
         label ='y Runge-Kutta 2 Orden')

plt.title('EDO: Solución con Runge-Kutta 2do Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()
Publicado en

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)
Publicado en

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