Ejercicio: 1Eva_2022PAOII_T1 Esfera flotando en agua

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

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

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

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

usando los valores dados para el ejercicio, r=1 y ρ_esfera = 200 Kg/m³ y ρ_agua    = 1000 kg/m³ se tiene que:

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

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

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

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

con lo que se puede plantear dos ecuaciones al despejar h

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

Iteraciones de la primera ecuación

itera = 0 ; h = h₀ = 0.5 ;

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

itera = 1 ; h = 1.1313 ;

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

itera = 2 ; h = 1.3086 ;

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

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

---

resultados.txt

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

Algoritmo en Python

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

import numpy as np
import matplotlib.pyplot as plt

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

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

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

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

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

# SALIDA
print('raiz:',respuesta)

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

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

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

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

literal a. gráfica de malla

literal b. Ecuaciones de diferencias divididas a usar

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

se agrupan las constantes,

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

literal d Determine el valor de λ

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

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

continuando luego con la ecuación general,

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

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

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

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

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

literal e iteraciones

iteración  i=1, j=0

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

iteración  i=2, j=0

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

iteración  i=3, j=0

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

literal f

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

O(h_x²) + O(h_t) = (1/3)² + 0.02 = 0,1311

literal g

Resultados usando el algoritmo en Python

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

Instrucciones en Python

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

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

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

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

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

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

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

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

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

Ejercicio: 2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

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

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

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

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

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

literal a

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

y'' = - 2y' - 5 y

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

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

se usa una tabla para llevar el registro de operaciones:

Se plantea las operaciones:

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

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

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

---

literal b

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

itera = 0

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

itera = 1

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

itera = 2

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

Resultados con el algoritmo en Python

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

---

Literal c

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

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

---

Literal d

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

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,7),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
    return(estimado)

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

t0 = 0
y0 = 2
z0 = -4

h = 0.01
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=4)
print('ti, yi, zi, K1y, K1z, K2y, K2z')
print(tabla)

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

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

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

Literal a. Integral con Simpson 1/3

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

El tamaño de paso se calcula como:

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

representados en una gráfica como

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

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

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

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

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

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

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

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

Literal b. Integral con Cuadratura de Gauss

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

xj= [0. 1.5 3. ]

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

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

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

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

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

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

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

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

Instrucciones integradas en Python

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

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

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

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

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

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

Ag = Ag1 + Ag2

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

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

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

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