s2Eva_2023PAOI_T3 EDP elíptica, placa rectangular con frontera variable

Ejercicio: 2Eva_2023PAOI_T3 EDP elíptica, placa rectangular con frontera variable

Dada la EDP elíptica,

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = \Big( x^2 + y^2 \Big) e^{xy} 0 \lt x \lt 1 0 \lt y \lt 0.5

Se convierte a la versión discreta usando diferencias divididas centradas, según se puede mostrar con la gráfica de malla.

literal b

EDP eliptica rectangular frontera variable

literal a


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

literal c

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

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

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

\lambda= \frac{\Delta y^2}{\Delta x^2} = 1

se simplifica el coeficiente en λ =1

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

Se agrupan los términosiguales

u[i-1,j]-4u[i,j]+ u[i+1,j] + u[i,j-1] +u[i,j+1] = \Big( x^2 + y^2 \Big) e^{xy}

Se desarrollan las iteraciones para tres rombos y se genera el sistema de ecuacioens a resolver.
para i=1,j=1

u[0,1]-4u[1,1]+ u[2,1] + u[1,0] +u[1,2] = \Big( x[1]^2 + y[1]^2 \Big) e^{x[1]y[1]} 1-4u[1,1]+ u[2,1] + 1 + \frac{1}{8} = \Big( 0.25^2 + 0.25^2 \Big) e^{(0.25) (0.25)} -4u[1,1]+ u[2,1] = \Big( 0.25^2 + 0.25^2 \Big) e^{(0.25) (0.25)} - \frac{1}{8}

para i=2, j=1

u[1,1]-4u[2,1]+ u[3,1] + u[2,0] +u[2,2] = \Big( x[2]^2 + y[1]^2 \Big) e^{x[2]y[1]} u[1,1]-4u[2,1]+ u[3,1] + 1 + 0.25 = \Big( 0.5^2 + 0.25^2 \Big) e^{(0.5)(0.25)} u[1,1]-4u[2,1]+ u[3,1] = \Big( 0.5^2 + 0.25^2 \Big) e^{(0.5)(0.25)} -1.25

para i=3, j=1

u[2,1]-4u[3,1]+ u[4,1] + u[3,0] +u[3,2] = \Big( x[3]^2 + y[1]^2 \Big) e^{x[3]y[1]} u[2,1]-4u[3,1]+ 0.25 + 0 + \frac{3}{8} = \Big( 0.75^2 + 0.25^2 \Big) e^{(0.75)(0.25)} u[2,1]-4u[3,1] = \Big( 0.75^2 + 0.25^2 \Big) e^{(0.75)(0.25)} - 0.25 - \frac{3}{8}

con lo que se puede crear un sistema de ecuaciones y resolver el sistema para cada punto desconocido

\begin{pmatrix} -4 & 1 & 0 & \Big| & 0.008061807364732415 \\ 1 & -4 & 1 & \Big| & -0.8958911084166168 \\0 & 1 & -4 &\Big| & 0.1288939058881129 \end{pmatrix}

se obtiene los resultados para:

u[1,1] = 0.05953113
u[2,1] = 0.24618634
u[3,1] = 0.02932311

>>> import numpy as np
>>> (0.25**2+0.25**2)*np.exp(0.25*0.25) - 1/8
0.008061807364732415
>>> (0.5**2+0.25**2)*np.exp(0.5*0.25) - 1.25
-0.8958911084166168
>>> (0.75**2+0.25**2)*np.exp(0.75*0.25) - 0.25 -3/8
0.1288939058881129
>>> A=[[-4,1,0],[1,-4,1],[0.0,1.0,-4.0]]
>>> B = [0.008061807364732415, -0.8958911084166168, 0.1288939058881129]
>>> np.linalg.solve(A,B)
array([0.05953113, 0.24618634, 0.02932311])

s2Eva_2023PAOI_T2 Péndulo vertical amortiguado

Ejercicio: 2Eva_2023PAOI_T2 Péndulo vertical amortiguado

literal a

\frac{d^2 \theta}{dt^2} = -\mu \frac{d\theta}{ dt}-\frac{g}{L}\sin (\theta)

Se simplifica su forma a:

\frac{d\theta}{dt}= z = f_t(t,\theta,z) \frac{d^2\theta }{dt^2}= z' = -\mu z -\frac{g}{L}\sin (\theta) = g_t(t,\theta,z)

se usan los valores dados: g = 9.81 m/s2, L = 2 m

f_t(t,\theta,z) = z g_t(t,\theta,z) = - 0.5 z -\frac{9.81}{2}\sin (\theta)

y los valores iniciales para la tabla: θ(0) = π/4 rad, θ’ (0) = 0 rad/s, se complementan los valores en la medida que se aplica el desarrollo.

ti θ(ti) θ'(ti)=z
0 π/4 0
0.2 0.7161 -0.6583
0.4 0.5267 -0.1156
0.6 0.2579 -0.1410

literal b

Iteración 1:  ti = 0 ; yi = π/4 ; zi = 0

K1y = h * ft(ti,yi,zi) 
    = 0.2*(0) = 0
K1z = h * gt(ti,yi,zi) 
    = 0.2*(-0.5(0) -(9.81/2)sin (π/4) = -0.6930
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(0-0.6930)= -0.1386
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-0.5(0-0.6930) -(9.81/2)sin(π/4-0) 
    = -0.6237

yi = yi + (K1y+K2y)/2 
   = π/4+ (0+(-0.1386))/2 = 0.7161
zi = zi + (K1z+K2z)/2 
   = 0+(-0.6930-0.6237)/2 = -0.6583
ti = ti + h = 0 + 0.2 = 0.2

estimado[i] = [0.2,0.7161,-0.6583]

Iteración 2:  ti = 0.2 ; yi = 0.7161 ; zi = -0.6583

K1y = h * ft(ti,yi,zi) 
    = 0.2*(-0.6583) = -0.1317
K1z = h * gt(ti,yi,zi) 
    = 0.2*(- 0.5 ( -0.6583) -(9.81/2)sin (0.7161) 
    = -0.5775
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-0.6583 -0.5775)= -0.2472
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(- 0.5 (-0.6583 -0.5775) -(9.81/2)sin(0.7161-0.1317) 
    = -0.4171

yi = yi + (K1y+K2y)/2 
   = 0.7161 + (-0.1317-0.2472)/2 = 0.5267
zi = zi + (K1z+K2z)/2 
   = -0.6583+(-0.5775-0.4171)/2 = -0.1156
ti = ti + h = 0.2 + 0.2 = 0.4

estimado[i] = [0.4,0.5267,-0.1156]

Iteración 3:  ti = 0.4 ; yi = 0.5267 ; zi = -1.156

K1y = h * ft(ti,yi,zi) 
    = 0.2*(-1.156) = -0.2311
K1z = h * gt(ti,yi,zi) 
    = 0.2*(- 0.5(-1.156) -(9.81/2)sin (0.5267) 
    = -0.3771
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-1.156 -0.3771)= -0.3065
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(- 0.5 ( -1.156 -0.3771) -(9.81/2)sin(0.5267-0.2311) 
    = -0.1322

yi = yi + (K1y+K2y)/2 
   = 0.5267 + (-0.2311-0.3065)/2 = 0.2579
zi = zi + (K1z+K2z)/2 
   = -1.156+(-0.3771-0.1322)/2 = -1.410
ti = ti + h = 0.4 + 0.2 = 0.6

estimado[i] = [0.6,0.2579,-1.410]

literal c

resultados del algoritmo:

[ t, 		 y, 	 dyi/dti=z,  K1y,	 K1z,	    K2y,      K2z]
[[ 0.000e+00  7.854e-01  0.000e+00  0.000e+00  0.000e+00  0.000e+00   0.000e+00]
 [ 2.000e-01  7.161e-01 -6.583e-01  0.000e+00 -6.930e-01 -1.386e-01  -6.237e-01]
 [ 4.000e-01  5.267e-01 -1.156e+00 -1.317e-01 -5.775e-01 -2.472e-01  -4.171e-01]
 [ 6.000e-01  2.579e-01 -1.410e+00 -2.311e-01 -3.771e-01 -3.065e-01  -1.322e-01]
 [ 8.000e-01 -3.508e-02 -1.377e+00 -2.820e-01 -1.089e-01 -3.038e-01    1.756e-01]
...

péndulo amortiguado 03

con h=0.2 se tienen 1/0.2 = 5 tramos por segundo, por lo que para 10 segundo serán 50 tramos. La cantidad de muestras = tramos + 1(valor inicial) = 51

con lo que se puede usar el algoritmo en EDO Runge-Kutta d2y/dx2

literal d

Se observa que la respuesta es oscilante y amortiguada en magnitud como se esperaba según el planteamiento. Con el tiempo se estabilizará en cero.

Instrucciones en Python

# 2Eva_2023PAOI_T2 Péndulo vertical amortiguado
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,z0]
    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)

# INGRESO theta = y
g = 9.8
L  = 2
ft = lambda t,y,z: z
gt = lambda t,y,z: -0.5*z +(-g/L)*np.sin(y)

t0 = 0
y0 = np.pi/4
z0 = 0
h = 0.2
muestras = 51

# PROCEDIMIENTO
tabla = rungekutta2_fg(ft,gt,t0,y0,z0,h,muestras)

# SALIDA
np.set_printoptions(precision=3)
print(' [ t, \t\t y, \t dyi/dti=z, K1y,\t K1z,\t K2y,\t K2z]')
print(tabla)

# Grafica
ti = np.copy(tabla[:,0])
yi = np.copy(tabla[:,1])
zi = np.copy(tabla[:,2])
plt.subplot(121)
plt.plot(ti,yi)
plt.grid()
plt.xlabel('ti')
plt.title('yi')
plt.subplot(122)
plt.plot(ti,zi, color='green')
plt.xlabel('ti')
plt.title('dyi/dti')
plt.grid()
plt.show()

s2Eva_2023PAOI_T1 Material para medalla de academia

Ejercicio: 2Eva_2023PAOI_T1 Material para medalla de academia

medalla area con integral numerico f(x) = 2-8\Big( \frac{1}{2} - x \Big)^2

0 \le x \lt 1 g(x) = -\Big( 1-x\Big)\ln \Big( 1- x \Big)

Para f(x) se usará Simpson de 1/3 que requiere al menos dos tramos para aplicar:

a. Realice el planteamiento de las ecuaciones para el ejercicio.

I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)]

b. Describa el criterio usado para determinar el número de tramos usado en cada caso.

h = \frac{b-a}{2} = \frac{1-0}{2} = 0.5

c. Desarrolle las expresiones completas del ejercicio para cada función.

I_{fx}\cong \frac{0.5}{3}[f(0)+4f(0.5) + f(1)] f(0) = 2-8\Big( \frac{1}{2} - (0) \Big)^2 = 0 f(0.5) = 2-8\Big( \frac{1}{2} - (0.5) \Big)^2 = 2 f(1) = 2-8\Big( \frac{1}{2} - (1) \Big)^2 = 0 I_{fx} = \frac{1}{6}[0+4(2) + 0] = \frac{8}{6} = \frac{4}{3} = 1.3333

cota de error O(h5) = O(0.55)= O(0.03125)

Para g(x) se usará Simpson de 3/8 que requiere al menos tres tramos para aplicar:

I\cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)] h = \frac{b-a}{3} = \frac{1-0}{3} = 0.3333 I_{gx}\cong \frac{3(0.3333)}{8}[f(0)+3f(0.3333) +3 f(0.6666)+f(1)] g(0) = -\Big( 1-0\Big)\ln \Big( 1- 0 \Big) = 0 g(0.3333) = -\Big( 1-0.3333\Big)\ln \Big( 1- 0.3333 \Big) = 0.2703 g(0.6666) = -\Big( 1-0.6666\Big)\ln \Big( 1- 0.6666 \Big) = 0.3662 g(0.9999) = -\Big( 1-0.9999\Big)\ln \Big( 1- 0.9999 \Big) = 0

para la evaluación numérica de 1 se usa un valor muy cercano desplazado con la tolerancia aplicada.

I_{gx}\cong \frac{3(0.3333)}{8}[0+3(0.2703) + 3(0.3662)+0] = 0.2387

d. Indique el resultado obtenido para el área requerida y la cota de error
Area = I_{fx} – I_{gx} = 1.3333 – 0.2387 = 1.0945

cota de error = O(0.03125) + O(0.00411) = 0.03536

e. Encuentre el valor del tamaño de paso si se requiere una cota de error de 0.00032

Si el factor de mayor error es de Simpson 1/3, se considera como primera aproximación que:

cota de error O(h5) = O(0.00032), h = (0.00032)(1/5) = 0.2
es decir el número de tramos es de al menos (b-a)/tramos = 0.2 , tramos = 5.
El número de tramos debe ser par en Simpson de 1/3, por lo que se toma el entero mayor tramos=6 y el tamaño de paso recomendado es al menos 1/6. EL error al aplicar 3 veces la formula es 3(O((1/6)5)) = 0.0003858.

Lo que podría indicar que es necesario al menos dos tramos adicionales con h=1/8 y error O(0,00012) que cumple con el requerimiento.

Se puede aplicar el mismo criterio para Simpson 3/8 y se combinan los errores para verificar que cumplen con el requerimiento.

Algoritmo con Python

Resultados

Ifx:  1.3333332933359998
Igx:  0.238779092876627
Area:  1.094554200459373

medalla area con integral numerico

Instrucciones en Python usando las funciones

# 2Eva_2023PAOI_T1 Material para medalla de academia
import numpy as np
import matplotlib.pyplot as plt

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

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)

# INGRESO
fx = lambda x: 2-8*(0.5-x)**2
gx = lambda x: -(1-x)*np.log(1-x)
a = 0
b = 1-1e-4
muestras1 = 2+1
muestras2 = 3+1

# PROCEDIMIENTO
xi1 = np.linspace(a,b,muestras1)
xi2 = np.linspace(a,b,muestras2)
fi = fx(xi1)
gi = gx(xi2)

Ifx = integrasimpson13_fi(xi1,fi)
Igx = integrasimpson38_fi(xi2,gi)
Area = Ifx - Igx

# SALIDA
print('Ifx: ', Ifx)
print('Igx: ', Igx)
print('Area: ', Area)

plt.plot(xi1,fi,'ob',label='f(x)')
plt.plot(xi2,gi,'or', label='g(x)')
plt.grid()
plt.legend()
plt.xlabel('xi')

# curvas suave con mas muestras (no en evaluación)
xi = np.linspace(a,b,51)
fxi = fx(xi)
gxi = gx(xi)
plt.fill_between(xi,fxi,gxi,color='navajowhite')
plt.plot(xi,fxi,color='blue',linestyle='dotted')
plt.plot(xi,gxi,color='red',linestyle='dotted')

plt.show()

s1Eva_2023PAOI_T3 Recoger los restos de sumergible

ejercicio: 1Eva_2023PAOI_T3 Recoger los restos de sumergible

literal a

Para realizar el planteamiento del ejercicio, se realiza la gráfica para observar mejor los puntos.

recorrido de polinomio

Se podría observar los puntos y plantear un primer recorrido como el siguiente

Rov Sumergible Polinomio01a

los puntos corresponden a los índices j = [0 ,2, 5, 6, 11]:

# INGRESO
xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
               3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
               2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])

# selecionando puntos
j = [0, 2, 5, 6, 11]
xi = xj[j]
fi = fj[j]

literal b

Partiendo del algoritmo de Diferencias divididas de Newton, y añadiendo vectores xj y fj para todos los puntos dados en el ejercicio.

Tabla Diferencia Dividida
[['i   ', 'xi  ', 'fi  ', 'F[1]', 'F[2]', 'F[3]', 'F[4]', 'F[5]']]
[[ 0.      0.2478  1.8108  2.569  -1.3163  0.3389 -0.0561  0.    ]
 [ 1.      1.3802  4.7199 -0.7915 -0.1861  0.09    0.      0.    ]
 [ 2.      2.8008  3.5955 -1.2014  0.111   0.      0.      0.    ]
 [ 3.      3.583   2.6558 -0.9928  0.      0.      0.      0.    ]
 [ 4.      4.6794  1.5673  0.      0.      0.      0.      0.    ]]
dDividida: 
[ 2.569  -1.3163  0.3389 -0.0561  0.    ]
polinomio: 
2.56896856234546*x - 0.0561488275125463*(x - 3.583)*(x - 2.8008)*(x - 1.3802)*(x - 0.2478) 
+ 0.338875729488637*(x - 2.8008)*(x - 1.3802)*(x - 0.2478) 
- 1.31628089036375*(x - 1.3802)*(x - 0.2478) + 1.17420959025079
polinomio simplificado: 
-0.0561488275125463*x**4 + 0.788728905753656*x**3 
- 3.98331084054791*x**2 + 7.41286537452727*x 
+ 0.206696844067185
>>>

usando la tabla de diferencias divididas se plantea la expresión para el polinomio:

p(x) = 1.8108 + 2.569 (x-0.2478) + + (-1.3163)(x-0.2478)(x-1.3802) + + 0.3389 (x-0.2478)(x-1.3802)(x-2.8008) + + (-0.0561)(x-0.2478)(x-1.3802)(x-2.8008)(x-3.583)

el algoritmo simplifica la expresión al siguiente polinomio,

p(x) = -0.0561 x^4 + 0.7887 x^3 - 3.9833 x^2 + 7.4128 x + 0.2066

literal c

Se usa el algoritmo Interpolación polinómica de Lagrange con Python para desarrollar el polinomio a partir de los puntos no usados en el literal anterior.

j = [3,4,7,8,9,10]

Para evitar oscilaciones grandes, se excluye el punto con índice 1, y se obtiene una opción mostrada:

Rov Sumergible Polinomio 02a

con el siguiente resultado:

    valores de fi:  [2.4529 2.3644 4.3657 4.1932 4.6998 4.7536]
divisores en L(i):  [-1.32578996  0.60430667 -0.02841115  0.02632723 -0.09228243  0.13986517]

Polinomio de Lagrange, expresiones
-1.85014223337671*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351) 
+ 3.9125829809921*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.4744) 
- 153.661523787291*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 2.7351)*(x - 2.4744) 
+ 159.272361555669*(x - 4.3757)*(x - 4.2796)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744) 
- 50.928437685792*(x - 4.3757)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744) 
+ 33.9870187307222*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744)

Polinomio de Lagrange: 
-9.26814043907703*x**5 + 163.708008152209*x**4 
- 1143.16428460497*x**3 + 3941.13583467614*x**2 
- 6701.74628786718*x + 4496.64364271972
>>> 

Para presentar la parte analítica puede seleccionar menos puntos, se busca mostrar la aplicación del algoritmo, no necesariamente cuan largo puede ser.

p(x) = + 2.4529\frac{(x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.7351) }{(2.4744 - 4.3757)(2.4744 - 4.2796)(2.4744 - 3.7213)(2.4744 - 3.6627)(2.4744 - 2.7351)} + 2.3644\frac{ (x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.4744) }{(2.7351 - 4.3757)(2.7351 - 4.2796)(2.7351 - 3.7213)(2.7351 - 3.6627)(2.7351 - 2.4744)} +4.3657\frac{(x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 2.7351)(x - 2.4744) }{(3.6627 - 4.3757)(3.6627 - 4.2796)(3.6627 - 3.7213)(3.6627 - 2.7351)(3.6627 - 2.4744)} + 4.1932 \frac{(x - 4.3757)(x - 4.2796)(x - 3.6627)(x - 2.7351)(x - 2.4744) }{(3.7213 - 4.3757)(3.7213 - 4.2796)(3.7213 - 3.6627)(3.7213 - 2.7351)(3.7213 - 2.4744)} +4.6998\frac{(x - 4.3757)(x - 3.7213)(x - 3.6627)(x - 2.7351)(x - 2.4744) }{(4.2796 - 4.3757)(4.2796 - 3.7213)(4.2796 - 3.6627)(4.2796 - 2.7351)(4.2796 - 2.4744)} + 4.7536 \frac{(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.7351)(x - 2.4744)}{(4.3757 - 4.2796)(4.3757 - 3.7213)(4.3757 - 3.6627)(4.3757 - 2.7351)(4.3757 - 2.4744)}

literal d

las gráficas se han presentado en el planteamiento para justificar los criterios usados.

literal e

Los errores de encuentran como la diferencia de los puntos no usados en el polinomio, y evaluando el polinomio en cada coordenada x.

Como las propuestas de polinomio pueden ser variadas se obtendrán diferentes respuestas para cada literal.

Se revisa la aplicación de conceptos en las propuestas.

 

Algoritmos usados

literal b. Diferencias divididas de Newton

# 1Eva_2023PAOI_T3 Recoger los restos de sumergible
# Polinomio interpolación
# Diferencias Divididas de Newton
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
               3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
               2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])

# selecionando puntos
j = [0, 2, 5, 6, 11]
xi = xj[j]
fi = fj[j]


# PROCEDIMIENTO

# Tabla de Diferencias Divididas Avanzadas
titulo = ['i   ','xi  ','fi  ']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias divididas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('F['+str(j-2)+']')

    # cada fila de columna
    i = 0
    paso = j-2 # inicia en 1
    while (i < diagonal):
        denominador = (xi[i+paso]-xi[i])
        numerador = tabla[i+1,j-1]-tabla[i,j-1]
        tabla[i,j] = numerador/denominador
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Divididas
# caso: puntos equidistantes en eje x
dDividida = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    factor = dDividida[j-1]
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

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

# SALIDA
np.set_printoptions(precision = 4)
print('Tabla Diferencia Dividida')
print([titulo])
print(tabla)
print('dDividida: ')
print(dDividida)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xj,fj,'o',color='red')
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')

plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.grid()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencias Divididas - Newton')
plt.show()

literal c. Polinomio de Lagrange

# 1Eva_2023PAOI_T3 Recoger los restos de sumergible
# 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
xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
               3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
               2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])

# selecionando puntos
#j = [0, 2, 5, 6, 11]
j = [3,4,7,8,9,10]
xi = xj[j]
fi = fj[j]

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

# 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(xj,fj,'o',color='red')
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.grid()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.show()

 

s1Eva_2023PAOI_T2 Productos en combo por subida de precio

Ejercicio: 1Eva_2023PAOI_T2 Productos en combo por subida de precio

literal a

Para el ejercicio solo es posible plantear tres ecuaciones, se muestran datos solo para tres semanas. Se tienen 4 incógnitas, por lo que tendría infinitas soluciones y no se podría resolver.

500 a + 600 b + 400 c + 90 d = 1660 800 a + 450 b + 300 c + 100 d = 1825 400 a + 300 b + 600 c + 80 d = 1430

Sin embargo desde el punto de vista práctico se puede usar una variable libre, considerando algunos criterios como:

– que en la nota se da el precio para la docena de huevos a 2.5 USD
– que la variable para huevos es la de menor cantidad se puede incurrir en un menor error si el precio ha variado en los últimos días respecto a lo que se podría haber tenido como referencia.

500 a + 600 b + 400 c = 1660 - 90 d 800 a + 450 b + 300 c = 1825 - 100 d 400 a + 300 b + 600 c = 1430 - 80 d

que al sustituir con d = 2.5

500 a + 600 b + 400 c = 1660 - 90 (2.5) =1435 800 a + 450 b + 300 c = 1825 - 100 (2.5) = 1575 400 a + 300 b + 600 c = 1430 - 80 (2.5) = 1230

Nota: el estudiante puede realizar una justificación diferente para el uso de otra variable libre.

literal b

La forma matricial del sistema de ecuaciones se convierte a:

\begin{pmatrix} 500 & 600 & 400 & \Big| & 1435 \\ 800 & 450 & 300 & \Big| & 1575 \\ 400 & 300 & 600 &\Big| & 1230 \end{pmatrix}

literal c

pivoteando la matriz aumentada, se encuentra que para la primera columna se pivotean la fila 1 y 2, por ser 800 el valor de mayor magnitud:

\begin{pmatrix} 800 & 450 & 300 & \Big| & 1575 \\ 500 & 600 & 400 & \Big| & 1435 \\ 400 & 300 & 600 &\Big| & 1230 \end{pmatrix}

luego, observando la segunda columna desde el elemento de la diagonal hacia abajo, se observa que el mayor valor en magnitud ya se encuentra en la diagonal. No se requieren intercambios de fila para la tercera columna.

literal d

Para el método iterativo de Gauss-Seidel se requiere el vector inicial, tomando las observaciones de la nota. Se da un precio de 50 USD por un quintal métrico de 100Kg, por lo que se estima que el precio por kilogramo ronda 50/100= 0.5 USD. Se indica además que los demás valores se encuentran en el mismo orden de magnitud, por lo que se tiene como vector inicial

X0 = [0.5, 0.5, 0.5]

Para un método iterativo se despejan las ecuaciones como:

a = \frac {1575 - 450 b - 300 c}{800} b = \frac {1435 -500 a - 400 c}{600} c = \frac {1230 - 400 a - 300 b}{600}

Para el cálculo del error se considera que se busca un precio, lo regular es considerar el centavo, es decir una tolerancia 10-2. Para el desarrollo se considera usar cuatro decimales por si se considera truncar o redondear.

itera = 0

a = \frac {1575 - 450 (0.5) - 300 (0.5)}{800} = 1.5 b = \frac {1435 -500(1.5) - 400 (0.5)}{600} = 0.8083 c = \frac {1230 - 400(1.5) - 300 (0.8083)}{600} =0.6458

errado = max|[1.5-0.5, 0.8083-0.5, 0.6458-0.5]|
errado = max|[1, 0.3083, 0.1458]| = 1

itera = 1

X1 = [1.5, 0.8083, 0.6458]

a = \frac {1575 - 450 (0.8083) - 300 (0.6458)}{800} = 1.2719 b = \frac {1435 -500(1.2719) - 400 (0.6458)}{600} = 0.9012 c = \frac {1230 - 400 (1.2719) - 300 (0.9012)}{600} = 0.7514

errado = max|[1.2719-1.5, 0.9012-0.8083, 0.7514-0.6458]|
errado = max|[-0.2281, 0.0929, 0.1056]| = 0.2281

el error disminuye entre iteraciones.

itera = 2

X2 = [1.2719 , 0.9012, 0.7514]

a = \frac {1575 - 450 (0.9012) - 300 (0.7514)}{800} = 1.1805 b = \frac {1435 -500 (1.1805) - 400 (0.7514)}{600} = 0.9069 c = \frac {1230 - 400 (1.1805) - 300 (0.9069)}{600} = 0.8095

errado = max|[ 1.1805-1.2719 , 0.9069-0.9012, 0.8095-0.7514]|
errado = max|[-0.0914, 0.0057, 0.1056]| = 0.1056

el error disminuye entre iteraciones.

literal e

Si el error entre iteraciones disminuye, se considera que el método converge.

usando el algoritmo Método de Gauss-Seidel con Python, se tiene como resultado en 13 iteraciones:

[1.5        0.80833333 0.64583333] 1.0
[1.271875   0.90121528 0.75147569] 0.2281249999999999
[1.18001302 0.90733869 0.80965531] 0.09186197916666661
[1.15475125 0.88960375 0.83536396] 0.025708648304880177
[1.1550864  0.87218536 0.84384972] 0.01741839593330019
[1.16170209 0.86101511 0.84502438] 0.011170246538965367
[1.16754486 0.85536303 0.84395525] 0.005842764350612484
[1.17112508 0.85309227 0.84270381] 0.0035802211655899807
[1.17287167 0.85247107 0.84185002] 0.0017465903731879173
[1.17354127 0.85248226 0.84139802] 0.0006695985845832642
[1.17370447 0.85264759 0.84120656] 0.00019146597344232852
[1.17368327 0.8527929  0.84114804] 0.00014530950399282982
[1.17362348 0.85288174 0.84114348] 8.884049018109685e-05
respuesta X: 
[[1.17362348]
 [0.85288174]
 [0.84114348]]
verificar A.X=B: 
[[1575.03861029]
 [1434.99817609]
 [1230.        ]]
iteraciones 13
>>> 

instrucciones en Python

Se debe recordar que las matrices y vectores deben ser tipo real (float) .

# Método de Gauss-Seidel
# solución de sistemas de ecuaciones
# por métodos iterativos

import numpy as np

# INGRESO
A = np.array([[800, 450, 300],
              [500, 600, 400],
              [400, 300, 600]], dtype=float)

B = np.array([1575, 1435,1230], dtype=float)

X0  = np.array([0.5,0.5,0.5])

tolera = 0.0001
iteramax = 100

# PROCEDIMIENTO

# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
#  valores iniciales
X = np.copy(X0)
diferencia = np.ones(n, dtype=float)
errado = 2*tolera

itera = 0
while not(errado<=tolera or itera>iteramax):
    # por fila
    for i in range(0,n,1):
        # por columna
        suma = 0 
        for j in range(0,m,1):
            # excepto diagonal de A
            if (i!=j): 
                suma = suma-A[i,j]*X[j]
        
        nuevo = (B[i]+suma)/A[i,i]
        diferencia[i] = np.abs(nuevo-X[i])
        X[i] = nuevo
    errado = np.max(diferencia)
    print(X,errado)
    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)
print('iteraciones',itera)

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