Categoría: Soluciones

Ejercicio: 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada

literal a

El planteamiento de las ecuaciones según se cita en la solución iterativa es el promedio de los nodos adjacentes:

a = \frac{100+40+b+c}{4} b = \frac{40+20+a+d}{4} c = \frac{100+80+a+d}{4} d = \frac{80+20+c+b}{4}

reordenando las ecuaciones para su forma maticial:

4a -b -c = 100+40 4b -a -d= 40+20 4c -a -d= 100+80 4d -c -b= 80+20

literal b

la forma matricial del sistema de ecuaciones es:

\begin{bmatrix} 4 & -1 &-1 & 0 \\ -1 & 4 & 0 &-1 \\ -1 & 0 & 4 &-1 \\ 0 & -1 &-1 & 4 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 140 \\ 60 \\ 180 \\ 100 \end{bmatrix}

matriz aumentada:

\left[ \begin{array} {cccc|c} 4 & -1 &-1 & 0 & 140\\ -1 & 4 & 0 &-1 & 60 \\ -1 & 0 & 4 &-1 & 180 \\ 0 & -1 &-1 & 4 & 100 \end{array}\right]

al revisar el procedimiento de pivoteo parcial por filas se muestra que ya se encuentra pivoteada la matriz. Por lo que no es necesario realizar cambios en las filas.

literal c

Las expresiones a partir de la matriz aumentada y pivoteada para usar en un método iterativo como Gauss-Seidel son:

a = 0.25(b+c+140) b = 0.25*(a+d+60) c = 0.25(a+d+180) d = 0.25(c+b+100)

el vector inicial debe considerar que las temperaturas serán medias entre los  bordes de la placa

X₀ = [ (100+40)/2, (40+20)/2, (100+80)/2, (80+20)/2 ]

X₀ = [ 70, 30, 90, 50 ]

literal d iteraciones

itera = 0

a = 0.25(30+90+140)= 65

b = 0.25*(65+50+60) = 43.75

c = 0.25(65+50+180) = 73.75

d = 0.25(73.75+43.75+100) = 54.375

error = max|[65-70, 43.75-30, 73.75-90, 54.375-50]|

error = max|[ -5, 13.75, -16.25, 4.375| = 16.25

X₁ = [ 65, 43.75, 73.75, 54.375 ]

itera = 1

a = 0.25(43.75+73.75+140)= 64.375

b = 0.25*(64.375+54.375+60) = 44.687

c = 0.25(64.375+54.375+180) = 74.687

d = 0.25(74.687+44.687+100) = 54.843

error = max|[64.375-65, 44.687-43.75, 74.687-73.75, 54.843-54.375]|

error = max|[-0.625, 0.937, 0.937, 0.468| = 0.937

X₂ = [ 64.375, 44.687, 74.687, 54.843 ]

itera = 2

a = 0.25(44.687+74.687+140)= 64.843

b = 0.25*(64.843 +54.843+60) = 44.921

c = 0.25(64.843+54.843+180) = 74.921

d = 0.25(74.921+44.921+100) = 54.960

error = max|[64.843-64.375, 44.921-44.687, 74.921-74.687, 54.960-54.843]|

error = max|[0.468, 0.234, 0.234, 0.117| = 0.468

X₃ = [ 64.843, 44.921, 74.921, 54.960 ]

literal e

El método converge pues los errores en cada iteración disminuyen. Los valores obtenidos en el resultado son acorde a las temperaturas esperadas en los nodos.

Algoritmos con Python

Los resultados obtenidos son:

X 0 = [65.    43.75  73.75  54.375]
X 1 = [64.375   44.6875  74.6875  54.84375]
X 2 = [64.84375   44.921875  74.921875  54.9609375]
X 3 = [64.9609375  44.98046875 74.98046875 54.99023438]
X 4 = [64.99023438 44.99511719 74.99511719 54.99755859]
X 5 = [64.99755859 44.9987793  74.9987793  54.99938965]
X 6 = [64.99938965 44.99969482 74.99969482 54.99984741]
X 7 = [64.99984741 44.99992371 74.99992371 54.99996185]
X 8 = [64.99996185 44.99998093 74.99998093 54.99999046]
X 9 = [64.99999046 44.99999523 74.99999523 54.99999762]
respuesta X: 
[[64.99999046]
 [44.99999523]
 [74.99999523]
 [54.99999762]]
verificar A.X=B: 
[[139.99997139]
 [ 59.99999285]
 [179.99999285]
 [100.        ]]
>>>

el algoritmo usado es:

# 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada
# Método de Gauss-Seidel
import numpy as np
# INGRESO
A = np.array([[4,-1,-1,0],
              [-1,4,0,-1],
              [-1,0,4,-1],
              [0,-1,-1,4]],dtype=float)
B = np.array([140,60,180,100],dtype=float)
X0  = np.array([70,30,90,50],dtype=float)

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
    print('X',itera,'=',X)
    errado = np.max(diferencia)
    itera = itera + 1

# Respuesta X en columna
X = np.transpose([X])
# revisa si NO converge
if (itera>iteramax):
    X=0
# revisa respuesta
verifica = np.dot(A,X)

# SALIDA
print('respuesta X: ')
print(X)
print('verificar A.X=B: ')
print(verifica)

 

Ejercicio: 1Eva_2024PAOI_T1 Vaciado de reservorio semiesférico

literal a. Planteamiento

A partir de la solución expresión propuesta y simplificada del ejercicio, se reemplaza los valores de R, K, a y g que son constantes. Como se da el valor de t=(15 min)(60 s/min), se unifica la unidad de medida de tiempo a segundos pues g=9.6 m/s².

---

-\frac{4}{3}R h^{3/2} + \frac{2}{5} h^{5/2} = -\frac{Ka}{\pi} t \sqrt{2g}

---

-\frac{4}{3}(3) h^{3/2} + \frac{2}{5} h^{5/2} = -\frac{(0.85)(0.01)}{\pi} (15)(60) \sqrt{2(9.8)}

la función para encontrar la raíz se reordena como f(h)=0

f(h) = -\frac{4}{3}(3) h^{3/2} + \frac{2}{5} h^{5/2} +\frac{(0.85)(0.01)}{\pi} (15)(60) \sqrt{2(9.8)} f(h) = -4h^{3/2} + 0.4 h^{5/2} +10.7805

Usando la gráfica proporcionada del reservorio semiesférico, el valor de h no puede superar la altura R. Al vaciar el reservorio h=0.

Por lo que el intervalo a usar es [0,3]. Se verifica la existencia de una raíz al con la gráfica y Python.

b. Método de Newton Raphson

El método requiere la derivada f'(h)

f(h) = -4h^{3/2} + 0.4 h^{5/2} +10.7805 f'(h) = -6\frac{ h^{3/2}}{h} + \frac{h^{5/2}}{h} f'(h) = -6 h^{1/2}+ h^{3/2}

El valor inicial puede ser por ejemplo, el centro del intervalo:

x_0= \frac{a+b}{2} = \frac{0+3}{2} = 1.5

itera=0

f(1.5) = -4(1.5)^{3/2} + 0.4 (1.5)^{5/2} +10.7805= 4.5343 f'(1.5) = -6\frac{ 1.5^{3/2}}{1.5} + \frac{1.5^{5/2}}{1.5} = -5.5113 x_1= x_0 - \frac{f(1.5)}{f'(1.5)}= 1.5-\frac{4.5343}{-5.5113} =2.3227 error_1 =|x1-x0| = |2.3227-1.5| =0.8227

itera=1

f(2.3227) = -4(2.3227)^{3/2} + 0.4 (2.3227)^{5/2} +10.7805= -0.09019 f'(2.3227) = -6\frac{ 2.3227^{3/2}}{2.3227} + \frac{2.3227^{5/2}}{2.3227} = -5.6043 x_2= 2.3227-\frac{-0.09019}{-5.6043} = 2.3066 error_2 = |2.3066-2.3227| =0.01611

itera=2

f(2.3227) = -4(2.3066)^{3/2} + 0.4 (2.3066)^{5/2} +10.7805= -0.0007104 f'(2.3227) = -6\frac{ 2.3227^{3/2}}{2.3227} + \frac{2.3227^{5/2}}{2.3227} = -5.6093 x_3= 2.3227-\frac{-0.0007104}{-5.6093} = 2.3066 error_3 = |2.3066-2.3066| =0.000007241

literal c, convergencia

considerando tolerancia de 10-3 (milimétrica), solo serán necesarias 3 iteraciones.

El método converge al mostrarse que los errores disminuyen en cada iteración.

El resultado se encuentra dentro del intervalo y es x = 2.3066

Algoritmo con Python

resultados:

fh:
         ____          ____                   
        /  3          /  5                    
- 4.0*\/  h   + 0.4*\/  h   + 10.7805172327811
dfh:
         ____          ____
        /  3          /  5 
  6.0*\/  h     1.0*\/  h  
- ----------- + -----------
       h             h     
['xi', 'xnuevo', 'tramo']
[[1.5000e+00 2.3227e+00 8.2272e-01]
 [2.3227e+00 2.3066e+00 1.6118e-02]
 [2.3066e+00 2.3066e+00 7.2412e-06]]
raiz en:  2.306612665792577
con error de:  7.241218404896443e-06

>>> f(1.5)
4.534318388683996
>>> df(1.5)
-5.5113519212621505
>>> f(2.3227)
-0.09019959114247378
>>> df(2.3227)
-5.604354799446854
>>> f(2.3066)
7.104683901282272e-05
>>> df(2.3066)
-5.609349350102558

instrucciones:

# 1Eva_2024PAOI_T1 Vaciado de reservorio semiesférico
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO
h = sym.Symbol('h')
t = 15*60
R = 3 
g = 9.8 
area = 0.01
K = 0.85 
fh = -(4/3)*R*sym.sqrt(h**3)+(2/5)*sym.sqrt(h**5)+(K*area/np.pi)*np.sqrt(2*g)*(t)

# Grafica de f(h)
a = 0
b = 3
muestras = 15

# PROCEDIMIENTO
hi = np.linspace(a,b,muestras)
# para evaluación numérica con numpy
f = sym.lambdify(h,fh)
fi = f(hi)

# derivada con sympy
dfh = sym.diff(fh,h,1)
df = sym.lambdify(h,dfh)

# SALIDA
print('fh:')
sym.pprint(fh)
print('dfh:')
sym.pprint(dfh)

plt.plot(hi,fi)
plt.axhline(0, color='black')
plt.xlabel('h')
plt.ylabel('f')
plt.title(fh)
plt.show()

## Método de Newton-Raphson

# INGRESO
fx  = f
dfx = df
x0 = (a+b)/2 # mitad del intervalo (ejemplo)
tolera = 0.001

# PROCEDIMIENTO
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(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)

 

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

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

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

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

Para el intervalo [0,3]

La función es una constante

literal a, b y c

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

Para el intervalo [3,5]

literal a

La función es una recta con pendiente negativa

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

para el integral con cuadratura de Gauss

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

literal b y c

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

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

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

Para el intervalo [5, 14.97]

literal a

Del polinomio obtenido en el ejercicio anterior para éste intervalo

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

literal b y c

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

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

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

literal d

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

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

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

literal e

Resultados con el algoritmo

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

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

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

Instrucciones en Python

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

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

xmuestras = 31

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

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

Añadir para graficar en 2D

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

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

Añadir para graficar en 3D

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

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

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

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

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

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

.

Ejercicio: 3Eva_2023PAOII_T2 perfil de un peón

literal a y b

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

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

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

p_1(x) = 15

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

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

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

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

simplificando la expresión con Python y Sympy:

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

literal c

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

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

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

literal d

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

Resultados para el intervalo [5, 14.97]

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

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

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

Instrucciones Python

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

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


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

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

# PROCEDIMIENTO
# Polinomio de Lagrange
n = len(xi)
x = sym.Symbol('x')
polinomio = 0
divisorL = np.zeros(n, dtype = float)
for i in range(0,n,1):
    
    # Termino de Lagrange
    numerador = 1
    denominador = 1
    for j  in range(0,n,1):
        if (j!=i):
            numerador = numerador*(x-xi[j])
            denominador = denominador*(xi[i]-xi[j])
    terminoLi = numerador/denominador

    polinomio = polinomio + terminoLi*fi[i]
    divisorL[i] = denominador

# simplifica el polinomio
polisimple = polinomio.expand()

# para evaluación numérica
px = sym.lambdify(x,polisimple)

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

# SALIDA
print('    valores de fi: ',fi)
print('divisores en L(i): ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(polisimple)

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

 

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

literal a y b

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

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

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

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

para la segunda:

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

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

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

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

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

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

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

Instrucciones en Python

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

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

a = 0.5
b = 16
muestras = 21

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

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

literal c

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

Para Newton-Raphson se usa:

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

(derivada obtenida con sympy)

itera = 0 ; xi = 3

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

itera = 1 ; xi = 4.55

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

itera = 2 ; xi = 4.98

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

literal d

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

Con el algoritmo se muestra que converge a x=5

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

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

Instrucciones en Python

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

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

x0 = 3
tolera = 0.0001 # 1e-4

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

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

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

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

Ejercicio: 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 w₀ = 1 000 lbs/ft, T₀. = 0.588345×10⁶.
dy(0)/dx = 0 y l_B=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: x₀ = 0, y₀=0, z₀ = 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

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

 

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 f_v:

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 h_f: 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):

x_is= [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(h⁵) = O(0.5666⁵) = 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 g_v:

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 h_g:

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

x_ic = [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

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

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

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 Δy²

\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])

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/s², 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]
...

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

Ejercicio: 2Eva_2023PAOI_T1 Material para medalla de academia

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.5⁵)= 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(h⁵) = 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)⁵)) = 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

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