Ejercicio: 3Eva_2024PAOII_T4 Sin()Cos() Integrar con Cuadratura de Gauss

A= \int_0^7 \Big( \sin(0.1t) \cos(0.7t) +3.7 \Big) \delta t

literal a

Para el planteamiento del integral es necesario observar la gráfica de la función a integrar dentro del intervalo.

La función tiene al menos dos "picos" y dos valles en el intervalo. Por lo que un solo tramo del integral podría aumentar el error de integración numérica con una figura trapezoidal equivalente como propone la cuadratura de Gauss.

Se plantea usar al menos dos tramos, y comparar el resultado con tres tramos para observar el error.

Para dos tramos se dispone de los segmentos entre los puntos
[0, 3.5, 7]

Para tres tramos se tiene los segmentos entre los puntos
[ 0, 7/3, 2(7/3), 7]

literal b

Si se usan dos tramos se tienen los segmentos entre los puntos [0,3.5,7]

tramo = [0, 3.5]

x_a = \frac{3.5+0}{2} - \frac{3.5-0}{2}\left(\frac{1}{\sqrt{3}} \right) = 0.7396 x_b = \frac{3.5+0}{2} + \frac{3.5-0}{2}\left(\frac{1}{\sqrt{3}} \right) = 2.7603 f(0.7396) =\sin(0.1(0.7396)) \cos(0.7(0.7396)) +3.7 =3.7642 f(2.7603) =\sin(0.1(2.7603)) \cos(0.7(2.7603)) +3.7 =3.6036 I \cong \frac{3.5-0}{2}(f(0.7396) + f(2.7603)) I \cong \frac{3.5-0}{2}(3.7642 + 3.6036) = 12.8937

tramo = [3.5, 7]

x_a = \frac{3.5+7}{2} - \frac{7-3.5}{2}\left(\frac{1}{\sqrt{3}} \right) = 4.2396 x_b = \frac{3.5+7}{2} + \frac{7-3.5}{2}\left(\frac{1}{\sqrt{3}} \right) = 6.2603 f(4.2396) =\sin(0.1(4.2396)) \cos(0.7(4.2396)) +3.7 =3.2948 f(6.2603) =\sin(0.1(6.2603)) \cos(0.7(6.2603)) +3.7 =3.5100 I \cong \frac{7-3.5}{2}(f(4.2396) + f(6.2603)) I \cong \frac{7-3.5}{2}(3.2948 + 3.5100) = 11.9085

literal c

Integral total : = 12.8937 + 11.9085 = 24.8022

Si de compara con 3 tramos, el error se estima como la diferencia entre los dos integrales calculados

[xa,xb,f(xa),f(xb)]
[0.49309135261210324, 1.8402419807212302, 3.7463820813248043, 3.75103137375189] 8.746982364256144
[xa,xb,f(xa),f(xb)]
[2.8264246859454367, 4.173575314054563, 3.5894184973574266, 3.304431611500099] 8.042825127000448
[xa,xb,f(xa),f(xb)]
[5.159758019278771, 6.506908647387897, 3.2601677912890605, 3.6049588227871885] 8.009314383088956
Integral:  24.79912187434555

Error usando 2 y 3 tramos, es del orden 10^(-3) :

>>> 24.802242263095337 - 24.79912187434555
0.003120388749788816

gráfica con dos tramos:

Ejercicio: 3Eva_2024PAOII_T3 EDO efecto Allee en poblaciones pequeñas

literal a

Dada la ecuación diferencial ordinaria, se reemplazan los valores de las constantes:

\frac{dx}{dt} = rx \left(\frac{x}{K}-1 \right) \left(1-\frac{x}{A} \right) \frac{dx}{dt} = 0.7 x \left(\frac{x}{10}-1 \right) \left(1-\frac{x}{50} \right)

siendo h = 0.2

K_1 = 0.2 f(t_i,x_i) = 0.2\left(0.7 x \left(\frac{x}{10}-1 \right) \left(1-\frac{x}{50} \right) \right) K_2 = 0.2 f(t_i+0.2, x_i + K_1) = 0.2\left(0.7(x+K_1) \left(\frac{x+K_1}{10}-1 \right) \left(1-\frac{x+K_1}{50} \right) \right) x_{i+1} = x_i + \frac{K_1+K_2}{2} t_{i+1} = t_i + 0.2

literal b

Para el desarrollo de las iteraciones se requieren valores iniciales. Para la variable independiente tiempo podría usar t = 0.

Para la variable dependiente población, según la descripción se encontraría entre el intervalo [10,50]. Si x₀<10 la población se extingue. Si la variable x₀>50 se encuentra saturada la capacidad del medio.

Por lo que se propone usar un valor mayor que el mínimo, por ejemplo x₀=11 y otro valor que se encuentre en el intervalo.

itera = 0 , t₀ = 0, x₀ = 11

K_1 = 0.2\left(0.7 (11) \left(\frac{11}{10}-1 \right) \left(1-\frac{11}{50} \right) \right) = 0.1201 K_2 = 0.2\left(0.7(11+0.1201) \left(\frac{11+0.1201}{10}-1 \right) \left(1-\frac{11+0.1201}{50} \right) \right) K_2= 0.1355 x_{i+1} = 11 + \frac{0.1201+0.1355}{2} = 11.1278 t_{i+1} = 0 + 0.2 = 0.2

itera = 1 , t₀ = 0.2, x₀ = 11.1278

K_1 = 0.2\left(0.7 (11.1278) \left(\frac{11.1278}{10}-1 \right) \left(1-\frac{11.1278}{50} \right) \right) = 0.1366 K_2 = 0.2\left(0.7(11.1278+0.1366) \left(\frac{11.1278+0.1366}{10}-1 \right) \left(1-\frac{11.1278+0.1366}{50} \right) \right) K_2= 0.1544 x_{i+1} = 11.1278 + \frac{0.1366+0.1544}{2} =11.2734 t_{i+1} = 0.2 + 0.2 = 0.4

itera = 2 , t₀ = 0.4, x₀ = 11.2734

K_1 = 0.2\left(0.7 (11.2734) \left(\frac{11.2734}{10}-1 \right) \left(1-\frac{11.2734}{50} \right) \right) = 0.1556 K_2 = 0.2\left(0.7(11.2734+0.1556) \left(\frac{11.2734+0.1556}{10}-1 \right) \left(1-\frac{11.2734+0.1556}{50} \right) \right) K_2 = 0.1763 x_{i+1} = 11.2734 + \frac{0.1556+0.1763}{2} = 11.4394 t_{i+1} = 0.4 + 0.2 = 0.6

literal c

Según los resultados de las tres iteraciones anteriores, la población crece. El resultado es acorde al concepto descrito en el enunciado para poblaciones mayores al valor mínimo de extinción. En el ejercicio se usó 11 como valor inicial.

Esto se puede comprobar usando el algoritmo y teniendo los resultados presentados en el literal d

literal d

Los resultados tabulados con el algoritmo son:

EDO con Runge-Kutta 2 Orden
 [ti, xi, K1, K2]
[[ 0.         11.          0.          0.        ]
 [ 0.2        11.12785958  0.12012     0.13559915]
 [ 0.4        11.2734037   0.13660392  0.15448432]
 [ 0.6        11.43943235  0.15566412  0.17639319]
 [ 0.8        11.629282    0.17778585  0.20191345]
 [ 1.         11.84695218  0.20356688  0.23177349]
...
[ 5.6        48.42689825  1.26594886  0.67489419]
 [ 5.8        49.04060051  0.81966583  0.4077387 ]
 [ 6.         49.41989811  0.5143157   0.2442795 ]]

La gráfica del ejercicio para 30 muestras es:

Instrucciones en Python: EDO dy/dx, Runge-Kutta con Python

# 3Eva_2024PAOII_T3 EDO efecto Allee en poblaciones pequeñas
# EDO. Método de RungeKutta 2do Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0, entrega tabla[xi,yi,K1,K2]
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras,
                vertabla=False,precision=6):
    '''solucion a EDO dy/dx, con Runge Kutta de 2do orden
    d1y es la expresi n dy/dx, tambien planteada como f(x,y),
    valores iniciales: x0,y0, tama o de paso h.
    muestras es la cantidad de puntos a calcular. 
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,2+2),dtype=float)
    tabla[0] = [x0,y0,0,0] # incluye el punto [x0,y0]
    
    xi = x0 # valores iniciales
    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
        
        tabla[i] = [xi,yi,K1,K2]
       
    if vertabla==True:
        np.set_printoptions(precision)
        print( 'EDO dy/dx con Runge-Kutta 2 Orden')
        print('i, [xi,     yi,     K1,    K2]')
        for i in range(0,tamano,1):
            print(i,tabla[i])

    return(tabla)

# PROGRAMA PRUEBA -------------------
# INGRESO
# d1y = y' = f, d2y = y'' = f'
r = 0.7
K = 10
A = 50
d1x = lambda t,x: r*x*(x/A-1)*(1-x/K)
t0 = 0
x0 = 11
h  = 0.2
muestras = 30

# PROCEDIMIENTO
tabla = rungekutta2(d1x,t0,x0,h,muestras)
xi = tabla[:,0]
yiRK2 = tabla[:,1]

# SALIDA
print( 'EDO con Runge-Kutta 2 Orden')
print(' [ti, xi, K1, K2]')
print(tabla)

# GRAFICA
import matplotlib.pyplot as plt
plt.plot(xi,yiRK2)
plt.plot(xi[0],yiRK2[0],
         'o',color='r', label ='[t0,x0]')
plt.plot(xi[1:],yiRK2[1:], 'o',color='m',
         label ='[ti,xi] Runge-Kutta 2 Orden')

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

 

Ejercicio: 3Eva_2024PAOII_T2 Interpolar trayectoria helicóptero

xi = [0.  , 2.50, 3.75, 5.00, 6.25, 7.5 ]
yi = [3.7 , 3.25, 4.05, 4.33, 2.95, 3.22]

literal a

Según la gráfica presentada, los puntos a considerar para todo el intervalo, deberían ser al menos el primero y el último. Para un polinomio de grado 3 se requieren usar 4 muestras, por lo que faltarían dos muestras dentro del intervalo. Los puntos adicionales estarían entre los puntos intermedios, tratando de mantener una distancia equidistante entre ellos y tratar de mantener los cambios de dirección.

Los puntos seleccionados para el ejercicio serán

xi = [0. , 2.50, 5.00, 7.5 ]
fi = [3.7 , 3.25, 4.33, 3.22]

Se puede construir el polinomio con los cualquiera de los métodos para interpolación dado que tienen tamaños de paso iguales entre tramos.

Desarrollando con el método de Lagrange, el polinomio se construye como:

término 1

L_{0} (x) = \frac{(x-2.5)(x-5.0)(x-7.5)}{(0-2.5)(0-5.0)(0-7.5)}

término 2

L_{1} (x) = \frac{(x-0)(x-5.0)(x-7.5)}{(2.5-0)(2.5-5.0)(2.5-7.5)}

término 3

L_{2} (x) = \frac{(x-0)(x-2.5)(x-7.5)}{(5-0)(5-2.5)(5-7.5)}

término 4

L_{3} (x) = \frac{(x-0)(x-2.5)(x-5)}{(7.5-0)(7.5-2.5)(7.5-5)}

se construye el polinomio usando la fórmula para f_n(x) para cada valor fi,

p_3(x) = 3.7 \frac{(x-2.5)(x-5.0)(x-7.5)}{(0-2.5)(0-5.0)(0-7.5)} + 3.25 \frac{(x-0)(x-5.0)(x-7.5)}{(2.5-0)(2.5-5.0)(2.5-7.5)} + 4.33 \frac{(x-0)(x-2.5)(x-7.5)}{(5-0)(5-2.5)(5-7.5)} + 3.22 \frac{(x-0)(x-2.5)(x-5)}{(7.5-0)(7.5-2.5)(7.5-5)}

simplificando con Sympy:

Polinomio de Lagrange, expresiones
0.104*x*(x - 7.5)*(x - 5.0) - 0.13856*x*(x - 7.5)*(x - 2.5) + 0.0343466666666667*x*(x - 5.0)*(x - 2.5) - 0.0394666666666667*(x - 7.5)*(x - 5.0)*(x - 2.5)

Polinomio de Lagrange: 
-0.03968*x**3 + 0.42*x**2 - 0.982*x + 3.7

p_3(x) = 3.7 - 0.982x + 0.42 x^2 -0.03968 x^3

literal b

Para verificar que el polinomio pasa por los puntos, se puede usar una gráfica o al menos dos puntos usados para crear el polinomio:

p_3(2.5) = 3.7 - 0.982(2.5) + 0.42 (2.5)^2 -0.03968 (2.5)^3 = 3.25 p_3(5) = 3.7 - 0.982(5) + 0.42 (5)^2 -0.03968 (5)^3 = 4.33

>>> polisimple
-0.03968*x**3 + 0.42*x**2 - 0.982*x + 3.7
>>> polisimple.subs(x,2.5)
3.25000000000000
>>> polisimple.subs(x,5)
4.33000000000000
>>>

Se comprueba que los valores obtenidos corresponden a las muestras, por lo que el polinomio cumple con los criterios básicos de interpolación. La gráfica permite verificar también el resultado.

literal c

Error en puntos no usados de las muestras

p_3(3.75) = 3.7 - 0.982(3.75) + 0.42 (3.75)^2 -0.03968 (3.75)^3 = 3.83125

error = |3.83125 - 4.05| = 0.218

p_3(6.25) = 3.7 - 0.982(6.25) + 0.42 (6.25)^2 -0.03968 (6.25)^3 = 4.28125

error = |4.28125 - 4.05| = 1.3312

literal d

Observando las gráficas de la trayectoria construida junto a las funciones descritas en el tema 1 se tiene que:

Un polinomio de grado 3 es insuficiente para describir la trayectoria, se debe aumentar el grado del polinomio para ajustar mejor la curva.

Por ejemplo, usando todos los puntos, la trayectoria y el polinomio son mas cercanas aunque no iguales.

literal e

Encontrar el error en P(5), como x=5 y es parte de los puntos de muestra, el error debería ser cero. Siempre y cuando x=5 sea parte de los puntos seleccionados.

p_3(5) = 3.7 - 0.982(5) + 0.42 (5)^2 -0.03968 (5)^3 = 4.33

error = | 4.33-4.33| = 0

literal f

Los resultados con el algoritmo de Lagrange se muestran como:

    valores de fi:  [3.7  3.25 4.33 3.22]
divisores en L(i):  [-93.75  31.25 -31.25  93.75]

Polinomio de Lagrange, expresiones
0.104*x*(x - 7.5)*(x - 5.0) - 0.13856*x*(x - 7.5)*(x - 2.5) + 0.0343466666666667*x*(x - 5.0)*(x - 2.5) - 0.0394666666666667*(x - 7.5)*(x - 5.0)*(x - 2.5)

Polinomio de Lagrange: 
-0.03968*x**3 + 0.42*x**2 - 0.982*x + 3.7

Algoritmo en Python. Interpolación polinómica de Lagrange

# 3Eva_2024PAOII_T2 Interpolar trayectoria helicóptero
# Interpolacion de Lagrange
# divisoresL solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
# todos los datos
xj = [0.  , 2.50, 3.75, 5.00, 6.25, 7.5 ]
fj = [3.7 , 3.25, 4.05, 4.33, 2.95, 3.22]

# datos seleccionados
xi = [0.  , 2.50, 5.00, 7.5 ]
fi = [3.7 , 3.25, 4.33, 3.22]

# trayectoria
gx = lambda t: 0.5*t +0 *t
gy = lambda t: np.sin(0.10*t)*np.cos(0.7*t)+ 3.7

ta = 0 ; tb = 15
muestrast = 7
muestrasj = 4*muestrast

# PROCEDIMIENTO
# trayectoria
tj = np.linspace(ta,tb,muestrasj)
gxj = gx(tj)
gyj = gy(tj)

# Interpolación
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)
# 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)

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

# Gráfica
plt.plot(gxj,gyj,color='orange', label = 'trayectoria')
plt.plot(xi,fi,'o', label = 'muestras')
plt.plot(pxi,pfi,color='green',linestyle='dashed', label = 'P(x)')
plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.grid()
plt.title('Interpolación Lagrange')
plt.show()

Ejercicio: 3Eva_2024PAOII_T1 Accidente entre aeronaves

literal a

La distancia entre las aeronaves se determina a partir de las ecuaciones proporcionadas en el enunciado.

d = \sqrt{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2}

La ecuación también mantiene la forma y el mínimo si se utiliza el cuadrado de la distancia:

D = d^2 = (x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2

Las diferencias de distancia por eje pueden ser convergentes hacia el punto de impacto, dado que se conoce que el accidente ya ocurrió. Por lo que también se pueden revisar las distancias entre ejes en lugar de la distancia total en 3D, simplificando un poco el ejercicio. Observe la gráfica proporcionada:

Avión	Helicóptero	distancia por eje
Ax(t) = 5.1	Hx(t) = 0.5t	|0.5t-5.1|
Ay(t) = 0.4t	Hy(t) = sin(0.1t)cos(0.7t)+3.7	|sin(0.1t)cos(0.7t) + 3.7-0.4|
Az(t) = 0.5 e^{-0.2t} + 0.3	Hz(t) = 0.36	|0.36 - (0.5 e^{-0.2t} + 0.3)|

D = (0.5t-5.1)^2 + (sin(0.1t)cos(0.7t)+3.7-0.4t)^2 + (0.36-(0.5 e^{-0.2t} + 0.3))^2 dx =0.5t-5.1 dy = sin(0.1t)cos(0.7t) + 3.7-0.4 dz = 0.36 - (0.5 e^{-0.2t} + 0.3)

Se realiza la gráfica para la distancia al cuadrado Di y las distancias por cada eje dxi, dyi, dzi :

Por lo que el resultado también se podría determinar usando por ejemplo el eje y. La ecuación a buscar la distancia mínima, punto de choque o cruce por cero podría ser:

f(t) =dy(t)= sin(0.1t)cos(0.7t) + 3.7-0.4t = 0

literal b

Por la gráfica se puede obtener un intervalo de búsqueda entre [8, 12], que usando solo las distancias en el eje y, se tiene:

dy(8) =1.0563 dy(12) =-1.5839

literal c

Se pide usar un método de búsqueda de raíces, pudiendo seleccionar Bisección en el intervalo encontrado en el literal b.

iteración 1

a = 8, b=12 c = \frac{a+b}{2} = \frac{8+12}{2} = 10 f(8) = 1.0563 f(10) =sin(0.1(10))cos(0.7(10)) + 3.7-0.4(10) = 0.3343 f(12) = -1.5839

cambio de signo a la derecha

a = 10, b = 12 tramo = |12-10| =2

iteración 2

a = 10, b=12 c = \frac{10+12}{2} = 11 f(11) =sin(0.1(11))cos(0.7(11)) + 3.7-0.4(11) = -0.5633

cambio de signo a la izquierda

a = 10, b = 11 tramo = |11-10| = 1

iteración 3

a = 10, b=11 c = \frac{10+11}{2} = 10.5 f(10.5) =sin(0.1(10.5))cos(0.7(10.5)) + 3.7-0.4(10.5) = -0.0811

cambio de signo a la izquierda

a = 10, b = 10.5 tramo = |10.5-10| = 0.5

literal d

La variable independiente en el ejercicio es tiempo 't', que podría ser segundos. Si la tolerancia se estima en milisegundos, tolera = 10-3 .

El error se ha calculado en cada iteración

literal e

El error disminuye en cada iteración, por lo que el método converge.

La distancia mínima entre las aeronaves no tiene que llegar a cero para que se produzca un accidente. Las dimensiones de las aeronaves muestran que se encuentran entre 70m y 4 m, por lo que se estima que debe existir una separación mayor a 70 metros en cualquiera de los ejes para evitar un accidente. Si las coordenadas se estiman en Km, la tolerancia sería de 0.070 Km al considerar más de 70 metros como la distancia segura.

literal f

Usando el algoritmo se encuentra la raíz en:

i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
   tramo
0 [8, 10.0, 12] [ 1.0564  0.3344 -1.584 ]
   2.0
1 [10.0, 11.0, 12] [ 0.3344 -0.5633 -1.584 ]
   1.0
2 [10.0, 10.5, 11.0] [ 0.3344 -0.0811 -0.5633]
   0.5
3 [10.0, 10.25, 10.5] [ 0.3344  0.1368 -0.0811]
   0.25
4 [10.25, 10.375, 10.5] [ 0.1368  0.0302 -0.0811]
   0.125
5 [10.375, 10.4375, 10.5] [ 0.0302 -0.0249 -0.0811]
   0.0625
6 [10.375, 10.40625, 10.4375] [ 0.0302  0.0028 -0.0249]
   0.03125
7 [10.40625, 10.421875, 10.4375] [ 0.0028 -0.011  -0.0249]
   0.015625
8 [10.40625, 10.4140625, 10.421875] [ 0.0028 -0.0041 -0.011 ]
   0.0078125
9 [10.40625, 10.41015625, 10.4140625] [ 0.0028 -0.0007 -0.0041]
   0.00390625
10 [10.40625, 10.408203125, 10.41015625] [ 0.0028  0.001  -0.0007]
   0.001953125
11 [10.408203125, 10.4091796875, 10.41015625] [ 0.001   0.0002 -0.0007]
   0.0009765625
raíz en:  10.4091796875

Por lo que las coordenadas de choque entre aeronaves será:

Avión	Helicóptero	distancia por eje
Ax(t) = 5.1	Hx(t) = 0.5 (10.40) = 5.2	|5.1-5.2| = 0.1
Ay(t) = 0.4(10.4) = 4.16	Hy(t) = sin(0.1(10.40))cos(0.7(10.40))+3.7 = 4.168	|4.16 - 4.168| = 0.008
Az(10.4) = 0.5 e^{-0.2(10.4)} + 0.3 = 0.3624	Hz(t) = 0.36	|0.3624 - 0.36|  = 0.0024

Instrucciones en Python para Gráfica de distancias

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO
# Avión Aterriza
Ax = lambda t: 5.1 + 0*t
Ay = lambda t: 0.4*t
Az = lambda t: 0.5*np.exp(-0.2*t) + 0.3
# helicóptero
Hx = lambda t: 0.5*t + 0*t
Hy = lambda t: np.sin(0.10*t)*np.cos(0.7*t)+ 3.7
Hz = lambda t: 0.36 + 0*t

a = 0
b = 6/0.5 # x entre[0,6] usando gx(t)= 6
muestras = 41

# PROCEDIMIENTO
# Distancia por ejes
dx = lambda t: Hx(t) - Ax(t)
dy = lambda t: Hy(t) - Ay(t)
dz = lambda t: Hz(t) - Az(t)
# Distancia 3D
distancia2 = lambda t: dx(t)**2 + dy(t)**2 + dz(t)**2

# Simulacion en segmento t=[a,b]
ti = np.linspace(a,b,muestras)
dxi = dx(ti)
dyi = dy(ti)
dzi = dz(ti)

Di = distancia2(ti)

# SALIDA
print(dy(8))
print(dy(12))

plt.plot(ti,Di, label='Di')
plt.plot(ti,dxi,label='dxi')
plt.plot(ti,dyi,label='dyi')
plt.plot(ti,dzi,label='dzi')
plt.xlabel('ti')
plt.ylabel('distancia2')
plt.grid()
plt.legend()
plt.show()

Instrucciones en Python - Algoritmo Bisección

# 3Eva_2024PAOII_T1 Accidente entre aeronaves
# 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

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

def biseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
    '''
    Algoritmo de Bisección
    Los valores de [a,b] son seleccionados
    desde la gráfica de la función
    error = tolera
    '''
    fa = fx(a)
    fb = fx(b)
    tramo = np.abs(b-a)
    itera = 0
    cambia = np.sign(fa)*np.sign(fb)
    if cambia<0: # existe cambio de signo f(a) vs f(b)
        if vertabla==True:
            print('método de Bisección')
            print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
            print('  ','tramo')
            np.set_printoptions(precision)
            
        while (tramo>=tolera and itera<=iteramax):
            c = (a+b)/2
            fc = fx(c)
            cambia = np.sign(fa)*np.sign(fc)
            if vertabla==True:
                print(itera,[a,c,b],np.array([fa,fc,fb]))
            if (cambia<0):
                b = c
                fb = fc
            else:
                a = c
                fa = fc
            tramo = np.abs(b-a)
            if vertabla==True:
                print('  ',tramo)
            itera = itera + 1
        respuesta = c
        # Valida respuesta
        if (itera>=iteramax):
            respuesta = np.nan

    else: 
        print(' No existe cambio de signo entre f(a) y f(b)')
        print(' f(a) =',fa,',  f(b) =',fb) 
        respuesta=np.nan
    return(respuesta)

# INGRESO
fx = lambda t: np.sin(0.10*t)*np.cos(0.7*t)+ 3.7 - 0.4*t 
a = 8
b = 12
tolera = 0.001

# PROCEDIMIENTO
respuesta = biseccion(fx,a,b,tolera,vertabla=True)
# SALIDA
print('raíz en: ', respuesta)

Ejercicio: 3Eva_2023PAOII_T3 Volumen por solido de revolución de un peón

El volumen se calcula a partir de la expresión:

V = \int_{a}^{b} \pi (f(x))^2 dx

xi=[ 0, 3, 5.  , 9.985 , 14.97 , 17.97, 40.04, 43.29, 51.6449, 60]
yi=[15,15,13.25,14.1552, 9.6768,  9.67,  4.64,  4.64,  8.9768, 0.]

Para el intervalo [0,3]

La función es una constante

literal a, b y c

f(x) = 15 V = \int_{0}^{3} \pi (15)^2 dx = \pi (15)^2 x \big|_{0}^{3} = \pi (15)^2 (3-0) = 2120.5750

Para el intervalo [3,5]

literal a

La función es una recta con pendiente negativa

f(x) = -0.875 x + 17.625 V = \int_{3}^{5} \pi (f(x))^2 dx V = \int_{3}^{5} \pi (-0.875 x+17.625)^2 dx

para el integral con cuadratura de Gauss

g(x) = \pi (-0.875*x+17.625)^2

literal b y c

x_a = \frac{5+3}{2} - \frac{5-3}{2}\frac{1}{\sqrt{3}} = 3.4226 x_b = \frac{5+3}{2} + \frac{5-3}{2}\frac{1}{\sqrt{3}} = 4.5773

con lo que el resultado aproximado del integral se convierte en:

I \cong \frac{5-3}{2}(g(3.4226) + g(4.5773)) = \frac{5-3}{2} (672.43+582.76) = 1255.1971

Para el intervalo [5, 14.97]

literal a

Del polinomio obtenido en el ejercicio anterior para éste intervalo

f(x) = -0.1083 x^2+1.8047 x + 6.9341 V = \int_{5}^{14.97} \pi (f(x))^2 dx V = \int_{5}^{14.97} \pi (-0.1083 x^2+1.8047 x + 6.9341)^2 dx [ g(x) = \pi (-0.1083 x^2+1.8047 x + 6.9341)^2

literal b y c

x_a = \frac{14.97+5}{2} - \frac{14.97-5}{2}\frac{1}{\sqrt{3}} = 7.1069 x_b = \frac{14.97+5}{2} + \frac{14.97-5}{2}\frac{1}{\sqrt{3}} = 12.8630

con lo que el resultado aproximado del integral se convierte en:

I \cong \frac{14.97-5}{2}(g(7.1069) + g(12.8630)) = \frac{14.97-5}{2} (641.5176+469.8124) = 5539.9805

literal d

El volumen total es la suma de los resultados de cada una de las secciones:

V_{total} = 2120.5750 + 1255.1971 + 5539.9805 + ...

Tarea: continuar con los otros intervalos para obtener el volumen total.

literal e

Resultados con el algoritmo

para f(x):
[xa,xb]= [0.6339745962155612, 2.366025403784439]
[f(xa),f(xb)]= [915.4418590078262, 760.106947830205]
Volumenfx:  2513.323210257047

para f(x):
[xa,xb]= [3.4226497308103743, 4.577350269189626]
[f(xa),f(xb)]= [672.4334354361756, 582.7637293668464]
Volumenfx:  1255.1971648030221

para f(x):
[xa,xb]= [7.106908908089714, 12.863091091910285]
[f(xa),f(xb)]= [641.5176069162055, 469.8124936184865]
Volumenfx:  5539.98055116544

Instrucciones en Python

# 3Eva_2023PAOII_T3 Volumen por solido de revolución de un peón
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# INGRESO
# intervalos eje x, funciones
#a=0 ; b=3 ; f = lambda x: np.pi*(15)**2 + 0*x
#a=3 ; b=5 ; f = lambda x: np.pi*(-0.875*x+17.625)**2
# con el polinomio del tema 2
a=5 ; b=14.97 ; f = lambda x: np.pi*(-0.1083*x**2+1.8047*x + 6.9341)**2
# con el circulo del tema 1
#a=5 ; b=14.97 ; f = lambda x: np.pi*(np.sqrt(6.7**2-(x-8.6)**2) + 7.6)**2

xmuestras = 31

# PROCEDIMIENTO
# Volumen para f(x) con Cuadratura de Gauss
xa = (b+a)/2 + (b-a)/2*(-1/np.sqrt(3)) 
xb = (b+a)/2 + (b-a)/2*(1/np.sqrt(3))
Volumenfx = (b-a)/2*(f(xa)+f(xb))
xab = [xa,xb]
fab = [f(xa),f(xb)]

# SALIDA
print("para f(x):")
print("[xa,xb]=", xab)
print("[f(xa),f(xb)]=", fab)
print("Volumenfx: ",Volumenfx)
print()

Añadir para graficar en 2D

# grafica 2D para el area de corte en eje y,x --------
# muestras en x
xi = np.linspace(a, b, xmuestras)
fi = f(xi)
f0 = np.zeros(xmuestras,dtype=float)

# grafica 2D
figura2D = plt.figure()
grafica2D = figura2D.add_subplot(111)
grafica2D.plot(xi,fi,color='blue',label='f(x)')
grafica2D.fill_between(xi,fi,f0,color='lightblue')
grafica2D.grid()
grafica2D.set_title('Area para sólido de revolución')
grafica2D.set_xlabel('x')
grafica2D.set_ylabel('f(x)')

Añadir para graficar en 3D

# Para grafica 3D -----------------------
# angulo w de rotación
w_a = 0
w_b = 2*np.pi
w_muestras = 31

# grafica 3D muestras en x y angulo w
wi = np.linspace(w_a, w_b, w_muestras)
X, W = np.meshgrid(xi, wi)
# proyeccion en cada eje 
Yf = f(xi)*np.cos(W)
Zf = f(xi)*np.sin(W)

# grafica 3D
figura3D = plt.figure()
grafica = figura3D.add_subplot(111, projection='3d')

grafica.plot_surface(X, Yf, Zf,
color='blue', label='f(x)',
alpha=0.6, rstride=6, cstride=12)

grafica.set_title('Sólido de revolución')
grafica.set_xlabel('x')
grafica.set_ylabel('y')
grafica.set_zlabel('z')
# grafica.legend()
eleva = 30
rota = -45
deltaw = 5
grafica.view_init(eleva, rota)

# rotacion de ejes
for angulo in range(rota, 360+rota, deltaw ):
grafica.view_init(eleva, angulo)
plt.draw()
plt.pause(.001)
plt.show()

.

Ejercicio: 3Eva_2023PAOII_T2 perfil de un peón

literal a y b

Para el planteamiento de los polinomios, es necesario observar la gráfica proporcionada para el ejercicio y la correspondiente tabla de datos:

xi	0	3	5	9.985	14.97	17.97	40.04	43.29	51.6456	60
yi	15	15	13.25	14.155	9.676	9.676	4.64	4.64	8.976	0

Para el intervalo [0,3], El perfil se describe con una constante, por lo tanto:

p_1(x) = 15

Para el intervalo [3,5] el perfil se describe como una recta con pendiente negativa.

p_2(x) =a_0 + a_1 x a_1 = \frac{ 13.25-15}{5-3} = -0.875 p_2(5) = 13.25 =a_0 -0.875 (5) a_0 = 13.25 +0.875 (5) = 17.625 p_2(x) = -0.875 x + 17.625

Para el con los puntos xi = [5, 9.985, 14.97] el perfil se describe con tres puntos, por lo que el polinomio de interpolación puede ser de grado 2. Usando el método de Lagrange:

p_3(x) =\frac{(x-9.985)( x-14.97)}{(5-9.985)( 5-14.97)} 13.25 + +\frac{(x-5)( x-14.97)}{(9.985-5)( 9.985-14.97)} 14.155 + +\frac{(x-5)( x-9.985)}{(14.97-5)( 14.97-9.985)} 9.676 +

simplificando la expresión con Python y Sympy:

P_3(x) = -0.1083 x^2 + 1.8047 x + 6.9341

literal c

El error para los polinomios acorde a la gráfica proporcionada es cero para los tramos [0,3] y [3,5], dado que el primero es una constante y el segundo es una recta con pendiente.

El error de mayor magnitud se presenta con la interpolación entre el círculo de la ecuación del tema 1 y el polinomio de grado 2 en el intervalo [5, 14.97]. Tomando como ejemplo el punto intermedio del tramo derecho en:

x_k = \frac{14.97-5}{4} + 9.985= 12.4775 p_3(12.4775) = -0.1083(12.4775)^2 + 1.8047(12.4775) + 6.9341 = 12.5889 f(12.4775) = \sqrt{6.7^2-(12.4775-8.6)^2} + 7.6 = 13.0639 errado = |p_3(12.4775) - f(12.4775)| = 0.4750

literal d

Se puede mejorar la aproximación del polinomio aumentando el grado y número de puntos a usar dentro del intervalo. El resultado se debería acercar mas al perfil superior del círculo de referencia del tema 1.

Resultados para el intervalo [5, 14.97]

    valores de fi:  [13.25   14.1552  9.6768]
divisores en L(i):  [ 49.70045  -24.850225  49.70045 ]

Polinomio de Lagrange, expresiones
0.266597183727713*(x - 14.97)*(x - 9.985) 
- 0.569620596996607*(x - 14.97)*(x - 5.0) 
+ 0.194702462452553*(x - 9.985)*(x - 5.0)

Polinomio de Lagrange: 
-0.108320950816341*x**2 + 1.80477420224566*x + 6.93415275918024

Instrucciones Python

# 3Eva_2023PAOII_T2 perfil de un peón
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
xj=[ 0, 3, 5.  , 9.985 , 14.97 , 17.97, 40.04, 43.29, 51.6449, 60]
yj=[15,15,13.25,14.1552, 9.6768,  9.67,  4.64,  4.64,  8.9768, 0.]


#  Datos de prueba
ia = 2
ib = 4+1
xi = np.array(xj[ia:ib])
fi = np.array(yj[ia:ib])

g = lambda x: np.sqrt(6.7**2-(x-8.6)**2) + 7.6
muestras = 21

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

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)
gi = g(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)

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

 

Ejercicio: 3Eva_2023PAOII_T1 Intersección con círculo
Encuentre las raíces de las ecuaciones simultaneas siguientes:

literal a y b

Use el enfoque gráfico para obtener los valores iniciales.

2y+1.75 x = 35.25 (y-7.6)^2 + (x-8.6)^2 = (6.7)^2

se despeja la variable dependiente de cada ecuación, para la primera:

f(x) = y = \frac{35.25}{2} - \frac{1.75}{2} x

para la segunda:

(y-7.6)^2 = (6.7)^2 - (x-8.6)^2 g(x) = y = \sqrt{(6.7)^2 - (x-8.6)^2} + 7.6

Al buscar la intersección entre f(x) y g(x) se puede encontrar con la raiz de:

distancia(x) = f(x) - g(x) distancia(x) = \Big( \frac{35.25}{2} - \frac{1.75}{2} x\Big) - \Big(\sqrt{(6.7)^2 - (x-8.6)^2} + 7.6\Big)

La primera ecuación es una recta, por lo que no aporta a las cotas de la gráfica.

La segunda ecuación es la general de un círculo centrado en (7.6, 8.6) y radio 6.7, por lo que se considera el intervalo para la gráfica entre:

[7.6 -6.7, 7.6+6.7] [0.9, 14.3]

Con lo que se puede formar la gráfica de la parte superior del círculo en Python:

Instrucciones en Python

# 3Eva_2023PAOII_T1 Intersección con círculo
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
f = lambda x: -(1.75/2)*x + (35.25/2)
g = lambda x: np.sqrt(6.7**2-(x-8.6)**2) + 7.6
distancia = lambda x: f(x)- g(x)

a = 0.5
b = 16
muestras = 21

# PROCEDIMIENTO
# literal a y b
xi = np.linspace(a,b,muestras)
fi = f(xi)
gi = g(xi)
dist_i = distancia(xi)

# SALIDA - Grafica
# literal a y b
plt.plot(xi,fi, label='f(x)')
plt.plot(xi,gi, label='g(x)')
plt.plot(xi,dist_i, label='f(x)-g(x)')
plt.axhline(0, color='black', linestyle='dashed')
plt.axvline(5, color='red', linestyle='dashed')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

literal c

Un punto inicial de búsqueda dentro del intervalo puede ser x0=3.

Para Newton-Raphson se usa:

f(x) = \Big( \frac{35.25}{2} - \frac{1.75}{2} x\Big) - \Big(\sqrt{(6.7)^2 - (x-8.6)^2} + 7.6\Big) \frac{d}{dx}f(x) = -\frac{8.6(0.1162-0.0135x)}{\sqrt{0.0609-(0.1162x-1)^2}}-0.875

(derivada obtenida con sympy)

itera = 0 ; xi = 3

f(3) = \Big( \frac{35.25}{2} - \frac{1.75}{2} (3)\Big) - \Big(\sqrt{(6.7)^2 - ((3)-8.6)^2} + 7.6\Big) \frac{d}{dx}f(3) = -\frac{8.6(0.1162-0.0135(3))}{\sqrt{0.0609-(0.1162(3)-1)^2}}-0.875 x_1 = x_0 - \frac{f(3)}{\frac{d}{dx}f(3)} = 4.55 tramo = |4.55-3| = 1.55

itera = 1 ; xi = 4.55

f(3) = \Big( \frac{35.25}{2} - \frac{1.75}{2} (4.55)\Big) - \Big(\sqrt{(6.7)^2 - ((4.55)-8.6)^2} + 7.6\Big) \frac{d}{dx}f(3) = -\frac{8.6(0.1162-0.0135(4.55))}{\sqrt{0.0609-(0.1162(4.55)-1)^2}}-0.875 x_2 = x_1 - \frac{f(4.55)}{\frac{d}{dx}f(4.55)} = 4.98 tramo = |4.98-4.55| = 0.43

itera = 2 ; xi = 4.98

f(3) = \Big( \frac{35.25}{2} - \frac{1.75}{2} (4.98)\Big) - \Big(\sqrt{(6.7)^2 - ((4.98)-8.6)^2} + 7.6\Big) \frac{d}{dx}f(3) = -\frac{8.6(0.1162-0.0135(4.98))}{\sqrt{0.0609-(0.1162(4.98)-1)^2}}-0.875 x_3 = x_2 - \frac{f(4.98)}{\frac{d}{dx}f(4.98)} = 4.99 tramo = |4.99-4.98| = 0.01

literal d

El error disminuye en cada iteración, por lo que el método converge.

Con el algoritmo se muestra que converge a x=5

dfg(x) = 
     8.6*(0.116279069767442 - 0.0135208220659816*x)          
- --------------------------------------------------- - 0.875
     ________________________________________________        
    /                                              2         
  \/  0.606949702541915 - (0.116279069767442*x - 1)          

['xi', 'xnuevo', 'tramo']
[[3.0000e+00 4.5524e+00 1.5524e+00]
 [4.5524e+00 4.9825e+00 4.3018e-01]
 [4.9825e+00 4.9995e+00 1.6999e-02]
 [4.9995e+00 4.9996e+00 2.3865e-05]]
raiz en:  4.9995611025201585
con error de:  2.3865455016647275e-05

Instrucciones en Python

# 3Eva_2023PAOII_T1 Intersección con círculo
import sympy as sym
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# forma algebraica con sympy
x = sym.Symbol('x')
f = -(1.75/2)*x + (35.25/2)
g = sym.sqrt(6.7**2-(x-8.6)**2) + 7.6
distancia = f - g

x0 = 3
tolera = 0.0001 # 1e-4

# PROCEDIMIENTO
# literal c
dfg = sym.diff(distancia,x)
# convierte a forma numerica con numpy
# Newton-Raphson
fx = sym.lambdify(x,distancia)
dfx = sym.lambdify(x,dfg)

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

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

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

Ejercicio: 3Eva_2021PAOI_T3 EDO 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: y₀(t)=0 , y’₀(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,t_n=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

t₀ = 0, y₀ = 0,  y’₀ = z₀ = -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()

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:

T_CA = T_CB

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

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 T_CB = 655.32162903
T_CB = 400

luego la primera fila:

0.81915204T_CA -0.25881905T_CB = -103.52761804

0.81915204T_CA = 0.25881905T_CB  -103.52761804

T_CA = 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 θ₁ y θ₂:

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