Autor: Edison Del Rosario

Ejercicio: 2Eva_2021PAOII_T2 EDO – Embudos cónicos para llenar botellas

literal a

La expresión dada en el enunciado para EDO, se reordena para definir la funcion a usar con Runge-Kutta:

\frac{\delta y(t)}{\delta t} + \frac{d^2}{4}\sqrt{2 g \text{ }y(t)}\Bigg[\frac{tan \theta}{y(t)} \Bigg]^2 = 0 \frac{\delta y(t)}{\delta t} = - \frac{d^2}{4}\sqrt{2 g \text{ }y(t)}\Bigg[\frac{tan \theta}{y(t)} \Bigg]^2

siendo h = 0.5,  con y(0) = 0.15 m y d= 0.01 m ajustando las unidades de medida.

\frac{\delta y(t)}{\delta t} = - \frac{0.01^2}{4}\sqrt{2 (9.8) \text{ }y(t)}\Bigg[\frac{tan (\pi/4)}{y(t)} \Bigg]^2 \frac{\delta y(t)}{\delta t} = - (1.1068e-4) \sqrt{ y(t)}\Bigg[\frac{1}{y(t)} \Bigg]^2 \frac{\delta y(t)}{\delta t}= - (1.1068e-4) \frac{\sqrt{ y(t)}}{y(t)^2}

literal b

se inicia el cálculo del siguiente punto de la tabla

i	t	y
0	0	0.15
1	0.5	0.1490
2	1	0.1480
3	1.5	0.1471

i = 0

K_1 = h\Bigg(- (1.1068e-4) \frac{\sqrt{ y(t)}}{y(t)^2} \Bigg) K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.15}}{0.15^2}\Bigg) = -9.5258e-04 K_2 = h\Bigg(- (1.1068e-4) \frac{\sqrt{ y(t)+K_1}}{(y(t)+K_1)^2} \Bigg) K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.15+-9.5258e-04}}{(0.15-9.5258e-04)^2} \Bigg) K_2 = -9.6173e-04 y_1 = y_0 + \frac{K_1 + K_2}{2} y_1 = 0.15 + \frac{-9.5258e-04 -9.6173e-04}{2} = 0.149

i = 1

K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.149}}{0.149^2}\Bigg) =-9.6177e-04 K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.149-9.7120e-04}}{(0.149-9.7120e-04)^2} \Bigg) K_2 = -9.7116e-04 y_2 = y_1 + \frac{-9.6177e-04 + -9.7116e-04}{2} = 0.1480

i = 2

K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.1480}}{0.1480^2}\Bigg) = -9.7120e-04 K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.1480-9.7120e-04}}{(0.1480-9.7120e-04)^2} \Bigg) K_2= -9.8084e-04 y_3 = y_2 + \frac{-9.7120e-04 + -9.8084e-04}{2} = 0.1471

---

literal c

Resultados usando Algoritmo, se encuentra que el embudo se vacia entre 3.15 y 3.20 segundos

 [ t , y , K1 , K2 ]
[[ 0.0000e+00  1.5000e-01  0.0000e+00  0.0000e+00]
 [ 5.0000e-01  1.4904e-01 -9.5258e-04 -9.6173e-04]
 [ 1.0000e+00  1.4808e-01 -9.6177e-04 -9.7116e-04]
 [ 1.5000e+00  1.4710e-01 -9.7120e-04 -9.8084e-04]
 [ 2.0000e+00  1.4611e-01 -9.8088e-04 -9.9078e-04]
 [ 2.5000e+00  1.4512e-01 -9.9083e-04 -1.0010e-03]
...
[ 3.1000e+01  2.8617e-02 -7.5631e-03 -1.0583e-02]
 [ 3.1500e+01  1.0620e-02 -1.1431e-02 -2.4563e-02]
 [ 3.2000e+01         nan -5.0566e-02         nan]
 [ 3.2500e+01         nan         nan         nan]

Instrucciones Python

# 2Eva_2021PAOII_T2 EDO – Embudos cónicos para llenar botellas
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,K1,K2]
    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] = [xi,yi,K1,K2]
    return(estimado)

# INGRESO
d = 0.01
theta = np.pi/4
g = 9.8
d1y = lambda t,y: -(d**2)/4*np.sqrt(2*g*y)*(np.tan(theta)/y)**2

t0 = 0
y0 = 0.15
h  = 0.5
muestras = 70

# PROCEDIMIENTO
tabla = rungekutta2(d1y,t0,y0,h,muestras)

# SALIDA
np.set_printoptions(precision=4)
print('[ t , y , K1 , K2 ]')
print(tabla)

Ejercicio: 2Eva_2021PAOII_T1 Promedio de precipitación de lluvia en área

Los datos de la tabla corresponden a las medidas de precipitación lluviosa en cada cuadrícula. El área observada es de 300×250, los tramos horizontales 6 y los verticales 4.

Los tamaños de paso son:
h = 300/6 = 50
k = 250/4 = 62.5

f(xi,yj)

i \ j	1	2	3	4	5	6
1	0.02	0.36	0.82	0.65	1.7	1.52
2	3.15	3.57	6.25	5	3.88	1.8
3	0.98	0.98	2.4	1.83	0.04	0.01
4	0.4	0.04	0.03	0.03	0.01	0.08

Para obtener el integral doble, usando la forma compuesta de Simpson, primero desarrolla por filas.

i=1

I_1 = \frac{3}{8}(50) ( 0.02+3(0.36)+3(0.82)+0.65) + \frac{1}{3}(50) (0.65+4(1.7)+1.52) = 228.43

i=2

I_2 = \frac{3}{8}(50) ( 3.15+3(3.57)+3(6.25)+5) + \frac{1}{3}(50) (5+4(3.88)+1.8) = 1077.18

i=3

I_3 = \frac{3}{8}(50) ( 0.98+3(0.98)+3(2.4)+1.83) + \frac{1}{3}(50) (1.83+4(0.04)+0.01) = 276.14

i=4

I_4 = \frac{3}{8}(50) ( 0.4+3(0.04)+3(0.03)+0.03) + \frac{1}{3}(50) (0.03+4(0.01)+0.08) = 14.5

con los resultados, luego se desarrolla por columnas:

columnas = [ 228.43,  1077.18,  276.14,  14.5 ]

I_{total} = \frac{3}{8}(62.5) ( 228.43+3(1077.18)+3(276.14)+14.5) I_{total} = 100850.09

para el promedio se divide el integral total para el área de rectángulo.

Precipitacion_{promedio} = \frac{100850.09}{(300)(250)} = 1.34

---

Solución con algoritmo

Resultados usando algoritmos con Python:

integral por fila: [ 228.4375     1077.1875      276.14583333   14.5       ]
integral total: 100850.09765625
promedio precipitación:  1.34466796875
>>>

Instrucciones con Python

# 2Eva_2021PAOII_T1 Promedio de precipitación de 
#  lluvia en área
import numpy as np

def simpson_compuesto(vector,h):
    m = len(vector)
    suma = 0
    j = 0
    while (j+3)<m:  # Simpson 3/8
        f = vector[j:j+4]
        s = (3*h/8)*(f[0]+3*f[1]+3*f[2]+f[3])
        suma = suma + s
        j = j + 3
    while (j+2)<m: # Simpson 1/3
        f = vector[j:j+3]
        s = (h/3)*(f[0]+4*f[1]+f[2])
        suma = suma + s
        j = j + 2
    while (j+1)<m: # Trapecio
        f = vector[j:j+2]
        s = (h/2)*(f[0]+f[1])
        suma = suma + s
        j = j + 1
    return(suma)
    

# INGRESO
A = [[0.02, 0.36, 0.82, 0.65, 1.7 , 1.52],
     [3.15, 3.57, 6.25, 5.  , 3.88, 1.8 ],
     [0.98, 0.98, 2.4 , 1.83, 0.04, 0.01],
     [0.4 , 0.04, 0.03, 0.03, 0.01, 0.08]]
base = 300
altura = 250

h = base/6
k = altura/4

# PROCEDIMIENTO
A = np.array(A)
[n,m] = np.shape(A)
columna = np.zeros(n,dtype=float)
for i in range(0,n,1):
    vector = A[i]
    integralfila = simpson_compuesto(vector,h)
    columna[i] = integralfila

total = simpson_compuesto(columna,k)
promedio = total/(base*altura)
    
# SALIDA
print('integral por fila:', columna)
print('integral total:', total)
print('promedio precipitación: ', promedio)

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)

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.

Según 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