s1Eva_2023PAOI_T1 Desacople de cohete de dos etapas

Ejercicio: 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas

literal a. Planteamiento

La ecuación a usar según el enunciado es y usando los valores dados es:

v = u \ln\Big(\frac{m_0}{m_0-qt}\Big) - gt 800 = 1870 \ln\Big(\frac{195000}{195000-2500 t}\Big) - 9.8 t

con lo que la función para buscar la raíz es:

f(t) = 1870 \ln\Big(\frac{195000}{195000-2500 t}\Big) - 9.8 t -800

Literal b. Intervalo de búsqueda

Para el intervalo de búsqueda se puede usar una gráfica e interpretar el punto a buscar alrededor de 800 m/s. Que de la gráfica se observa que un intervalo alrededor de 35 sería válido para el método de la Bisección. Para otros métodos abiertos, también es posible deducir un punto t0.

La validación se muestra con la primera iteración al evaluar f(30) que es negativo y f(40) que es de signo positivo.

Desacople de cohete de dos etapasInstrucciones en Python

# 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
v = lambda t: 1870*np.log(195000/(195000-2500*t))-9.8*t
a = 0
b = 40
tramos = 51

# PROCEDIMIENTO
ti = np.linspace(a,b,tramos)
vi = v(ti)

# SALIDA
plt.plot(ti,vi)
plt.xlabel('ti')
plt.ylabel('vi')
plt.title('Velocidad vertical vs tiempo')
plt.grid()
plt.show()

literal c. Desarrollo con algoritmo de Bisección

itera =0, intervalo [30,40]

f(30) = 1870 \ln\Big(\frac{195000}{195000-2500 (30)}\Big) - 9.8 (30) -800 = -186.100 f(40) = 1870 \ln\Big(\frac{195000}{195000-2500 (40)}\Big) - 9.8 (40) -800 = 152.759 c= \frac{a+b}{2}= \frac{30+40}{2} = 35 f(35) = 1870 \ln\Big(\frac{195000}{195000-2500 (35)}\Big) - 9.8 (35) -800 = -29.399

error = tramo = |40-30| = 10

como f(a) y f(c) son del mismo signo, el intervalo nuevo será [35,40]

itera =1 , intervalo [35,40]

c= \frac{35+40}{2} = 37.5 f(37.5) = 1870 \ln\Big(\frac{195000}{195000-2500 (37.5)}\Big) - 9.8 (37.5) -800 = 58.111

error = tramo = |40-35| = 5

como f(c) y f(b) son del mismo signo, el intervalo nuevo será [35,37.5]

itera =2 , intervalo [35,37.5]

c= \frac{37.5+35}{2} = 36.25 f(36.25) = 1870 \ln\Big(\frac{195000}{195000-2500 (36.25)}\Big) - 9.8 (36.25) -800 = 13.518

error = tramo = |37.5-35| = 2.5

como f(c) y f(b) son del mismo signo, el intervalo nuevo será [35,36.25]

[ i, a,    c,     b,     f(a),    f(c),   f(b),   tramo]
  1 30.000 35.000 40.000 -186.100 -29.399 152.759 10.000 
  2 35.000 37.500 40.000 -29.399   58.111 152.759  5.000 
  3 35.000 36.250 37.500 -29.399   13.518 58.111   2.500 
  4 35.000 35.625 36.250 -29.399   -8.144 13.518   1.250 
  5 35.625 35.938 36.250 -8.144     2.635 13.518   0.625 
  6 35.625 35.781 35.938 -8.144    -2.767 2.635    0.312 
  7 35.781 35.859 35.938 -2.767    -0.069 2.635    0.156 
  8 35.859 35.898 35.938 -0.069     1.282 2.635    0.078 
raiz:  35.8984375
>>>

literal d. tolerancia y errores

La tolerancia depende de la escala a la que se mide y el instrumento de medición. Si consideramos décimas de segundo la tolerancia será de 10-1. ó 0.1

Los errores entre iteraciones se muestran en el literal anterior.

literal e. convergencia

Los errores en cada iteración disminuye, lo que muestra que el método converge. Luego de 8 iteraciones se encuentra el tiempo ti a usar como 35.8 s.

Se adjunta el algoritmo de la bisección ajustado para el ejercicio. La gráfica es complementaria a la presentada en el literal b, que puede ser incorporada al mismo algoritmo para la presentación.

Instrucciones en Python

# Algoritmo de Bisección
# 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas
import numpy as np

# INGRESO
fx = lambda t: 1870*np.log(195000/(195000-2500*t))-9.8*t -800
a = 30
b = 40
tolera = 0.1

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

s2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

Ejercicio: 2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

\frac{\partial^2 u}{\partial x^2} = b \frac{\partial u}{\partial t}

2Eva_2022PAOII_T3 EDP Parabolica Malla

\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = b\frac{u_{i,j+1}-u_{i,j}}{\Delta t}

agrupando variables

\frac{\Delta t}{b} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = \frac{\Delta t}{b}b\frac{u_{i,j+1}-u_{i,j}}{\Delta t} λ = \frac{\Delta t}{b(\Delta x)^2} λ = \frac{0.002}{2(0.2)^2} =0.025

como λ<0.5 el método converge.

\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] = u[i,j+1]-u[i,j]

por el método explícito:

u[i,j+1] =\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] + u[i,j] u[i,j+1] =\lambda u[i+1,j]+(1-2\lambda)u[i,j]+\lambda u[i-1,j]

iteración i=1, j=0

u[1,1] =\lambda u[0,0]+(1-2\lambda)u[1,0]+\lambda u[2,0] u[1,1] =0.025(1) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.2\Big)+0.025 \cos \Big( \frac{3π}{2}0.4\Big)

iteración i=2, j=0

u[2,1] =\lambda u[1,0]+(1-2\lambda)u[2,0]+\lambda u[3,0] u[2,1] =0.025\cos \Big( \frac{3π}{2}0.2\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.4\Big) +0.025 \cos \Big( \frac{3π}{2}0.6\Big)

iteración i=3, j=0

u[3,1] =\lambda u[2,0]+(1-2\lambda)u[3,0]+\lambda u[4,0] u[3,1] =0.025\cos \Big( \frac{3π}{2}0.4\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.6\Big) +0.025 \cos \Big( \frac{3π}{2}0.8\Big)

iteración i=4, j=0

u[4,1] =\lambda u[3,0]+(1-2\lambda)u[4,0]+\lambda u[5,0] u[4,1] =0.025\cos \Big( \frac{3π}{2}0.6\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.8\Big) +0.025 (0)

continuar con las iteraciones en el algoritmo

Resultados con el algoritmo

Tabla de resultados
[[ 1.    1.    1.    1.    1.    1.    1.    1.    1.    1.  ]
 [ 0.59  0.58  0.56  0.55  0.54  0.53  0.53  0.52  0.51  0.5 ]
 [-0.31 -0.3  -0.3  -0.29 -0.28 -0.28 -0.27 -0.27 -0.26 -0.26]
 [-0.95 -0.93 -0.91 -0.89 -0.88 -0.86 -0.84 -0.82 -0.81 -0.79]
 [-0.81 -0.79 -0.78 -0.76 -0.74 -0.73 -0.71 -0.7  -0.68 -0.67]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]]

2Eva2022PAOII_T3 EDP Parabolica 02

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito,usando diferencias divididas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 1
Tb = 0
#T0 = 25
fx = lambda x: np.cos(3*np.pi/2*x)
# longitud en x
a = 0
b = 1
# Constante K
K = 2
# Tamaño de paso
dx = 0.2
dt = dx/100
# iteraciones en tiempo
n = 10

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
fi = fx(xi)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot): # igual con lazo for
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
    
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

s2Eva_2022PAOII_T2 EDO – población de protestantes en una sociedad

Ejercicio: 2Eva_2022PAOII_T2 EDO – población de protestantes en una sociedad

\frac{\delta}{\delta t}x(t) = b x(t) - d (x(t))^2 \frac{\delta}{\delta t}y(t) = b y(t) - d (y(t))^2 +r b (x(t)-y(t))

literal a

simplificando la nomenclatura

x' = b x - d x^2 y' = b y - d y^2 +r b (x-y)

sustituyendo constantes, y considerando x(0)=1 ; y(0)=0.01 ; h=0.5

x' = 0.02 x - 0.015 x^2 y' = 0.02 y - 0.015 y^2 +0.1(0.02) (x-y)

el planteamiento de Runge Kutta se hace junto a la primera iteración, además de encontrarse en las instrucciones con Python.

literal b

Se describen 3 iteraciones usando los resultados de la tabla con el algoritmo, para mostrar la comprensión del algoritmo.

t = 0
K1x = 0.5 (0.02 (1) – 0.015 (1)2 = 0.0025
K1y = 0.5(0.02 (0.01) – 0.015 (0.01)2 +0.1(0.02) (1-0.01)= 0.001089

K2x = 0.5 (0.02 (1+0.0025) – 0.015 (1+0.0025)2= 0.00248
K2y = 0.5(0.02 (0.01+0.00108) – 0.015 (0.01+0.00108)2+0.1(0.02) ((1+0.0025)-(0.01+0.00108)) = 0.001101

x1 = 1 + (1/2)(0.0025+0.00248) = 1.0025
y1 = 0.01 + (1/2)(0.001089+0.001101) = 0.01109
t1 = 0 + 0.5 =0.5

t=0.5
K1x = 0.5 (0.02 (1.0025) – 0.015 (1.0025)2 = 0.002487
K1y = 0.5(0.02 (0.01109) – 0.015 (0.01109)2 +0.1(0.02) (1.0025-0.01109)= 0.001101

K2x = 0.5 (0.02 (1.0025+ 0.002487) – 0.015 (1.0025+ 0.002487)2= 0.002474
K2y = 0.5(0.02 (0.01109+0.001101) – 0.015 (0.01109+0.001101)2+0.1(0.02) ((1.0025+ 0.002487)-(0.01109+0.001101)) = 0.001113

x2 = 1.0025 + (1/2)(0.002487+0.002474) = 1.0050
y2 = 0.01109 + (1/2)(0.001101+0.001113) = 0.01220
t2 = 0.5 + 0.5 = 1

t=1
K1x = 0.5 (0.02 (1.0050) – 0.015 (1.0050)2 = 0.002474
K1y = 0.5(0.02 (0.01220) – 0.015 (0.01220)2 +0.1(0.02) (1.0050-0.01220)= 0.001113

K2x = 0.5 (0.02 (1.0050+ 0.002474) – 0.015 (1.0050+ 0.002474)2= 0.002462
K2y = 0.5(0.02 (0.01220+0.001113) – 0.015 (0.01220+0.001113)2+0.1(0.02) ((1.0050+ 0.002474)-(0.01220+0.001113)) = 0.001126

x3 = 1.0050 + (1/2)(0.002474+0.002462) = 1.0074
y3 = 0.01220 + (1/2)(0.001113+0.001126) = 0.01332
t3 = 1 + 0.5 = 1.5

Resultado con el algoritmo

Para obtener los datos de las iteraciones, primero se ejecuta el algoritmo para pocas iteraciones.
Para la pregunta sobre 200 años, se incrementa las iteraciones a 2 por año y las condiciones iniciales, es decir 401 muestras.

 [ ti, xi, yi]
 [ ti, xi, yi]
[[0.0000e+00 1.0000e+00 1.0000e-02 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00]
 [5.0000e-01 1.0025e+00 1.1095e-02 2.5000e-03 1.0892e-03 2.4875e-03 1.1014e-03]
 [1.0000e+00 1.0050e+00 1.2203e-02 2.4875e-03 1.1014e-03 2.4749e-03 1.1136e-03]
 [1.5000e+00 1.0074e+00 1.3323e-02 2.4749e-03 1.1137e-03 2.4623e-03 1.1260e-03]
 [2.0000e+00 1.0099e+00 1.4455e-02 2.4624e-03 1.1260e-03 2.4497e-03 1.1384e-03]
 [2.5000e+00 1.0123e+00 1.5600e-02 2.4498e-03 1.1384e-03 2.4371e-03 1.1509e-03]
 [3.0000e+00 1.0148e+00 1.6757e-02 2.4371e-03 1.1509e-03 2.4245e-03 1.1634e-03]
 [3.5000e+00 1.0172e+00 1.7926e-02 2.4245e-03 1.1635e-03 2.4118e-03 1.1761e-03]
 [4.0000e+00 1.0196e+00 1.9109e-02 2.4118e-03 1.1761e-03 2.3991e-03 1.1888e-03]
...
 [1.9950e+02 1.3252e+00 1.1561e+00 ... 1.7202e-03 8.1217e-05 1.7059e-03]
 [2.0000e+02 1.3252e+00 1.1578e+00 ... 1.7060e-03 8.0418e-05 1.6918e-03]
 [2.0050e+02 1.3253e+00 1.1595e+00 ... 1.6919e-03 7.9628e-05 1.6778e-03]]
>>> 

Observación: La población identificada como protestante, continua creciendo, mientras que la proporción de «conformistas» se reduce según los parámetros indicados en el ejercicio. Los valores de natalidad y defunción cambian con el tiempo mucho más en años por otras variables, por lo que se deben realizar ajustes si se pretende extender el modelo.

2Eva2022PAOII_T2 poblacion protestantes
Instrucciones en Python

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

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

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

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
b = 0.02
d = 0.015
r = 0.1

# Ecuaciones
f = lambda t,x,y : (b-d*x)*x
g = lambda t,x,y : (b-d*y)*y + r*b*(x-y)

# Condiciones iniciales
t0 = 0
x0 = 1
y0 = 0.01

# parámetros del algoritmo
h = 0.5
muestras = 401

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, xi, yi, K1x, K1y, K2x, K2y]')
print(tabla)

# Grafica tiempos vs población
import matplotlib.pyplot as plt

plt.plot(ti,xi, label='xi poblacion')
plt.plot(ti,yi, label='yi protestante')

plt.title('población y protestantes')
plt.xlabel('t años')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

# gráfica xi vs yi
plt.plot(xi,yi)

plt.title('población y protestantes [xi,yi]')
plt.xlabel('x población')
plt.ylabel('y protestantes')
plt.grid()
plt.show()

s2Eva_2022PAOII_T1 Altura de cohete en 30 segundos

Ejercicio: 2Eva_2022PAOII_T1 Altura de cohete en 30 segundos

literal a

v = u \ln\Big(\frac{m_0}{m_0-qt}\Big) - gt v = 1800 \ln\Big(\frac{160000}{160000-2500t}\Big) - 9.8t

Seleccionando el método de Simpson de 3/8, se requieren al menos 3 tramos o segmentos para usarlo, que generan 4 muestras. El vector de tiempo se obtiene como:

v = lambda t: 1800*np.log(160000/(160000-2500*t))-9.8*t
a = 0
b = 30
tramos = 3
h = (b-a)/tramos
ti = np.linspace(a,b,tramos+1)
vi = v(ti)

siendo los vectores:

ti = [ 0. 10. 20. 30.]
vi = [ 0. 207.81826623 478.44820899 844.54060574]

la aplicación del método de Simpson de 3/8 es:

I = \frac{3}{8}(10) \Bigg(1800 \ln\Big(\frac{160000}{160000-2500(0)}\Big) - 9.8(0) +3(1800 \ln\Big(\frac{160000}{160000-2500(10)}\Big) - 9.8(10)) +3(1800 \ln\Big(\frac{160000}{160000-2500(20)}\Big) - 9.8(20)) +1800 \ln\Big(\frac{160000}{160000-2500(30)}\Big) - 9.8(30) \Bigg) = I = \frac{3}{8}(10) \Big(v(0)+3(v(10))+3(v(20))+v(30) \Big) I = \frac{3}{8}(10) \Big(0+3(207.81)+3(478.44)+844.54 \Big) I = 10887.52

literal b

para el primer segmento se usa t entre [0,10]

x_a = \frac{0+10}{2} + \frac{1}{\sqrt{3}}\frac{10-0}{2} = 7.88 x_b = \frac{0+10}{2} - \frac{1}{\sqrt{3}}\frac{10-0}{2} = 2.11 I = \frac{10-0}{2}\Big(v(7.88)+v(2.11)\Big)=995.79

para el 2do segmento se usa t entre [10,20]

x_a = \frac{10+20}{2} + \frac{1}{\sqrt{3}}\frac{20-10}{2} = 17.88 x_b = \frac{10+20}{2} - \frac{1}{\sqrt{3}}\frac{20-10}{2} = 12.11 I = \frac{20-10}{2}\Big(v(17.88)+v(12.11)\Big) =3368.42

para el 3er segmento se usa t entre [20,30]

x_a = \frac{20+30}{2} + \frac{1}{\sqrt{3}}\frac{30-20}{2} = 27.88 x_b = \frac{20+30}{2} - \frac{1}{\sqrt{3}}\frac{30-20}{2} = 22.11 I = \frac{30-20}{2}\Big(v(27.88)+v(22.11)\Big) = 6515.23 Altura = 995.79+ 3368.42 + 6515.23 = 10879.44

literal c

el error es la diferencia entre los métodos
error_entre = |10887.52-10879.44| = 8.079

Resultados con algoritmo

Método de Simpon 3/8
ti
[ 0. 10. 20. 30.]
vi
[ 0. 207.81826623 478.44820899 844.54060574]
Altura con Simpson 3/8 : 10887.52511781406
segmento Cuad_Gauss :    [995.792, 3368.421, 6515.231]
Altura Cuadratura Gauss: 10879.445437288954
diferencia s3/8 y Cuad_Gauss: 8.079680525106596
>>>

Instrucciones en Python

# 2Eva_2022PAOII_T1 Altura de cohete en 30 segundos
import numpy as np

# INGRESO
v = lambda t: 1800*np.log(160000/(160000-2500*t))-9.8*t
a = 0
b = 30
tramos = 3

# PROCEDIMIENTO literal a
def integrasimpson38_fi(xi,fi,tolera = 1e-10):
    ''' sobre muestras de fi para cada xi
        integral con método de Simpson 3/8
        respuesta es np.nan para tramos desiguales,
        no hay suficientes puntos.
    '''
    n = len(xi)
    i = 0
    suma = 0
    while not(i>=(n-3)):
        h  = xi[i+1]-xi[i]
        h1 = (xi[i+2]-xi[i+1])
        h2 = (xi[i+3]-xi[i+2])
        dh = abs(h-h1)+abs(h-h2)
        if dh<tolera:# tramos iguales
            unS38 = fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]
            unS38 = (3/8)*h*unS38
            suma = suma + unS38
        else:  # tramos desiguales
            suma = 'tramos desiguales'
        i = i + 3
    if (i+1)<n: # incompleto, tramos por calcular
        suma = 'tramos incompletos, faltan '
        suma = suma +str(n-(i+1))+' tramos'
    return(suma)

h = (b-a)/tramos
ti = np.linspace(a,b,tramos+1)
vi = v(ti)
altura = integrasimpson38_fi(ti,vi)

# SALIDA
print('Método de Simpon 3/8')
print('ti')
print(ti)
print('vi')
print(vi)
print('Altura con Simpson 3/8 :',altura)

# PROCEDIMIENTO literal b
# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    return(area)

area = 0
area_i =[]
for i in range(0,tramos,1):
    deltaA = integraCuadGauss2p(v,ti[i],ti[i+1])
    area = area + deltaA
    area_i.append(deltaA)
# SALIDA
print('segmento Cuad_Gauss :   ', area_i)
print('Altura Cuadratura Gauss:', area)

print('diferencia s3/8 y Cuad_Gauss:',altura-area)

import matplotlib.pyplot as plt
plt.plot(ti,vi)
plt.plot(ti,vi,'o')
plt.title('v(t)')
plt.xlabel('t (s)')
plt.ylabel('v (m/s)')
plt.grid()
plt.show()

s1Eva_2022PAOII_T3 Trayectoria de dron con polinomios

Ejercicio: 1Eva_2022PAOII_T3 Trayectoria de dron con polinomios

La variable independiente para la trayectoria es tiempo, con datos en el vector de ti.

ti  = [0, 1, 2, 3, 4]
xti = [2, 1, 3, 4, 2]
yti = [0, 1, 5, 1, 0]

Considerando que los puntos marcan posiciones por donde debe pasar el dron y se define la trayectoria, se usarán todos los puntos. Cada polinomio será de grado 4 al incluir los 5 puntos disponibles para cada eje.

Nota: podría usar polinomios de menor grado, siempre que considere que se debe completar la trayectoria y regresar al punto de salida.

px(t) = 2\frac{(t-1)(t-2)(t-3)(t-4)}{(0-1)(0-2)(0-3)(0-4)} + 1 \frac{(t-0)(t-2)(t-3)(t-4)}{(1-0)(1-2)(1-3)(1-4)} + 3 \frac{(t-0)(t-1)(t-3)(t-4)}{(2-0)(2-1)(2-3)(2-4)} + 4 \frac{(t-0)(t-1)(t-2)(t-4)}{(3-0)(3-1)(3-2)(3-4)} + 2 \frac{(t-0)(t-1)(t-2)(t-3)}{(4-0)(4-1)(4-2)(4-3)}

simplificando con el algoritmo:

px(t) = \frac{1}{12}t^4 - \frac{7}{6}t^3 + \frac{53}{12}t^2 - \frac{13}{3}t + 2

Realizando lo mismo con el algoritmo para polinomio de Lagrange se obtiene:

py(t) = \frac{11}{12}t^4 - \frac{22}{3}t^3 + \frac{205}{12}t^2 - \frac{29}{3}t

se muestra la gráfica de trayectorias por cada eje vs tiempo

Observaciones: La trayectoria usada tiene el mismo punto de salida como de retorno. La trayectoria presenta lóbulos que podrían ser reducidos y minimizar uso de recursos como bateria. Considere usar trazadores cúbicos y observe la misma gráfica de trayectorias x(t) vs y(t).

Resultado con el algoritmo:

Polinomio de Lagrange x: 
x**4/12 - 7*x**3/6 + 53*x**2/12 - 13*x/3 + 2
Polinomio de Lagrange y: 
11*x**4/12 - 22*x**3/3 + 205*x**2/12 - 29*x/3

Algoritmo en Python

# 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
ti  = [0,1,2,3,4]
xti = [2,1,3,4,2]
yti = [0,1,5,1,0]

# PROCEDIMIENTO
x = sym.Symbol('x')

def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if (j!=i):
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

# para ejex
polinomiox = interpola_lagrange(ti,xti)
polisimplex = polinomiox.expand()
px = sym.lambdify(x,polisimplex)

# para ejey
polinomioy = interpola_lagrange(ti,yti)
polisimpley = polinomioy.expand()
py = sym.lambdify(x,polisimpley)

# Puntos para la gráfica
muestras = 101
a = np.min(ti)
b = np.max(ti)
ti = np.linspace(a,b,muestras)
pxi = px(ti)
pyi = py(ti)

# SALIDA
print('Polinomio de Lagrange x: ')
print(polisimplex)
print('Polinomio de Lagrange y: ')
print(polisimpley)

# Gráfica
figura, enplano = plt.subplots()
plt.scatter(xti,yti, color='red')
plt.plot(pxi,pyi)
plt.ylabel('y(t)')
plt.xlabel('x(t)')
plt.title('trayectoria 2D')
plt.grid()

figura, entiempo = plt.subplots()
plt.plot(ti,pxi, label = 'px')
plt.plot(ti,pyi, label = 'py')
plt.legend()
plt.title('posicion en tiempo')
plt.xlabel('t')
plt.ylabel('p(t)')
plt.grid()

plt.show()

s1Eva_2022PAOII_T2 Admisión universitaria – cupos por recursos

Ejercicios: 1Eva_2022PAOII_T2 Admisión universitaria – cupos por recursos

Se requiere determinar la distribución de cupos en base a los costos relativos al promedio por estudiante para docencia, infraestructura y servicios mostrados en la tabla.

Costo referencial /carrera Mecatrónica Computación Civil Matemáticas
Docencia 1.5 0.9 0.6 0.7
Infraestructura 0.8 1.4 0.4 0.5
Servicios 0.45 0.55 1.1 0.5

Con los datos del total de recursos relativos al promedio por estudiante disponibles son docencia 271, infraestructura 250 y servicios 230.

1.5 a + 0.9 b + 0.6 c + 0.7 d = 271 0.8 a + 1.4 b + 0.4 c + 0.5 d = 250 0.45 a + 0.55 b + 1.1 c + 0.5 d = 230

se indica que en carreras como matemáticas de baja demanda, se establece el cupo de 10,

1.5 a + 0.9 b + 0.6 c + 0.7 (10) = 271 0.8 a + 1.4 b + 0.4 c + 0.5(10) = 250 0.45 a + 0.55 b + 1.1 c + 0.5(10) = 230

el sistema se convierte en:

1.5 a + 0.9 b + 0.6 c = 271 - 0.7 (10) 0.8 a + 1.4 b + 0.4 c = 250 - 0.5(10) 0.45 a + 0.55 b + 1.1 c = 230 - 0.5(10)

Para usar un método iterativo se convierte a matriz aumentada:

\begin{pmatrix} 1.5 & 0.9 & 0.6 & \Big| & 264 \\ 0.8 & 1.4 & 0.4 & \Big| & 245 \\ 0.45 & 0.55 & 1.1 &\Big| & 225 \end{pmatrix}

con pivoteo parcial por filas, la matriz aumentada se mantiene igual, pues los valores de la diagonal ya son los mayores posibles según el algoritmo.

Para un método iterativo se despeja una ecuación por cada incógnita.

a = \frac{1}{1.5}(264 - 0.9 b - 0.6 c) b = \frac{1}{1.4}(245 - 0.8 a - 0.4 c) c = \frac{}{1.1}(225 -0.45 a - 0.55 b)

Para los valores iniciales se consideran números mayores que cero, pues existen recursos para los cupos. No se admiten cupos negativos.

X_0 = [50,50,50]

Las iteraciones para el método iterativo de Gauss-Seidel

itera = 0

a = \frac{1}{1.5}(264 - 0.9 (50) - 0.6 (50)) = 126 b = \frac{1}{1.4}(245 - 0.8 (126) - 0.4 (50)) = 88.714 c = \frac{}{1.1}(225 -0.45 (126) - 0.55 (88.714)) =108.642 diferencia = [126-50, 88.714-50, 108.42-50] diferencia = [76, 38.714, 58.642] errado = max|[76, 38.714, 58.642]| =76 X = [126, 88.714, 108.42]

itera = 1

a = \frac{1}{1.5}(264 - 0.9 (88.714) - 0.6 (108.42)) = 79.314 b = \frac{1}{1.4}(245 - 0.8 (79.314) - 0.4 (108.42)) = 98.637 c = \frac{}{1.1}(225 -0.45 (79.314) - 0.55 (98.637)) =122.780 diferencia = [79.314-126, 88, 98.637-88.714, 122.780-108.42] diferencia = [46.685,9.922, 14.137] errado = max| [46.685,9.922, 14.137] | = 46.685

el error disminuye en la iteración

X = [79.314 , 98.637, 122.780]

itera = 2

a = \frac{1}{1.5}(264 - 0.9 (79.314) - 0.6 (122.780)) = 67.705 b = \frac{1}{1.4}(245 - 0.8 (67.705) - 0.4 (122.780)) = 101.230 c = \frac{}{1.1}(225 -0.45 (67.705) - 0.55 (101.230)) =126.232 diferencia = [67.705-79.314, 101.230-98.637, 126.232-122.780] diferencia = [-11.608, 2.594, 3.451] errado = max| [-11.608, 2.594, 3.451] | = 11.608

el error disminuye en la iteración, se considera que el método converge

X = [67.705 , 101.230, 126.232]

con el algoritmo se tiene como resultado:

[126.          88.71428571 108.64285714]
[76.         38.71428571 58.64285714]

[ 79.31428571  98.63673469 122.78033395]
[46.68571429  9.92244898 14.13747681]

[ 67.7058256  101.23086138 126.23218611]
[11.60846011  2.59412669  3.45185216]

[ 64.76860873 101.92302755 127.08769174]
[2.93721688 0.69216617 0.85550564]

[ 64.01110677 102.11145563 127.30336487]
[0.75750196 0.18842808 0.21567312]

[ 63.81178067 102.16373537 127.3587675 ]
[0.1993261  0.05227973 0.05540263]

[ 63.75825178 102.17849398 127.37328637]
[0.05352889 0.01475862 0.01451887]

[ 63.74358906 102.18272443 127.37716953]
[0.01466272 0.00423045 0.00388316]

[ 63.73949753 102.18395297 127.37822907]
[0.00409153 0.00122854 0.00105954]

[ 63.73833659 102.18431364 127.37852366]
[0.00116094 0.00036067 0.0002946 ]

[ 63.73800235 102.18442047 127.37860699]
[3.34240232e-04 1.06824300e-04 8.33224905e-05]

respuesta X: 
[[ 63.73800235]
 [102.18442047]
 [127.37860699]]
verificar A.X=B: 
[[264.00014614]
 [245.00003333]
 [225.        ]]
>>>

se interpreta la respuesta como la parte entera de la solución:

cupos = [ 63, 102 , 127]

s1Eva_2022PAOII_T1 Esfera flotando en agua

Ejercicio: 1Eva_2022PAOII_T1 Esfera flotando en agua

Según el principio de Arquímedes, la fuerza de flotación o empuje es igual al peso de el fluido desplazado por la porción sumergida de un objeto.

F_{empuje} = F_{peso} \rho_{agua} V_{sumergido} \text{ } g = \rho_{esfera}V_{esfera} \text{ } g V_{sumergido} = \frac{\rho_{esfera}}{\rho_{agua}}V_{esfera} V_{esfera} - V_{sobreagua} = \frac{\rho_{esfera}}{\rho_{agua}}V_{esfera} V_{sobreagua} = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) V_{esfera} V_{esfera} = \frac{4}{3}\pi r^3 \frac{\pi h^2}{3}(3r-h) = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) \frac{4}{3}\pi r^3 h^2(3r-h) = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) 4 r^3

El planteamiento para la búsqueda de raíces es f(x) = 0, que para este caso será:

f(h) = h^2(3r-h) - \Big( 1 - \frac{\rho_{esfera}}{\rho_{agua}}\Big) 4 r^3 = 0

usando los valores dados para el ejercicio, r=1 y ρesfera = 200 Kg/m3 y ρagua    = 1000 kg/m3 se tiene que:

f(h) = h^2(3-h) - \Big( 1 - \frac{200}{1000}\Big) 4 f(h) = h^2(3-h) - \frac{16}{5}

Se observa la gráfica de f(h) en el intervalo de h entre[0,2] interpretado como totalmente sumergida y totalmente flotando sobre el agua, confirmando que existe una raíz

Para el caso de aplicar el método del punto fijo se plantea que x=g(x),

h = g(h) h^2(3-h) = \frac{16}{5}

con lo que se puede plantear dos ecuaciones al despejar h

h = \sqrt{ \frac{16}{5(3-h)}} h = 3-\frac{16}{5 h^2}

Iteraciones de la primera ecuación

itera = 0 ; h = h0 = 0.5 ;

g(h) = \sqrt{ \frac{16}{5(3-0.5)}} = 1.1313 tramo = |1.1313-0.5|=0.6313

itera = 1 ; h = 1.1313 ;

g(h) = \sqrt{ \frac{16}{5(3-1.1313)}} = 1.3086 tramo = |1.3086-1.1313| = 0.1772

itera = 2 ; h = 1.3086 ;

g(h) = \sqrt{ \frac{16}{5(3-1.3086)}} = 1.3754 tramo = |1.3754-1.3086| = 0.0668

Observando los errores o tramos en cada iteración se tiene que se reduce, el método converge.


resultados.txt

x,g(x),tramo
0.5 1.131370849898476 0.631370849898476
1.131370849898476 1.308619626317284 0.17724877641880799
1.308619626317284 1.3754802083033437 0.06686058198605971
1.3754802083033437 1.4035002223557855 0.02802001405244181
1.4035002223557855 1.4157629993958152 0.012262777040029649
1.4157629993958152 1.4212317895316 0.005468790135784829
1.4212317895316 1.4236912066694054 0.0024594171378053975
1.4236912066694054 1.424801422465215 0.0011102157958096104
1.424801422465215 1.4253034412081806 0.0005020187429656264
raiz: 1.4253034412081806

Algoritmo en Python

# Algoritmo de punto fijo
# [a,b] intervalo de búsqueda
# error = tolera

import numpy as np
import matplotlib.pyplot as plt

def puntofijo(gx,a,tolera, iteramax = 15):
    i = 1 # iteración
    b = gx(a)
    tramo = abs(b-a)
    print('x,g(x),tramo')
    print(a,b,tramo)
    while(tramo>=tolera and i<=iteramax ):
        a = b
        b = gx(a)
        tramo = abs(b-a)
        print(a,b,tramo)
        i = i + 1
    respuesta = b
    
    # Validar respuesta
    if (i>=iteramax ):
        respuesta = np.nan
    return(respuesta)

# PROGRAMA ---------
# INGRESO
fx = lambda h: h**2*(3-h)-16/5
gx = lambda h: np.sqrt(16/(5*(3-h)))

#fx = lambda h: h**2*(3-h)-16/5
#gx = lambda h: 3-16/(5*(h**2))

x0 = 0.5
tolera = 0.001
iteramax = 50  # itera máximo
a = 0     # intervalo
b = 2
muestras = 51  # gráfico

# PROCEDIMIENTO
respuesta = puntofijo(gx,x0,tolera)

# SALIDA
print('raiz:',respuesta)

hi = np.linspace(a,b,muestras)
fi = fx(hi)
gi = gx(hi)
plt.plot(hi,fi,label='f(h)')
plt.plot(hi,gi,label='g(h)')
plt.plot(hi,hi,label='Identidad')
plt.axhline(0,color='grey')
plt.grid()
plt.xlabel('h')
plt.ylabel('f(h)')
plt.title('esfera sumergida')
plt.legend()
plt.show()

s2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro

Ejercicio: 2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro

Para la ecuación dada con Δx = 1/3, Δt = 0.02, en una revisíón rápida para cumplir la convergencia dt<dx/10, condición que debe verificarse con la expresión obtenida para λ al desarrollar el ejercicio.

\frac{\partial U}{\partial t} - \frac{1}{9} \frac{\partial ^2 U}{\partial x^2} = 0 0 \leq x \leq 2, t>0

literal a. grafica de malla

literal b. Ecuaciones de diferencias divididas a usar

\frac{\partial U}{\partial t} - \frac{1}{9} \frac{\partial ^2 U}{\partial x^2} = 0 \frac{\partial ^2 U}{\partial x^2} = 9 \frac{\partial U}{\partial t} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = 9 \frac{u_{i,j+1}-u_{i,j}}{\Delta t}

se agrupan las constantes,

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

literal d Determine el valor de λ

\lambda = \frac{\Delta t}{9(\Delta x)^2} =\frac{0.02}{9(1/3)^2} = 0.02

valor de λ que es menor que 1/2, por lo que el método converge.

continuando luego con la ecuación general,

\lambda \Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) = u[i,j+1]-u[i,j] \lambda u[i-1,j]-2 \lambda u[i,j] + \lambda u[i+1,j] \Big) = u[i,j+1]-u[i,j]

literal c. Encuentre las ecuaciones considerando las condiciones dadas en el problema.

\lambda u[i-1,j]+(1-2 \lambda ) u[i,j] + \lambda u[i+1,j] = u[i,j+1]

el punto que no se conoce su valor es u[i,j+1] que es la ecuación buscada.

u[i,j+1] = \lambda u[i-1,j]+(1-2 \lambda ) u[i,j] + \lambda u[i+1,j]

literal e iteraciones

iteración  i=1, j=0

u[1,1] = \lambda u[0,0]+(1-2 \lambda ) u[1,0] + \lambda u[2,0] u[1,1] =0.02 \cos \Big( \frac{\pi}{2}(0-3)\Big) + (1-2(0.02) ) \cos \Big( \frac{\pi}{2}\big(\frac{1}{3}-3\big)\Big) + 0.02 \cos \Big( \frac{\pi}{2}\big( \frac{2}{3}-3\big) \Big) u[1,1] =0.02(0)+(0.96)(-0.5)+0.02(-0.8660)=-0.4973

iteración  i=2, j=0

u[2,1] = \lambda u[1,0]+(1-2 \lambda ) u[2,0] + \lambda u[3,0] u[2,1] = 0.02 \cos \Big( \frac{\pi}{2}(\frac{1}{3}-3)\Big) + (1-2(0.02) ) \cos \Big( \frac{\pi}{2}(\frac{2}{3}-3)\Big)+ + 0.02 \cos \Big( \frac{\pi}{2}\big(\frac{3}{3}-3\big)\Big) u[2,1] = 0.02 (-0.5) + (0.96 ) (-0.866025) + 0.02 (-1) =-0.8614

iteración  i=3, j=0

u[3,1] = \lambda u[2,0]+(1-2 \lambda ) u[3,0] + \lambda u[4,0] u[3,1] = 0.02 \cos \Big( \frac{\pi}{2}\big( \frac{2}{3}-3\big)\Big)+(1-2 (0.02) ) \cos \Big( \frac{\pi}{2}(1-3)\Big) + + 0.02 \cos \Big( \frac{\pi}{2}\big(\frac{4}{3}-3\big)\Big) u[3,1] = 0.02 (-0.866025)+(0.96 ) (-1) + 0.02 (-0,866025) = -0,9946

literal f

la cotas de errores de truncamiento en la ecuación corresponden a segunda derivada O(hx2) y el de primera derivada O(ht), al reemplazar los valores será la suma}

O(hx2) + O(ht) = (1/3)2 + 0.02 = 0,1311

literal g

Resultados usando el algoritmo en Python

Tabla de resultados
[[ 0.      0.      0.      0.      0.      0.      0.      0.      0.       0.    ]
 [-0.5    -0.4973 -0.4947 -0.492  -0.4894 -0.4867 -0.4841 -0.4815 -0.479   -0.4764]
 [-0.866  -0.8614 -0.8568 -0.8522 -0.8476 -0.8431 -0.8385 -0.8341 -0.8296  -0.8251]
 [-1.     -0.9946 -0.9893 -0.984  -0.9787 -0.9735 -0.9683 -0.9631 -0.9579  -0.9528]
 [-0.866  -0.8614 -0.8568 -0.8522 -0.8476 -0.8431 -0.8385 -0.8341 -0.8296  -0.8251]
 [-0.5    -0.4973 -0.4947 -0.492  -0.4894 -0.4867 -0.4841 -0.4815 -0.479   -0.4764]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.       0.    ]]

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito, usando diferencias finitas
# 2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 0
Tb = 0
T0 = lambda x: np.cos((np.pi/2)*(x-3))
# longitud en x
a = 0.0
b = 2.0
# Constante K
K = 9
# Tamaño de paso
dx = 1/3
dt = 0.02
tramos = int(np.round((b-a)/dx,0))
muestras = tramos + 1
# iteraciones en tiempo
n = 10

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.linspace(a,b,muestras)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,j]= Ta
u[1:ultimox,j] = T0(xi[1:ultimox])
u[ultimox,j] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot):
    u[0,j+1] = Ta
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    u[m-1,j+1] = Tb
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

s2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

Ejercicio: 2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

La corriente del inductor y(t) para t≥0 se deriva para tener la expresión solo derivadas:

\frac{\delta}{\delta t}y(t) + 2 y(t) + 5 \int_{-\infty}^t y(\tau) \delta \tau = 10 \mu(t)

Para t>0 que es donde transcurre el experimento, el escalón es una constante, se tiene que:

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

tomando las condiciones iniciales dadas para t=0, y(0)=2, y'(0)=–4

literal a

EL resultadoes perado es el planteamiento del problema. Se reescribe la ecuación con la nomenclatura simplificada y se resordena segun el modelo del método:

y'' = - 2y' - 5 y

luego se sustituye la variable y se convierte a las ecuaciones:

z =y' = f_x(t,y,z) z' = - 2z - 5 y = g_z(t,y,z)

se usa una tabla para llevar el registro de operaciones:

Se plantea las operaciones:

K1y = h * f(ti,yi,zi)
K1z = h * g(ti,yi,zi)

K2y = h * f(ti+h, yi + K1y, zi + K1z)
K2z = h * g(ti+h, yi + K1y, zi + K1z)

yi = yi + (K1y+K2y)/2
zi = zi + (K1z+K2z)/2
ti = ti + h

literal b

El resultado esperado es la aplicación correcta de los valores en las expresiones para al menos tres iteraciones usando h=0.01

itera = 0

K1y = 0.01 y'(0) = 0.01(-4) = -0.04 K1z = 0.01 (- 2z(0) - 5 y(0)) = 0.01(- 2(-4) - 5 (2)) = -0.02 K2y = 0.01 (-4-0.02) = -0.0402 K2z = 0.01 (-2(-4-0.02)-5(2-0.04)) = -0.0176 yi = yi + \frac{K1y+K2y}{2} = 2+\frac{-0.04-0.0402} {2} = 1.9599 zi = zi + \frac{K1z+K2z}{2} = -4 +\frac{-0.02-0.0176}{2} = -4.0188 ti = ti + h = 0+0.01 = 0.01

itera = 1

K1y = 0.01(-4.0188) = -0.040188 K1z = 0.01(- 2(-4.0188) - 5 (1.9599)) = -0.0176 K2y = 0.01 (-4.0188-0.0176) = -0.0403 K2z = 0.01 (-2(-4.0188-0.0176)-5(1.9599-0.040188)) = -0.0152 yi = 1.9599 +\frac{-0.040188-0.0403} {2} = 1.9196 zi = -4.0188 +\frac{-0.0176-0.0152}{2} = -4.0352 ti = ti + h = 0.01+0.01 = 0.02

itera = 2

K1y = 0.01(-4.0352) = -0.040352 K1z = 0.01(- 2(-4.0352) - 5 (1.9196)) = -0.0152 K2y = 0.01 (-4.0352-0.0152) = -0.0405 K2z = 0.01 (-2(-4.0352-0.0152)-5(1.9196-0.040352)) = -0.0129 yi = 1.9196 +\frac{-0.040352-0.0405} {2} =1.8792 zi = -4.0352 +\frac{-0.0152-0.0129}{2} = -4.0494 ti = ti + h = 0.02+0.01 = 0.03

Resultados con el algoritmo en Python

   ti,   yi,    zi,      K1y,    K1z,    K2y,     K2z
[[ 0.00  2.0000 -4.0000  0.0000  0.0000  0.0000   0.0000]
 [ 0.01  1.9599 -4.0188 -0.0400 -0.0200 -0.0402  -0.0176]
 [ 0.02  1.9196 -4.0352 -0.0401 -0.0176 -0.0403  -0.0152]
 [ 0.03  1.8792 -4.0494 -0.0403 -0.0152 -0.0405  -0.0129]
...

Literal c

Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

por lo que el error está en el orden de (0.01)3 = 0.000001


Literal d

Se requiere presentar el resultado para el intervalo t entre [0,5]. Siendo el tamaño de paso h=0.01 que es pequeño, se requieren realizar (5-0)/0.01=500 iteraciones, que es más práctico realizarlas usando el algoritmo.

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,7),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0,0,0,0,0]
    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,K1y,K1z,K2y,K2z]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -2*z -5*y + 0

t0 = 0
y0 = 2
z0 = -4

h = 0.01
tn = 5
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=4)
print('ti, yi, zi, K1y, K1z, K2y, K2z')
print(tabla)

# GRAFICA
plt.plot(ti,yi, color = 'orange', label='y_RK(t)')
plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('y(t) con Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Nota: En el curso TELG1001 Señales y Sistemas, la solución se realiza con Transformadas de Laplace

s2Eva_2022PAOI_T1 Comparar integrales numéricos Simpson y Cuadratura de Gauss

Ejercicio: 2Eva_2022PAOI_T1 Comparar integrales numéricos Simpson y Cuadratura de Gauss

Literal a. Integral con Simpson 1/3

Para la ecuación en el intervalo x entre [0,3] aplicando dos veces la fórmula en el intervalo se requieren al menos dos tramos cada Simpson de 1/3. Por lo que la cantidad de tramos es 4 (dos por cada formula, y dos fórmulas) que corresponden a 5 muestras.

El tamaño de paso se calcula como:

h = \frac{b-a}{tramos}=\frac{3-0}{4} = \frac{3}{4} = 0.75

representados en una gráfica como

A = \int_0^3 \frac{e^x \sin(x)}{1+x^2} \delta x

con lo que se define la función del integral f(x)

f(x) = \frac{e^x \sin(x)}{1+x^2}

Con lo que aplicando la fórmula se puede obtener los valores de las muestras:

xi= [0. 0.75   1.5    2.25   3.    ]
fi= [0. 0.9235 1.3755 1.2176 0.2834]

Nota: realizar las expresiones completas para las fórmulas si no adjunta el algoritmo en Python

Aplicando Simpson de 1/3 en cada tramo se tiene:

A_s= \frac{1}{3} \Big( \frac{3}{4} \Big ) \Big( 0 + 4(0.9235) + 1.3755 \Big) + + \frac{1}{3} \Big( \frac{3}{4} \Big ) \Big( 1.3755 + 4(1.2176) +0.2834 \Big) A_s = 2.8998

Literal b. Integral con Cuadratura de Gauss

Para usar dos veces el método de Cuadratura de Gauss se usan dos intervalos, con lo que las muestras en x serán:

xj= [0. 1.5 3. ]

se calculan los valores para el tramo [0, 1.5]:

x_{a1} = \frac{0+1.5}{2} - \frac{1}{\sqrt{3}}\frac{1.5-0}{2} = 0.3169 x_{b1} = \frac{0+1.5}{2} + \frac{1}{\sqrt{3}}\frac{1.5-0}{2} = 1.1830 A_{g1} =\frac{1.5-0}{2} \Big( f(0.3169)+f(1.1830)\Big) = 1.2361

se calculan los valores para el tramo [1.5, 3]:

x_{a2} = \frac{1.5+3}{2} - \frac{1}{\sqrt{3}}\frac{3-1.5}{2} = 1.8169 x_{b2} = \frac{1.5+3}{2} + \frac{1}{\sqrt{3}}\frac{3-1.5}{2} = 2.6830 A_{g2} =\frac{3-1.5}{2} \Big( f(1.8169)+f(2.6830)\Big) = 1.6329

El total del integral para el intervalo  [0,3]

A_g = A_{g1} + A_{g2} = 2.8691

Al comparar los resultados entre los métodos del literal a y b

errado = |A_s - A_g| = 2.8998 - 2.8691 = 0.0307

Instrucciones integradas en Python

# 2Eva_2022PAOI_T1
# Comparar integrales numéricos Simpson
#   y Cuadratura de Gauss
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: (np.exp(x)*np.sin(x))/(1+x**2)
a = 0
b = 3

# PROCEDIMIENTO
# Aplicando Simpson 1/3
tramos = 4
muestras = tramos+1
xi = np.linspace(a,b,muestras)
fi = fx(xi)
hs = xi[1]-xi[0]
As = (1/3)*(3/4)*(fi[0]+4*fi[1]+fi[2])
As = As + (1/3)*(3/4)*(fi[2]+4*fi[3]+fi[4])
erradoS = 2*(hs**5)

# Aplicando Cuadratura de Gauss
tramosG = 2
muestrasG = tramosG+1
xj = np.linspace(a,b,muestrasG)
hg = xj[1]-xj[0]

xa1 = (xj[0]+xj[1])/2 - (1/np.sqrt(3))*(xj[1]-xj[0])/2
xb1 = (xj[0]+xj[1])/2 + (1/np.sqrt(3))*(xj[1]-xj[0])/2
Ag1 = (hg/2)*(fx(xa1)+fx(xb1))

xa2 = (xj[1]+xj[2])/2 - (1/np.sqrt(3))*(xj[2]-xj[1])/2
xb2 = (xj[1]+xj[2])/2 + (1/np.sqrt(3))*(xj[2]-xj[1])/2
Ag2 = (hg/2)*(fx(xa2)+fx(xb2))

Ag = Ag1 + Ag2

# error entre métodos
errado = np.abs(As-Ag)

# SALIDA
print('xi=',xi)
print('fi=',fi)
print('A Simpson =', As)
print('error Truncamiento Simpson 2*(h^5):',
      erradoS)
print('Cuadratura de Gauss xa,xb por tramos:',
      [xa1,xb1,xa2,xb2])
print('  fx(xa),fx(xb) por tramos:',
      [fx(xa1),fx(xb1),fx(xa2),fx(xb2)])
print('  integral de cada tramo:', [Ag1,Ag2])
print('A Gauss =', Ag)
print('errado entre Simpson y Gauss',errado)

# Grafica con mejor resolucion
xk = np.linspace(a,b,5*tramos+1)
fk = fx(xk)

plt.plot(xk,fk)
plt.plot(xi,fi,'o')
for unx in xi:
    plt.axvline(unx,color='red',
                linestyle='dotted')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('(np.exp(x)*np.sin(x))/(1+x**2)')
plt.grid()
plt.show()