s3Eva_2023PAOII_T3 Volumen por solido de revolución de un peón

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

volumen de un Peon 2D de revolucion

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

un peon 3D

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

.

s3Eva_2023PAOII_T2 perfil de un peón

Ejercicio: 3Eva_2023PAOII_T2 perfil de un peón

literal a y b

perfil de un peon 01Para 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

interpola peon 01

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

 

s3Eva_2023PAOII_T1 Intersección con círculo

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:

interseccion con circulo

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)

s2Eva_2023PAOII_T2 Cable cuelga entre apoyos A y B

Ejercicio: 2Eva_2023PAOII_T2 Cable cuelga entre apoyos A y B

Literal a

La ecuación diferencial a resolver es:

\frac{d^2y}{dx^2} = \frac{w_0}{T_0} \Big[ 1+ \sin \Big(\frac{\pi x}{2l_B} \Big) \Big]

donde w0 = 1 000 lbs/ft, T0. = 0.588345×106.
dy(0)/dx = 0 y lB=200 de la gráfica presentada.

\frac{d^2y}{dx^2} = \frac{1000}{0.588345×10^6} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big]

Para usar Runge Kutta para segunda derivada:

z= y' = f(x,y,z) z' = (y')' = 0z + \frac{1000}{0.588345×10^6} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big] g(x,y,z) = \frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi x}{2(200)} \Big) \Big]

los valores iniciales para el ejercicio acorde al enunciado son: x0 = 0, y0=0, z0 = 0, con h=0.5

literal b

para itera 0

K1y = h*z = 0.5*0 = 0 K1z = (0.5)\frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi (0)}{2(200)} \Big) \Big] = 0.0008498 K2y = h*(z+K1z) = (0.5) (0+0.00084984) = 0.0004249 K2z = (0.5)\frac{1}{0.588345×10^3} \Big[ 1+ \sin \Big(\frac{\pi (0+0.5)}{2(200)} \Big) \Big] =0.0008531 y = 0+\frac{0+0.0004249}{2} = 0.0002124 z = 0+\frac{0.0008498+0.0008531}{2} = 0.0008515 x = 0 + 0.5 = 0.5

Desarrollar dos iteraciones adicionales como tarea.

Para las primeras iteraciones de un total de 400+1, los valores con Python y en resultados.txt :

estimado[xi,yi,zi,K1y,K2y,K1z,K2z]
[0.0000 0.0000e+00 0.0000e+00 
0.0000e+00 0.0000e+00 
0.0000e+00 0.0000e+00]

[0.5000 2.124603761398499073e-04 8.515101601627471980e-04 
0.000000000000000000e+00 4.249207522796998146e-04 
8.498415045593996292e-04 8.531788157660947667e-04]

[1.0000 8.515101601627470896e-04 1.706357605799455881e-03 
4.257550800813735990e-04 8.523444879644209281e-04 
8.531788157660947667e-04 8.565160755073224913e-04]

[1.5000 1.918817981939305653e-03 2.564542259712321911e-03 
8.531788028997279406e-04 1.281436840653389295e-03 
8.565160755073224913e-04 8.598532323184093513e-04]

[2.0000 3.416052419875068892e-03 3.426063993239661116e-03 
1.282271129856160955e-03 1.712197746015365740e-03 
8.598532323184093513e-04 8.631902347362692754e-04]
...

literal c

resultado en archivo.txt al ejecutar el algoritmo.

literal d

cable entre apoyos A y B

Algoritmo con Python

# 2Eva_2023PAOII_T2 Cable cuelga entre apoyos A y B
import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3+4),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,K2y,K1z,K2z]
    return(estimado)

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

# INGRESO
T0 = 0.588345e6
LB = 200
f = lambda x,y,z: z
g = lambda x,y,z: (1000/T0)*(1+np.sin(np.pi*x/(2*LB)))
x0 = 0
y0 = 0
z0 = 0
h  = 0.5
muestras = 401

# PROCEDIMIENTO
puntosRK2 = rungekutta2_fg(f,g,x0,y0,z0,h,muestras)
xi = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# SALIDA
np.set_printoptions(precision=4)
print('estimado[xi,yi,zi,K1y,K2y,K1z,K2z]')
print(puntosRK2)
np.savetxt("tablaRk2.txt",puntosRK2)


# Gráfica
import matplotlib.pyplot as plt

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

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

 

s2Eva_2023PAOII_T1 Volumen por solido de revolución

Ejercicio: 2Eva_2023PAOII_T1 Volumen por solido de revolución

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

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

literal a y c

Para el volumen con f(x) con al menos 3 tramos y un método de Simpson, directamente se puede usar 3/8. Por lo que se Se reemplaza en la fórmula de volumen del sólido de revolución f(x) con:

f(x) = \sqrt{\sin (x/2)}

obteniendo:

V_{fx} = \int_{a}^{b} \pi \Big(\sqrt{\sin (x/2)} \Big)^2 dx = \int_{a}^{b} \pi \sin (x/2) dx

La expresión dentro del integral se denomina como fv:

f_v (x)= \pi \sin (x/2)

en el intervalo [0.1, 1.8],  con al menos 3 tramos, se requieren 4 muestras con tamaño de paso hf: y truncando a 4 decimales los resultados calculados con Python.

h_f =\frac{b-a}{tramos} = \frac{1.8-0.1}{3} = 0.5666

los puntos de muestra quedan np.linspace(0.1,1.8,3+1):

xis= [0.1, 0.6666, 1.2333, 1.8 ]

El integral se calcula con los puntos de muestra,

V_{fx} = \frac{3}{8} (0.5666) \Big( f_v(0.1) +3 f_v(0.6666) + + 3 f_v(1.2333)+ f_v(1.8)\Big)

recordando que se usa en radianes,

V_{fx} = \frac{3}{8} (0.5666) \Bigg( \pi \sin \Big(\frac{0.1}{2}\Big) +3 \pi \sin \Big(\frac{0.6666}{2}\Big) + + 3 \pi \sin\Big(\frac{1.2333}{2}\Big)+ \pi \sin \Big(\frac{1.8}{2}\Big)\Bigg) = \frac{3}{8} (0.5666) \Big( 0.1570+3 (1.0279) + + 3 (1.8168)+ 2.4608\Big)

literal d. el volumen generado por f(x) tiene como resultado:

V_{fx} = 2.3698

la cota de error para fx es el orden de O(h5) = O(0.56665) = O(0.05843), queda como tarea completar la cota de error total.

literal b y c

Para el volumen con g(x) con al menos 2 tramos y Cuadratura de Gauss de dos puntos, se reemplaza en la fórmula de volumen de sólido de revolución:

g(x) = e^{x/3} - 1 V_{gx} = \int_{a}^{b} \pi (e^{x/3} - 1)^2 dx

La expresión dentro del integral se denomina como gv:

g_v = \pi (e^{x/3} - 1)^2

en el intervalo [0.1, 1.8],  con al menos 2 tramos, se requieren 3 muestras con tamaño de paso hg:

h_g =\frac{b-a}{tramos} = \frac{1.8-0.1}{2} = 0.85

xic = [0.1, 0.95, 1.8 ]

tramo 1: [0.1, 0.95] , a = 0.1 , b= 0.95, truncando a 4 decimales

x_a = \frac{0.95+0.1}{2} - \frac{0.95-0.1}{2}\frac{1}{\sqrt{3}} = 0.2796 x_b = \frac{0.95+0.1}{2} + \frac{0.95-0.1}{2}\frac{1}{\sqrt{3}} = 0.7703 g_v(0.2796) = \pi (e^{0.2796/3} - 1)^2 = 0.02998 g_v(0.7703) = \pi (e^{0.7703/3} - 1)^2 = 0.2692 V_{c1} = \frac{0.95-0.1}{2}(g_v(0.2796) + g_v(0.7703)) V_{c1} = \frac{0.95-0.1}{2}(0.02998 + 0.2692) V_{c1} = 0.1271

tramo 2: [0.95, 1.8] , a = 0.95 , b= 1.8

x_a = \frac{1.8+0.95}{2} - \frac{1.8-0.95}{2}\frac{1}{\sqrt{3}} = 1.1296 x_b = \frac{1.8+0.95}{2} - \frac{1.8-0.95}{2}\frac{1}{\sqrt{3}} = 1.6203 g_v(1.1296) = \pi (e^{1.1296/3} - 1)^2 = 0.6567 g_v(1.6203) = \pi (e^{1.6203/3} - 1)^2 = 1.6115 V_{c2} = \frac{1.8-0.95}{2}(g_v(1.1296) + g_v(1.6203)) V_{c2} = \frac{1.8-0.95}{2}(0.6567 + 1.6115) V_{c2} = 0.9640

literal d. volumen generado por g(x)

V_{gx} = V_{c1} + V_{c2} = 0.1271 + 0.9640 = 1.0912

completar la cota de error para cuadratura de Gauss de dos puntos.

literal e. El volumen de revolución se genera como la resta del volumen de f(x) y volumen g(x)

V = V_{fx} - V_{gx} = 2.3698 - 1.0912 = 1.2785

Algoritmo con Python

Los resultados usando el algoritmo con las operaciones usadas en el planteamiento son:

para f(x):
xis= [0.1        0.66666667 1.23333333 1.8       ]
fiv= [0.15701419 1.02791246 1.81684275 2.46089406]
Volumenfx:  2.369836948864926

para g(x):
Por tramos: [0.1  0.95 1.8 ]
xab= [0.2796261355944091, 0.770373864405591,
      1.129626135594409, 1.620373864405591]
gab= [0.02998177327492598, 0.26928904479784566,
      0.6567986343358181, 1.6115494735531555]
Vc1= 0.12719009768092793  ; Vc2= 0.964047945852814
Volumengx:  1.0912380435337419

Volumen solido revolucion: 1.2785989053311841

Considerando realizar los cálculos para cada sección:

# 2Eva_2023PAOII_T1 Volumen por solido de revolución
import numpy as np

# INGRESO
fx = lambda x: np.sqrt(np.sin(x/2))
gx = lambda x: np.exp(x/3)-1
a = 0.1
b = 1.8
tramosSimpson = 3
tramosCGauss = 2

# PROCEDIMIENTO
# Volumen para f(x) con Simpson
fv = lambda x: np.pi*np.sin(x/2)
hs = (b-a)/tramosSimpson
xis = np.linspace(a,b,tramosSimpson +1)
fiv = fv(xis)
Vs = (3/8)*hs*(fiv[0]+3*fiv[1]+3*fiv[2]+ fiv[3])

# Volumen para g(x) con Cuadratura de Gauss
gv = lambda x: np.pi*(np.exp(x/3)-1)**2
hc = (b-a)/tramosSimpson
xic = np.linspace(a,b,tramosCGauss +1)
# tramo 1
ac = xic[0]
bc = xic[1]
xa = (bc+ac)/2 + (bc-ac)/2*(-1/np.sqrt(3)) 
xb = (bc+ac)/2 + (bc-ac)/2*(1/np.sqrt(3))
Vc1 = (bc-ac)/2*(gv(xa)+gv(xb))
xab = [xa,xb]
gab = [gv(xa),gv(xb)]
# tramo 2
ac = xic[1]
bc = xic[2]
xa = (bc+ac)/2 + (bc-ac)/2*(-1/np.sqrt(3)) 
xb = (bc+ac)/2 + (bc-ac)/2*(1/np.sqrt(3))
Vc2 = (bc-ac)/2*(gv(xa)+gv(xb))
Vc = Vc1+Vc2
xab.append(xa)
xab.append(xb)
gab.append(gv(xa))
gab.append(gv(xb))

# Volumen solido revolucion
Volumen = Vs - Vc

# SALIDA
print("para f(x):")
print("xis=", xis)
print("fiv=", fiv)
print("Volumenfx: ",Vs)
print()
print("para g(x):")
print("Por tramos:",xic)
print("xab=", xab)
print("gab=", gab)
print("Vc1=",Vc1," ; Vc2=",Vc2) 
print("Volumengx: ",Vc)
print()
print("Volumen solido revolucion:",Volumen)

para la gráfica presentada en el enunciado (no requerida) , se complementa con las instrucciones:

# para grafica -------------------
import matplotlib.pyplot as plt
muestras = 21 # grafica
xi = np.linspace(a,b,muestras)
fi = fx(xi)
gi = gx(xi)
xig = np.linspace(a,b,tramosCGauss+1)
fis = fx(xis)
gig = gx(xig)

# grafica
plt.plot(xi,fi, label="f(x)")
plt.plot(xi,gi, label="g(x)")
plt.plot([0.0,2.0],[0,0], marker=".", color="blue")
plt.fill_between(xi,fi,gi,color="lightgreen")
plt.axhline(0)
plt.axvline(a, linestyle="dashed")
plt.axvline(b, linestyle="dashed")
plt.xlabel('x')
plt.ylabel('f(x), g(x)')
plt.legend()
plt.plot(xis,fis,'.b')
plt.plot(xig,gig,'.r')
plt.grid()
plt.show()

Gráfica de sólido de revolución en 3D

sólido de revolución 3D

Instrucciones en Python

# 2Eva_2023PAOII_T1 Volumen por solido de revolución
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d

# INGRESO
f = lambda x: np.sqrt(np.sin(x/2))
g = lambda x: np.exp(x/3)-1

# eje x
xa = 0.1
xb = 1.8
xmuestras = 31
# angulo w de rotación
w_a = 0
w_b = 2*np.pi
w_muestras = 31

# PROCEDIMIENTO
# muestreo en x y angulo w
xi = np.linspace(xa, xb, xmuestras)
wi = np.linspace(w_a, w_b, w_muestras)
X, W = np.meshgrid(xi, wi)

# evalua f(x) en 3D
Yf = f(xi)*np.cos(W)
Zf = f(xi)*np.sin(W)

# evalua g(x) en 3D
Yg = g(xi)*np.cos(W)
Zg = g(xi)*np.sin(W)

# SALIDA

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

grafica.plot_surface(X, Yf, Zf,
                     color='blue', label='f(x)',
                     alpha=0.3, rstride=6, cstride=12)
grafica.plot_surface(X, Yg, Zg,
                     color='orange', label='g(x)',
                     alpha=0.3, rstride=6, cstride=12)

grafica.set_title('Solidos de revolución')
grafica.set_xlabel('x')
grafica.set_ylabel('y')
grafica.set_zlabel('z')
# grafica.legend()
eleva = 45
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()

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)