Autor: Edison Del Rosario

Ejercicio: 2Eva_IIT2019_T4 Integrar con Cuadratura de Gauss

f(x) = x ln(x)

1 ≤x≤4

se requiere:

I = \int_1^4 x ln(x) dx

literal a. Usando el método de Cuadratura de Gauss con 2 términos

x_a = \frac{b+a}{2} + \frac{b-a}{2}x_0 = \frac{4+1}{2} + \frac{4-1}{2}\Big(\frac{-1}{\sqrt{3}} \Big)

x_a =1.6339745962155612

x_b = \frac{b+a}{2} + \frac{b-a}{2}x_1 = \frac{4+1}{2} + \frac{4-1}{2}\Big(\frac{1}{\sqrt{3}} \Big)

x_b =3.366025403784439

I \cong \frac{b-a}{2}(f(x_a) + f(x_b)) I \cong \frac{4-1}{2}(x_a ln(x_a) + x_b ln(x_b))

I = 7.33164251999249

literal b.  De la fórmula , despejar el valor del error<0.0001

\Big|\frac{(b-a)}{180}h^4 f^{(4)} (\xi)\Big| <0.0001; \xi \in[a,b] h^4 <0.0001\frac{180}{(4-1)}\frac{1}{f^{(4)} (\xi)} h^4 < 0.006\frac{1}{f^{(4)} (\xi)} h <\Big(0.006\frac{1}{f^{(4)} (\xi)}\Big)^{1/4}

obteniendo la 4ta derivada de la función:

f(x) = x ln(x) f'(x) = ln(x) + x\Big(\frac{1}{x} \Big) = ln(x) +1 f''(x) = \frac{1}{x} f'''(x) = -\frac{1}{x^2} f^{(4)}(x) = 2\frac{1}{x^3}

se tiene que:

h <\Big(0.006\frac{1}{f^{(4)} (\xi)}\Big)^{1/4} h <\Big(0.006\frac{1}{2\frac{1}{\xi^3}}\Big)^{1/4} h <\Big(0.003\xi^3\Big)^{1/4} h <(0.003)^{1/4}\xi^{3/4}

en el peor de los casos, se toma el valor menor de ξ =1

h <(0.003)^{1/4} h<0.2340347319320716

 

Ejercicio: 2Eva_IIT2019_T3 EDP elíptica, placa en (1,1)

dada la ecuación del problema:

\frac{\delta ^2 u}{\delta x^2} + \frac{\delta ^2 u}{\delta y^2} = \frac{x}{y} + \frac{y}{x}

1 <  x < 2
1 <  y < 2

Se convierte a la versión discreta usando diferencias divididas centradas:

---

\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} = \frac{x_i}{y_j} + \frac{y_j}{x_i}

---

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) = =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)

---

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

---

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

---

u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)

---

Iteraciones

que permite plantear las ecuaciones para cada punto en posición [i,j]

---

i=1, j=1

u[0,1]-4u[1,1]+u[2,1] + + u[1,0]+u[1,2] =(0.25)^2\Big( \frac{0.25}{0.25} + \frac{0.25}{0.25}\Big) 0.25ln(0.25)-4u[1,1]+u[2,1] + + 0.25ln(0.25)+u[1,2] = 0.125 -4u[1,1]+u[2,1] +u[1,2] = 0.125 - 2(0.25)ln(0.25) -4u[1,1]+u[2,1] +u[1,2] = 0.8181

---

i=2, j=1

u[1,1]-4u[2,1]+u[3,1]+ +0.5 ln(0.5)+u[2,2]=0.15625 u[1,1]-4u[2,1]+u[3,1]+u[2,2]=0.15625-0.5 ln(0.5) u[1,1]-4u[2,1]+u[3,1]+u[2,2]=0.8493

---

i=3, j=1

u[2,1]-4u[3,1]+u[4,1] + + u[3,0]+u[3,2] =(0.25)^2\Big( \frac{0.75}{0.25} + \frac{0.25}{0.75}\Big) u[2,1]-4u[3,1]+u[4,1]+u[3,2] = 0.20833333-0.75 ln(0.75) u[2,1]-4u[3,1]+u[4,1]+u[3,2] =0.4240

---

Tarea: continuar con el ejercicio hasta plantear todo el sistema de ecuaciones.

A = np.array([
[-4, 1, 0, 1, 0, 0, 0, 0, 0],
[ 1,-4, 1, 0, 1, 0, 0, 0, 0],
[ 0, 1,-4, 0, 0, 1, 0, 0, 0],
[ 1, 0, 0,-4, 1, 0, 1, 0, 0],
[ 0, 1, 0, 1,-4, 1, 0, 1, 0],
[ 0, 0, 1, 0, 1,-4, 0, 0, 1],
[ 0, 0, 0, 1, 0, 0,-4, 1, 0],
[ 0, 0, 0, 0, 1, 0, 1,-4, 1],
[ 0, 0, 0, 0, 0, 1, 0, 1,-4]])
B = np.array(
[0.125 - 2(0.25)ln(0.25),
 0.15625 - 0.5ln(0.5),
 0.20833 - 0.75 ln(0.75),
 ...])

B = [0.8181,
     0.8493,
     0.4240,
     ...]

---

Algoritmo con Python

Con valores para la matriz solución:

iteraciones:  15
error entre iteraciones:  6.772297286980838e-05
solución para u: 
[[0.        0.2789294 0.6081976 0.9793276 1.3862943]
 [0.2789294 0.6978116 1.1792239 1.7127402 2.2907268]
 [0.6081976 1.1792239 1.8252746 2.5338403 3.2958368]
 [0.9793276 1.7127402 2.5338403 3.4280053 4.3846703]
 [1.3862943 2.2907268 3.2958368 4.3846703 5.5451774]]
>>>

y algoritmo detallado:

# 2Eva_IIT2019_T2 EDO, problema de valor inicial
# método iterativo
import numpy as np

# INGRESO
# longitud en x
a = 1
b = 2
# longitud en y
c = 1
d = 2
# tamaño de paso
dx = 0.25
dy = 0.25
# funciones en los bordes de la placa
abajo     = lambda x,y: x*np.log(x)
arriba    = lambda x,y: x*np.log(4*(x**2))
izquierda = lambda x,y: y*np.log(y)
derecha   = lambda x,y: 2*y*np.log(2*y)
# función de la ecuación
fxy = lambda x,y: x/y + y/x

# control de iteraciones
maxitera = 100
tolera = 0.0001

# PROCEDIMIENTO
# tamaño de la matriz
n = int((b-a)/dx)+1
m = int((d-c)/dy)+1
# vectores con valore de ejes
xi = np.linspace(a,b,n)
yj = np.linspace(c,d,m)
# matriz de puntos muestra
u = np.zeros(shape=(n,m),dtype=float)

# valores en los bordes
u[:,0]   = abajo(xi,yj[0])
u[:,m-1] = arriba(xi,yj[m-1])
u[0,:]   = izquierda(xi[0],yj)
u[n-1,:] = derecha(xi[n-1],yj)

# valores interiores la mitad en intervalo x,
# mitad en intervalo y, para menos iteraciones
mx = int(n/2)
my = int(m/2)
promedio = (u[mx,0]+u[mx,m-1]+u[0,my]+u[n-1,my])/4
u[1:n-1,1:m-1] = promedio

# método iterativo
itera = 0
converge = 0
while not(itera>=maxitera or converge==1):
    itera = itera +1
    # copia u para calcular errores entre iteraciones
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            # usar fórmula desarrollada para algoritmo
            fij = (dy**2)*fxy(xi[i],yj[j])
            u[i,j]=(u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-fij)/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
    if (erroru<tolera):
        converge=1

# SALIDA
print('iteraciones: ',itera)
print('error entre iteraciones: ',erroru)
print('solución para u: ')
print(u)

# Gráfica
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

# matrices de ejes para la gráfica 3D
X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()

grafica = Axes3D(figura)
grafica.plot_surface(X, Y, U, rstride=1, cstride=1,
                     cmap=cm.Reds)

plt.title('EDP elíptica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Ejercicio: 2Eva_IIT2019_T2 EDO, problema de valor inicial

la ecuación del problema planteado es:

y'(t) = f(t,y) = \frac{y}{2t^3}

0 ≤ t ≤ 1
y(0.5) = 1.5

---

literal a

La solución empieza usando la Serie de Taylor por ejemplo para tres términos:

y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h E = \frac{h^3}{3!} y'''(z) = O(h^3)

Se observa que se tiene el valor inicial y la primera derivada, si usamos tres términos se puede usar la segunda derivada.

y'(t) = f(t,y) = \frac{y}{2t^3} y''(t) = f'(t,y) = \frac{y'}{2t^3} + y \Big(\frac{-3}{2t^4}\Big) y''(t) = \frac{1}{2t^3}\Big( \frac{y}{2t^3} \Big) - \frac{3y}{2t^4} y''(t) = \frac{y}{4t^6} -\frac{3y}{2t^4}

por lo que la ecuación de Taylor a usar queda de la siguiente forma:

y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i y_{i+1} = y_{i} + h\frac{y}{2t^3}+\frac{h^2}{2} \Big( \frac{y}{4t^6} -\frac{3y}{2t^4}\Big)

que es la ecuacion que se usará con un error de O(h³)

Reemplazando los valores en la fórmula se obtiene la siguiente tabla:

estimado
 [ti,      yi,      d1yi,    d2yi]
[[  0.5    1.5      6.      -12.        ]
 [  0.6    2.04     4.7222  -12.68004115]
 [  0.7    2.4488   3.5697  -10.09510219]
 [  0.8    2.7553   2.6907  -7.46259891]
 [  0.9    2.9870   2.0487  -5.42399041]
 [  1.     3.1648   0.       0.        ]]

cuya gráfica es:

---

literal b

Para desarrollar Runge-Kutta de 2do orden se dispone de los siguientes datos:

y'(t) = f(t,y) = \frac{y}{2t^3}

t₀ = 0.5, y₀ = 1.5, h = 0.1

pasos del algoritmo,

K_1 = h * y'(t_i) K_2 = h * y'(t_i+h, y_i + K_1) y_{i+1} = y_i + \frac{K_1+K_2}{2} t_i = t_i + h</p> <p>

iteración 1: i = 0

K_1 = 0.1 * y'(0.5) = 0.6 K_2 = 0.1 * y'(0.5+0.1, 1.5 + 0.6) = 0.4861 y_{1} = 1.5 + \frac{0.6+0.4861}{2} = 2.0430 t_1 = 0.5 + 0.1 = 0.6</p> <p>

iteración 2: i = 1

K_1 = 0.1 * y'(0.6) = 0.4729 K_2 = 0.1 * y'(0.6+0.1, 2.0430 + 0.4729) = 0.3667 y_{1} = 2.0430 + \frac{0.4729+0.3667}{2} = 2.4629 t_1 = 0.6 + 0.1 = 0.7</p> <p>

iteración 3: i = 2

K_1 = 0.1 * y'(0.7) = 0.3590 K_2 = 0.1 * y'(0.7+0.1, 2.4629 + 0.3590) = 0.3667 y_{1} = 2.4629 + \frac{0.3590+0.3667}{2} = 2.7802 t_1 = 0.7 + 0.1 = 0.8</p> <p>

obteniendo la siguiente tabla:

estimado
 [ti,   yi,     K1,     K2]
[[0.5   1.5     0.      0.    ]
 [0.6   2.0430  0.6     0.4861]
 [0.7   2.4629  0.4729  0.3667]
 [0.8   2.7802  0.3590  0.2755]
 [0.9   3.0206  0.2715  0.2093]
 [1.    3.2048  0.2071  0.1613]]
Diferencias entre Taylor y Runge-Kutta2
[ 0.   -0.0030 -0.0140 -0.0248 -0.0335 -0.0400]

---

Algoritmo en Python

# EDO. Método de Taylor 3 términos
# Runge-Kutta de 2 Orden
# 2Eva_IIT2019_T2 EDO, problema de valor inicial
import numpy as np

# Funciones desarrolladas en clase
def edo_taylor3t(d1y,d2y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        estimado[i-1,2:]= [d1y(x,y),d2y(x,y)]
        y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i,0:2] = [x,y]
    return(estimado)
def rungekutta2(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)
        yi = yi + (K1+K2)/2
        xi = xi + h
        estimado[i] = [xi,yi,K1,K2]
    return(estimado)

# PROGRAMA PRUEBA
# INGRESO
# d1y = y', d2y = y''
d1y = lambda t,y: y/(2*t**3)
d2y = lambda t,y: y/(4*t**6)-(3/2)*y/(t**4)
t0 = 0.5
y0 = 1.5
h = 0.1
muestras = 5

# PROCEDIMIENTO
# Edo con Taylor
puntos = edo_taylor3t(d1y,d2y,t0,y0,h,muestras)
ti = puntos[:,0]
yi = puntos[:,1]

# Runge-Kutta
puntosRK2 = rungekutta2(d1y,t0,y0,h,muestras)
# ti = puntosRK2[:,0] # lo mismo del anterior
yiRK2 = puntosRK2[:,1]

# diferencias
diferencia = yi-yiRK2

# SALIDA
print('estimado[ti, yi, d1yi, d2yi]')
print(puntos)

print('estimado[ti, yi, K1, K2]')
print(puntosRK2)

print('Diferencias entre Taylor y Runge-Kutta2')
print(diferencia)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(ti[0],yi[0],'o', color='r',
         label ='[t0,y0]')
plt.plot(ti[1:],yi[1:],'o', color='g',
         label ='y con Taylor 3 términos')
plt.plot(ti[1:],yiRK2[1:],'o', color='m',
         label ='y Runge-Kutta 2 Orden')
plt.title('EDO: Solución numérica')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show())

Ejercicio: 2Eva_IIT2019_T1 Canteras y urbanizaciones

Literal a. área de cantera

Canteras– frontera superior

xi	55	85	195	305	390	780	1170
f(xi)	752	825	886	1130	1086	1391	1219

Para proceder se calculan los tamaños de paso, h, en cada intervalo:

dx – tamaños de paso

dxi

30

110

110

85

390

390

___

Ics = \frac{30}{2}(752+825) + \frac{110}{2}(825+886) + + \frac{110}{2}(1130+886) + \frac{85}{2}(1130+1086) + \frac{390}{2}(1086+1391) +\frac{390}{2}(1391+1219)

que tiene como resultado: Ics = 1342435.0

Canteras– frontera inferior

xi	55	…	705	705	850	850	1010	1170
f(xi)	260	…	260	550	741	855	855	1055

Para proceder se calculan los tamaños de paso, h, en cada intervalo:

dx – tamaños de paso

dxi

650

…

0.0

145

0.0

160

160

____

Se observa que existen rectángulos en los intervalos, por lo que se simplifica la fórmula.

Ici = (650)(260) + \frac{145}{2}(741+550) + + (160)(855) + \frac{160}{2}(1055+855)

cuyo resultado es: Ici =552197.5

El área correspondiente a la cantera es:

I_cantera = Ics -Ici =1342435.0 – 552197.5 = 790237.5

---

Literal b. área de urbanización
la frontera inferior está referenciada a la eje x con g(x)=0, por lo que solo es necesario realizar el integral para la frontera superior. El valor de la integral de la frontera inferior de la urbanización es cero.

Urbanización – frontera superior

xi	720	800	890	890	1170	1220
g(xi)	527	630	630	760	760	533

dx – tamaños de paso

dxi

80

90

0.0

280

50

____

Ius = \frac{80}{2}(527+630) + (90)(630) + + (280)(760) + \frac{50}{2}(760+533)

El valor del área de la urbanización es:

I_u = I_us – I_ui = 348105.0 – 0 = 348105.0

---

literal c

Se pude mejora la precisión para los intervalos donde el tamaño de paso es igual, sin necesidad de aumentar o quitar puntos.

Observando los tramaños de paso en cada sección se sugiere usar el método de Simpson de 1/3 donde existen dos tamaños de paso iguales y de forma consecutiva.

Cantera – frontera superior: en el intervalo xi= [85,195,305] donde h es= 110

Cantera – frontera inferior: en el intervalo xi = [850,110,1170] donde h es= 160

---

Algoritmo con Python

Para trapecios en todos los intervalos. Considera que si es un rectángulo, la fórmula del trapecio también funciona.

# 2Eva_IIT2019_T1 Canteras y urbanizaciones
import numpy as np
import matplotlib.pyplot as plt

# Funciones para integrar realizadas en clase
def itrapecio (xi,fi):
    n=len(fi)
    integral=0
    for i in range(0,n-1,1):
        h = xi[i+1]-xi[i]
        darea = (h/2)*(fi[i]+fi[i+1])
        integral = integral + darea 
    return(integral)

# INGRESO
# Canteras - frontera superior
xcs = [  55.,  85, 195,  305,  390,  780, 1170]
ycs = [ 752., 825, 886, 1130, 1086, 1391, 1219]
# Canteras - frontera inferior
xci = [ 55., 705, 705, 850, 850, 1010, 1170]
yci = [260., 260, 550, 741, 855,  855, 1055]

# Urbanización - frontera superior
xus = [720., 800, 890, 890, 1170, 1220]
yus = [527., 630, 630, 760,  760,  533]
# Urbanización - frontera inferior
xui = [720., 1220]
yui = [  0.,    0]

# PROCEDIMIENTO

# Area de cantera
Ics = itrapecio(xcs,ycs)
Ici = itrapecio(xci,yci)
Icantera = Ics-Ici

# Area de urbanización
Iurb = itrapecio(xus,yus)

# SALIDA
print('Area canteras: ',Icantera)
print('Area urbanización: ', Iurb)

# Gráfica canteras
plt.plot(xcs,ycs,color='brown')
plt.plot(xci,yci,color='brown')
plt.plot([xci[0],xcs[0]],[yci[0],ycs[0]],color='brown')
plt.plot([xci[-1],xcs[-1]],[yci[-1],ycs[-1]],color='brown')

# Gráfica urbanizaciones
plt.plot(xus,yus, color='green')
plt.plot(xui,yui, color='green')
plt.plot([xui[0],xus[0]],[yui[0],yus[0]], color='green')
plt.plot([xui[-1],xus[-1]],[yui[-1],yus[-1]], color='green')
plt.show()