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

Ejercicio: 1Eva_IIT2019_T1 Ecuación Recursiva

la ecuación recursiva es:

x_n = g(x_{n-1}) = \sqrt{3 + x_{n-1}}

literal a y b

g(x) es creciente en todo el intervalo, con valor minimo en g(1) = 2, y máximo en g(3) =2.449. Por observación de la gráfica, la pendiente g(x) es menor que la recta identidad en todo el intervalo

Verifique la cota de g'(x)

g(x) = \sqrt{3 + x} =(3+x)^{1/2} g'(x) =\frac{1}{2}(3+x)^{-1/2} g'(x) =\frac{1}{2\sqrt{3+x}}

Tarea: verificar que g'(x) es menor que 1 en todo el intervalo.

Literal c

Usando el algoritmo del punto fijo, iniciando con el punto x₀=2
y tolerancia de 10^-4, se tiene que:

Iteración 1: x₀=2

g(x_0) = \sqrt{3 + 2} = 2.2361

error = |2.2361 - 2| = 0.2361

Iteración 2: x₁ = 2.2361

g(x_1) = \sqrt{3 + 2.2361} = 2.2882

error = |2.2882 - 2.2361| = 0.0522

Iteración 3: x₂ = 2.2882

g(x_2) = \sqrt{3 + 2.2882} = 2.2996

error = |2.2996 - 2.28821| = 0.0114

Iteración 4: x₃ = 2.2996

g(x_3) = \sqrt{3 + 2.2996} = 2.3021

error = |2.3021- 2.2996| = 0.0025

Iteración 5: x₄ = 2.3021

g(x_4) = \sqrt{3 + 2.3021} = 2.3026

error = |2.3021- 2.2996| = 5.3672e-04

con lo que determina que el error en la 5ta iteración es de 5.3672e-04 y el error se reduce en casi 1/4 entre iteraciones. El punto fijo converge a 2.3028

Se muestra como referencia la tabla resumen.

[[ x ,   g(x), 	 tramo  ] 
 [1.      2.      1.    ]
 [2.      2.2361  0.2361]
 [2.2361  2.2882  0.0522]
 [2.2882  2.2996  0.0114]
 [2.2996  2.3021  0.0025]
 [2.3021  2.3026  5.3672e-04]
 [2.3026  2.3027  1.1654e-04]
 [2.3027  2.3028  2.5305e-05]
raiz:  2.3027686193257098

con el siguiente comportamiento de la funcion:

literal e

Realizando el mismo ejercicio para el método de la bisección, se requiere cambiar a la forma f(x)=0

x = \sqrt{3 + x} 0 = \sqrt{3 + x} -x f(x) = \sqrt{3 + x} -x

tomando como intervalo el mismo que el inicio del problema [1,3], al realizar las operaciones se tiene que:

a = 1 ; f(a) = 1
b = 3 ; f(b) = -0.551
c = (a+b)/2 = (1+3)/2 = 2
f(c) = f(2) = (3 + 2)^(.5) +2 = 0.236
Siendo f(c) positivo, y tamaño de paso 2, se reduce a 1

a = 2 ; f(a) = 0.236
b = 3 ; f(b) = -0.551
c = (a+b)/2 = (2+3)/2 = 2.5
f(c) = f(2.5) = (3 + 2.5)^(.5) +2.5 = -0.155
Siendo fc(c) negativo y tamaño de paso 1, se reduce a .5

a = 2
b = 2.5
...

Siguiendo las operaciones se obtiene la siguiente tabla:

[ i, a,   c,   b,    f(a),  f(c),  f(b), paso]
 1 1.000 2.000 3.000 1.000  0.236 -0.551 2.000 
 2 2.000 2.500 3.000 0.236 -0.155 -0.551 1.000 
 3 2.000 2.250 2.500 0.236  0.041 -0.155 0.500 
 4 2.250 2.375 2.500 0.041 -0.057 -0.155 0.250 
 5 2.250 2.312 2.375 0.041 -0.008 -0.057 0.125 
 6 2.250 2.281 2.312 0.041  0.017 -0.008 0.062 
 7 2.281 2.297 2.312 0.017  0.005 -0.008 0.031 
 8 2.297 2.305 2.312 0.005 -0.001 -0.008 0.016 
 9 2.297 2.301 2.305 0.005  0.002 -0.001 0.008 
10 2.301 2.303 2.305 0.002  0.000 -0.001 0.004 
11 2.303 2.304 2.305 0.000 -0.001 -0.001 0.002 
12 2.303 2.303 2.304 0.000 -0.000 -0.001 0.001 
13 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
14 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
15 2.303 2.303 2.303 0.000 -0.000 -0.000 0.000 
16 2.303 2.303 2.303 0.000  0.000 -0.000 0.000 
raiz:  2.302764892578125

Donde se observa que para la misma tolerancia de 10^-4, se incrementan las iteraciones a 16. Mientra que con punto fijo eran solo 8.

Nota: En la evaluación solo se requeria calcular hasta la 5ta iteración. Lo presentado es para fines didácticos

Ejercicio: 1Eva_IIT2019_T3 Circuito eléctrico

Las ecuaciones del problema son:

55 I_1 - 25 I_4 =-200 -37 I_3 - 4 I_4 =-250 -25 I_1 - 4 I_3 +29 I_4 =100 I_2 =-10

Planteo del problema en la forma A.X=B

A = [[ 55.0, 0,  0, -25],
     [  0  , 0,-37,  -4],
     [-25  , 0, -4,  29],
     [  0  ,  1, 0,   0]]

B = [-200,-250,100,-10]

El ejercicio se puede simplificar con una matriz de 3x3 dado que una de las corrientes I₂ es conocida con valor -10, queda resolver el problema para
[I₁ ,I₃ ,I₄ ]

A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

conformando la matriz aumentada

[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]

que se pivotea por filas para acercar a matriz diagonal dominante:

[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]

Literal a

Para métodos directos se aplica el método de eliminación hacia adelante.

Usando el primer elemento en la diagonal se convierten en ceros los números debajo de la posición primera de la diagonal

[[  55.    0.  -25.         -200.      ]
 [   0.  -37.   -4.         -250.      ]
 [   0.   -4.   17.636363      9.090909]]

luego se continúa con la segunda columna:

[[  55.    0.  -25.         -200.      ]
 [   0.  -37.   -4.         -250.      ]
 [   0.    0.   18.068796     36.117936]]

y para el método de Gauss se emplea sustitución hacia atrás
se determina el valor de I₄

18.068796 I_4 = 36.11793612 I_4 =\frac{36.11793612}{18.068796}= 1.99891216 -37 I_3 -4 I_4 = -250 -37 I_3= -250 + 4 I_4 I_3=\frac{-250 + 4 I_4}{-37} I_3=\frac{-250 + 4 (1.99891216)}{-37} = 6.54065815

y planteando se obtiene el último valor

55 I_1 +25 I_4 = -200 55 I_1 = -200 -25 I_4 I_1 = \frac{-200 -25 I_4}{55} I_1 = \frac{-200 -25(1.99891216)}{55} = -2.7277672

con lo que el vector solución es:

[-2.7277672   6.54065815  1.99891216]

sin embargo, para verificar la respuesta se aplica A.X=B

verificar que A.X = B, obteniendo nuevamente el vector B.
[-200.  -250.  100.]]

---

literal b

La norma de la matriz infinito se determina como:

||x|| = max\Big[ |x_i| \Big]

considere que en el problema el término en A de magnitud mayor es 55.
El vector suma de filas es:

[[| 55|+|  0|+|-25|],    [[80],
 [|  0|+|-37|+| -4|],  =  [41],
 [[-25|+| -4|+| 29|]]     [58]]

por lo que la norma ∞ ejemplo ||A||∞ 
es el maximo de suma de filas: 80

para revisar la estabilidad de la solución, se observa el número de condición

>>> np.linalg.cond(A)
4.997509004325602

En éste caso no está muy alejado de 1. De resultar alejado del valor ideal de uno,  la solución se considera poco estable. Pequeños cambios en la entrada del sistema generan grandes cambios en la salida.

Tarea: Matriz de transición de Jacobi

---

Literal c

En el método de Gauss-Seidel acorde a lo indicado, se inicia con el vector cero. Como no se indica el valor de tolerancia para el error, se considera tolera = 0.0001

las ecuaciones para el método son:

I_1 =\frac{-200 + 25 I_4}{55} I_3 = \frac{-250+ 4 I_4}{-37} I_4 =\frac{100 +25 I_1 + 4 I_3}{29}

Como I₂ es constante, no se usa en las iteraciones

I_2 =-10

teniendo como resultados de las iteraciones:

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X]
      diferencia,errado
0 [0. 0. 0.] 2e-05
1 [-3.6363636  6.7567568  1.2454461]
   [3.6363636 6.7567568 1.2454461] 6.756756756756757
2 [-3.0702518  6.6221139  1.7149021]
   [0.5661119 0.1346428 0.469456 ] 0.5661118513783094
3 [-2.8568627  6.5713619  1.891858 ]
   [0.2133891 0.050752  0.1769559] 0.2133891067340583
4 [-2.7764282  6.5522316  1.9585594]
   [0.0804345 0.0191304 0.0667014] 0.08043447732439457
5 [-2.7461094  6.5450206  1.9837016]
   [0.0303188 0.007211  0.0251423] 0.030318816370094925
6 [-2.7346811  6.5423025  1.9931787]
   [0.0114283 0.0027181 0.0094771] 0.011428316023939011
7 [-2.7303733  6.541278   1.996751 ]
   [0.0043078 0.0010246 0.0035723] 0.004307767346479974
8 [-2.7287495  6.5408918  1.9980975]
   [0.0016238 0.0003862 0.0013465] 0.001623761494915943
9 [-2.7281375  6.5407462  1.9986051]
   [0.0006121 0.0001456 0.0005076] 0.0006120575185017962
10 [-2.7279068  6.5406913  1.9987964]
   [2.3070778e-04 5.4871039e-05 1.9131760e-04] 0.00023070777766820427
11 [-2.7278198  6.5406707  1.9988685]
   [8.6962544e-05 2.0682983e-05 7.2114885e-05] 8.696254366213907e-05
12 [-2.727787   6.5406629  1.9988957]
   [3.2779493e-05 7.7962038e-06 2.7182845e-05] 3.277949307367578e-05
13 [-2.7277747  6.5406599  1.998906 ]
   [1.2355839e-05 2.9386860e-06 1.0246249e-05] 1.235583874370505e-05
14 [-2.72777    6.5406588  1.9989098]
   [4.6573860e-06 1.1077026e-06 3.8622013e-06] 4.6573859666665385e-06
numero de condición: 4.997509004325604
respuesta con Gauss-Seidel
[-2.72777    6.5406588  1.9989098]
>>>

con lo que el vector resultante es:

respuesta con Gauss-Seidel
[-2.72777 6.5406588 1.9989098]

que para verificar, se realiza la operación A.X
observando que el resultado es igual a B

[[-200.00002751]
 [-249.9999956 ]
 [ 100.0000125 ]]

---

---

Solución alterna

---

Usando la matriz de 4x4, los resultados son iguales para las corrientes
[I₁ ,I₂ , I₃ ,I₄ ]. Realizando la matriz aumentada,

[[  55.    0.    0.  -25. -200.]
 [   0.    0.  -37.   -4. -250.]
 [ -25.    0.   -4.   29.  100.]
 [   0.    1.    0.    0.  -10.]]

que se pivotea por filas para acercar a matriz diagonal dominante:

[[  55.    0.    0.  -25. -200.]
 [   0.    1.    0.    0.  -10.]
 [   0.    0.  -37.   -4. -250.]
 [ -25.    0.   -4.   29.  100.]]

Literal a

Para métodos directos se aplica el método de eliminación hacia adelante.

Usando el primer elemento  en la diagonal.

[[  55.     0.     0.   -25.         -200.        ]
 [   0.     1.     0.     0.          -10.        ]
 [   0.     0.   -37.    -4.         -250.        ]
 [   0.     0.    -4.    17.63636364    9.09090909]]

para el segundo no es necesario, por debajo se encuentran valores cero.
Por lo que se pasa al tercer elemento de la diagonal

[[  55.     0.     0.     -25.         -200.        ]
 [   0.     1.     0.      0.          -10.        ]
 [   0.     0.   -37.     -4.         -250.        ]
 [   0.     0.     0.     18.06879607    36.11793612]]

y para el método de Gauss se emplea sustitución hacia atras.
para x₄:

18.06879607 x_4 = 36.11793612 x_4 = 1.99891216

para x₃:

-37 x_3 -4 x_3 = -250 37 x_3 = 250-4 x_4 = 250-4(1.99891216) x_3 = 6.54065815

como ejercicio, continuar con x₁, dado que x₂=-10

55 x_1 + 25 x_4 = -200

El vector solución obtenido es:

el vector solución X es:
[[ -2.7277672 ]
 [-10.        ]
 [  6.54065815]
 [  1.99891216]]

sin embargo, para verificar la respuesta se aplica A.X=B.

[[-200.]
 [-250.]
 [ 100.]
 [ -10.]]

Se revisa el número de condición de la matriz:

>>> np.linalg.cond(A)
70.21827416891405

Y para éste caso, el número de condición se encuentra alejado del valor 1, contrario a la respuesta del la primera forma de solución con la matriz 3x3. De resultar alejado del valor ideal de uno, la solución se considera poco estable. Pequeños cambios en la entrada del sistema generan grandes cambios en la salida.

---

Algoritmo en Python

Presentado por partes para revisión:

Para el método de Gauss, los resultados del algoritmo se muestran como:

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Elimina hacia adelante:
 fila 0 pivote:  55.0
   factor:  0.0  para fila:  1
   factor:  -0.45454545454545453  para fila:  2
 fila 1 pivote:  -37.0
   factor:  0.10810810810810811  para fila:  2
 fila 2 pivote:  18.06879606879607
[[  55.            0.          -25.         -200.        ]
 [   0.          -37.           -4.         -250.        ]
 [   0.            0.           18.06879607   36.11793612]]
solución: 
[-2.7277672   6.54065815  1.99891216]

Instrucciones en Python usando las funciones creadas en la unidad:

# 1Eva_IIT2019_T3 Circuito eléctrico
# Método de Gauss
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B
import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

            pivoteado = pivoteado + 1
            if vertabla==True:
                print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print(AB)
    return(AB)

def gauss_eliminaAdelante(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss elimina hacia adelante
    tarea: verificar términos cero
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    if vertabla==True:
        print('Elimina hacia adelante:')
    for i in range(0,n,1):
        pivote = AB[i,i]
        adelante = i+1
        if vertabla==True:
            print(' fila',i,'pivote: ', pivote)
        for k in range(adelante,n,1):
            if (np.abs(pivote)>=casicero):
                factor = AB[k,i]/pivote
                AB[k,:] = AB[k,:] - factor*AB[i,:]
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
    if vertabla==True:
        print(AB)
    return(AB)

def gauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss sustituye hacia atras
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    # Sustitución hacia atras
    X = np.zeros(n,dtype=float) 
    ultfila = n-1
    ultcolumna = m-1
    for i in range(ultfila,0-1,-1):
        suma = 0
        for j in range(i+1,ultcolumna,1):
            suma = suma + AB[i,j]*X[j]
        X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
    return(X)

# INGRESO
A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_sustituyeAtras(AB,vertabla=True)

# SALIDA
print('solución: ')
print(X)

---

literal c

Resultados con el algoritmo de Gauss Seidel

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X]
      diferencia,errado
0 [0. 0. 0.] 2e-05
1 [-3.6363636  6.7567568  1.2454461]
   [3.6363636 6.7567568 1.2454461] 6.756756756756757
2 [-3.0702518  6.6221139  1.7149021]
   [0.5661119 0.1346428 0.469456 ] 0.5661118513783094
3 [-2.8568627  6.5713619  1.891858 ]
   [0.2133891 0.050752  0.1769559] 0.2133891067340583
4 [-2.7764282  6.5522316  1.9585594]
   [0.0804345 0.0191304 0.0667014] 0.08043447732439457
5 [-2.7461094  6.5450206  1.9837016]
   [0.0303188 0.007211  0.0251423] 0.030318816370094925
6 [-2.7346811  6.5423025  1.9931787]
   [0.0114283 0.0027181 0.0094771] 0.011428316023939011
7 [-2.7303733  6.541278   1.996751 ]
   [0.0043078 0.0010246 0.0035723] 0.004307767346479974
8 [-2.7287495  6.5408918  1.9980975]
   [0.0016238 0.0003862 0.0013465] 0.001623761494915943
9 [-2.7281375  6.5407462  1.9986051]
   [0.0006121 0.0001456 0.0005076] 0.0006120575185017962
10 [-2.7279068  6.5406913  1.9987964]
   [2.3070778e-04 5.4871039e-05 1.9131760e-04] 0.00023070777766820427
11 [-2.7278198  6.5406707  1.9988685]
   [8.6962544e-05 2.0682983e-05 7.2114885e-05] 8.696254366213907e-05
12 [-2.727787   6.5406629  1.9988957]
   [3.2779493e-05 7.7962038e-06 2.7182845e-05] 3.277949307367578e-05
13 [-2.7277747  6.5406599  1.998906 ]
   [1.2355839e-05 2.9386860e-06 1.0246249e-05] 1.235583874370505e-05
14 [-2.72777    6.5406588  1.9989098]
   [4.6573860e-06 1.1077026e-06 3.8622013e-06] 4.6573859666665385e-06
numero de condición: 4.997509004325604
respuesta con Gauss-Seidel
[-2.72777    6.5406588  1.9989098]
>>> 

Instrucciones en Python

# 1Eva_IIT2019_T3 Circuito eléctrico
# Algoritmo Gauss-Seidel
# solución de matrices
# métodos iterativos
# Referencia: Chapra 11.2, p.310,
#      Rodriguez 5.2 p.162
import numpy as np

def gauss_seidel(A,B,X0,tolera, iteramax=100,
                 vertabla=False, precision=4):
    ''' Método de Gauss Seidel, tolerancia, vector inicial X0
        para mostrar iteraciones: vertabla=True
    '''
    A = np.array(A, dtype=float)
    B = np.array(B, dtype=float)
    X0 = np.array(X0, dtype=float)
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = 2*tolera*np.ones(n, dtype=float)
    errado = np.max(diferencia)
    X = np.copy(X0)

    itera = 0
    if vertabla==True:
        print('Iteraciones Gauss-Seidel')
        print('itera,[X]')
        print('      diferencia,errado')
        print(itera, X, errado)
        np.set_printoptions(precision)
    while (errado>tolera and itera<iteramax):
        for i in range(0,n,1):
            xi = B[i]
            for j in range(0,m,1):
                if (i!=j):
                    xi = xi-A[i,j]*X[j]
            xi = xi/A[i,i]
            diferencia[i] = np.abs(xi-X[i])
            X[i] = xi
        errado = np.max(diferencia)
        itera = itera + 1
        if vertabla==True:
            print(itera, X)
            print('  ', diferencia, errado)
    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    return(X)

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas, entrega matriz aumentada AB
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    
    if vertabla==True:
        print('Matriz aumentada')
        print(AB)
        print('Pivoteo parcial:')
    
    # Pivoteo por filas AB
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    # Para cada fila en AB
    pivoteado = 0
    for i in range(0,n-1,1):
        # columna desde diagonal i en adelante
        columna = np.abs(AB[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(AB[i,:])
            AB[i,:] = AB[dondemax+i,:]
            AB[dondemax+i,:] = temporal

            pivoteado = pivoteado + 1
            if vertabla==True:
                print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
    if vertabla==True:
        if pivoteado==0:
            print('  Pivoteo por filas NO requerido')
        else:
            print('AB')
    return(AB)

# INGRESO
A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

X0  = [0.,0.,0.]

tolera = 0.00001
iteramax = 100
verdecimal = 7

# PROCEDIMIENTO
# numero de condicion
ncond = np.linalg.cond(A)

AB = pivoteafila(A,B,vertabla=True)
# separa matriz aumentada en A y B
[n,m] = np.shape(AB)
A = AB[:,:m-1]
B = AB[:,m-1]

respuesta = gauss_seidel(A,B,X0, tolera,
                         vertabla=True, precision=verdecimal)

# SALIDA
print('numero de condición:', ncond)
print('respuesta con Gauss-Seidel')
print(respuesta)

.

Ejercicio: 1Eva_IIT2019_T2 Proceso Termodinámico

la ecuación para el problema se describe como:

f(x)=e^{-0.5x}

ecuación que se usa para describir los siguientes puntos:

x	0	1	2	3	4
f(x)	1	0.60653065	0.36787944	0.22313016	0.13533528

Como el polinomio es de grado 2, se utilizan tres puntos. Para cubrir el intervalo los puntos seleccionados incluyen los extremos y el punto medio.

literal a

Con los puntos seleccionados se escriben las ecuaciones del polinomio:

p_2(x)= a_0 x^2 + a_1 x + a_2

usando los valores de la tabla:

p_2(0)=a_0 (0)^2 + a_1 (0) + a_2 = 1 p_2(2)=a_0 (2)^2 + a_1 (2) + a_2 = 0.36787944 p_2(4)=a_0 (4)^2 + a_1 (4) + a_2 = 0.13533528

con la que se escribe la matriz Vandermonde con la forma A.x=B

A= [[ 0.,  0.,  1.,]
    [ 4.,  2.,  1.,]
    [16.,  4.,  1.,]]

B= [[1.        ],
    [0.36787944],
    [0.13533528]]) 

matriz aumentada

[[ 0.,  0.,  1.,  1.        ]
 [ 4.,  2.,  1.,  0.36787944]
 [16.,  4.,  1.,  0.13533528]]

matriz pivoteada

[[16.,  4.,  1.,  0.13533528]
 [ 4.,  2.,  1.,  0.36787944]
 [ 0.,  0.,  1.,  1.        ]]

Resolviendo por algún método directo, la solución proporciona los coeficientes del polinomio

Tarea: escribir la solución del método directo, semejante a la presentada en el tema 3

[ 0.04994705 -0.41595438  1.        ]

con lo que el polinomio de interpolación es:

p_2(x) = 0.04994705 x^2 - 0.41595438 x + 1.0

en el enunciado se requiere la evaluación en x=2.4

p_2(2.4) = 0.04994705 (2.4)^2 - 0.41595438 (2.4) + 1.0 f(2.4)=e^{-0.5(2.4)} error = |f(2.4)-p_2(2.4)|

Evaluando en X1:  2.4
Evaluando p(x1):  0.2894044975129779
Error en x1:      0.011789714399224216
 Error relativo:  0.039143230291095066

La diferencia entre la función y el polinomio de interpolación se puede observar en la gráfica:

---

literal b

Tarea: Encontrar la cota de error con f(1.7)

---

Algoritmo en Python

Resultado con el algoritmo

Matriz Vandermonde: 
[[ 0.  0.  1.]
 [ 4.  2.  1.]
 [16.  4.  1.]]
los coeficientes del polinomio: 
[ 0.04994705 -0.41595438  1.        ]
Polinomio de interpolación: 
0.049947050111716*x**2 - 0.415954379637711*x + 1.0

 formato pprint
                   2                            
0.049947050111716*x  - 0.415954379637711*x + 1.0

Evaluando en X1:  2.4
Evaluando p(x1):  0.2894044975129779
Error en x1:      0.011789714399224216
 Error relativo:  0.039143230291095066

Evaluando en X2:  1.7
Evaluando p(x2):  0.2894044975129779
Error en x2:      0.011789714399224216
 Error relativo:  0.039143230291095066

Presentado por secciones, semejante a lo desarrollado en clases

# 1Eva_IIT2019_T2 Proceso Termodinámico
# El polinomio de interpolación
import numpy as np
import sympy as sym

# INGRESO
fx = lambda x: np.exp(-0.5*x)
xi =np.array([0,2,4],dtype=float)

# determina vector
fi= fx(xi)

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xi,dtype=float)
fi = np.array(fi,dtype=float)

B = fi
n = len(xi)

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype=float)
for i in range(0,n,1):
    for j in range(0,n,1):
        potencia = (n-1)-j # Derecha a izquierda
        D[i,j] = xi[i]**potencia

# Aplicar métodos Unidad03. Tarea
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)

# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i in range(0,n,1):
    potencia = (n-1)-i   # Derecha a izquierda
    termino = coeficiente[i]*(x**potencia)
    polinomio = polinomio + termino

# Polinomio a forma Lambda x:
# para evaluación con vectores de datos xin
muestras = 21
px = sym.lambdify(x,polinomio)

# SALIDA
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)


# literal b
x1 = 2.4
px1 = px(x1)
fx1 = fx(x1)
errorx1 = np.abs(px1-fx1)
errorx1rel = errorx1/fx1
x2 = 1.7
px2 = px(x1)
fx2 = fx(x1)
errorx2 = np.abs(px1-fx1)
errorx2rel = errorx1/fx1
print()
print('Evaluando en X1: ',x1)
print('Evaluando p(x1): ',px1)
print('Error en x1:     ',errorx1)
print(' Error relativo: ', errorx1rel)
print()
print('Evaluando en X2: ',x2)
print('Evaluando p(x2): ',px2)
print('Error en x2:     ',errorx2)
print(' Error relativo: ', errorx2rel)


# GRAFICA
import matplotlib.pyplot as plt
a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)

# Usando evaluación simbólica
##yin = np.zeros(muestras,dtype=float)
##for j in range(0,muestras,1):
##    yin[j] = polinomio.subs(x,xin[j])

plt.plot(xi,fi,'o', label='[xi,fi]')
plt.plot(xin,yin, label='p(x)')
plt.plot(xin,fx(xin), label='f(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.title(polinomio)
plt.show()

 

Ejercicio: 1Eva_IIT2019_T4 Concentración de químico

formula a usar:

C = C_{ent} ( 1 - e^{-0.04t})+C_{0} e^{-0.03t}

Se sustituyen los valores dados con:
C₀ = 4, C_ent = 10, C = 0.93 C_ent.

0.93(10) = 10 ( 1 - e^{-0.04t}) + 4 e^{-0.03t}

igualando a cero para forma estandarizada del algoritmo,

10( 1 - e^{-0.04t}) + 4 e^{-0.03t} - 9.3 = 0 0.7 - 10 e^{-0.04t} + 4 e^{-0.03t} = 0

se usas las funciones f(t) y f'(t) para el método de Newton-Raphson,

f(t) = 0.7 - 10 e^{-0.04t} + 4 e^{-0.03t} f'(t) = - 10(-0.04) e^{-0.04t} + 4(-0.03) e^{-0.03t} f'(t) = 0.4 e^{-0.04t} - 0.12 e^{-0.03t}

con lo que se pueden realizar los cálculos de forma iterativa.

t_{i+1} = t_i -\frac{f(t_i)}{f'(t_i)}

De no disponer de la gráfica de f(t), y se desconoce el valor inicial para t₀ se usa 0. Como no se indica la tolerancia, se estima en 10^-4

Iteración 1

t₀ = 0

f(0) = 0.7 - 10 e^{-0.04(0)} + 4 e^{-0.03(0)} = 5.3 f'(0) = 0.4 e^{-0.04(0)} - 0.12 e^{-0.03(0)} = -0.28 t_{1} = 0 -\frac{5.3}{-0.28} = 18.92 error = |18.92-0| = 18.92

Iteración 2

f(18.92) = 0.7 - 10 e^{-0.04(18.92)} + 4 e^{-0.03(18.92)} = -1.72308 f'(18.92) = 0.4 e^{-0.04(18.92)} - 0.12 e^{-0.03(18.92)} = 0.119593 t_{2} = 18.92 -\frac{-1.723087}{0.119593} = 33.3365 error = |33.3365 - 18.92| = 14.4079

Iteración 3

f(33.3365) = 0.7 - 10 e^{-0.04(33.3365)} + 4 e^{-0.03(33.3365)} = -0.4642 f'(33.3365) = 0.4 e^{-0.04(33.3365)} - 0.12 e^{-0.03(33.3365)} = 0.06128 t_{3} = 33.3365 -\frac{-0.46427}{-5.8013} = 40.912 error = |40.912 - 33.3365| = 7.5755

Observando que los errores disminuyen entre cada iteración, se encuentra que el método converge.

y se forma la siguiente tabla:

['xi' ,  'xnuevo', 'error']
[ 0.      18.9286  18.9286]
[18.9286  33.3365  14.4079]
[33.3365  40.912    7.5755]
[40.912   42.654     1.742]
[42.654   42.7316   0.0776]
[4.2732e+01 4.2732e+01 1.4632e-04]
raiz en:  42.731721341402796

Observando la gráfica de la función puede observar el resultado:

---

Algoritmo en Python

# 1Eva_IIT2019_T4
# Método de Newton-Raphson

import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda t: 0.7-10*np.exp(-0.04*t)+4*np.exp(-0.03*t)
dfx = lambda t:0.40*np.exp(-0.04*t)-0.12*np.exp(-0.03*t)

x0 = 0
tolera = 0.001

a = 0
b = 60
muestras = 21

# 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

tabla = np.array(tabla)
n=len(tabla)

# para la gráfica
xk = np.linspace(a,b,muestras)
fk = fx(xk)

# SALIDA
print(['xi', 'xnuevo', 'error'])
np.set_printoptions(precision = 4)
for i in range(0,n,1):
    print(tabla[i])
print('raiz en: ', xi)

# grafica
plt.plot(xk,fk)
plt.axhline(0, color='black')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.show()

Ejercicio: 3Eva_IT2019_T3 EDP Difusión en sólidos

Siguiendo el procedimiento planteado en la sección EDP parabólicas, se plantea la malla del ejercicio:

Para plantear la ecuación en forma discreta:

\frac{\phi_{i,j+1}-\phi_{i,j}}{\Delta t}=D\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{(\Delta x)^2}

y resolver usando el método explícito para ecuaciones parabólicas, obteniendo el siguiente resultado:

\phi_{i,j+1}-\phi_{i,j}=D\frac{\Delta t }{\Delta x^2}(\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}) \lambda = D\frac{\Delta t }{\Delta x^2} \phi_{i,j+1}-\phi_{i,j}=\lambda (\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}) \phi_{i,j+1} =\lambda \phi_{i+1,j}-2\lambda\phi_{i,j}+\lambda\phi_{i-1,j}+\phi_{i,j} \phi_{i,j+1} =\lambda \phi_{i+1,j}(1-2\lambda)\phi_{i,j}+\lambda\phi_{i-1,j} \phi_{i,j+1} =P \phi_{i+1,j}+Q\phi_{i,j}+R\phi_{i-1,j}

siendo:
P = λ = 0.16 (Δx/100)/Δx² = 0.0016/Δx = 0.0016/0.02=0.08
Q = 1-2λ = 1-2*(0.08) = 0.84
R = λ =0.08

\phi_{i,j+1} =0.08 \phi_{i+1,j}+ 0.84\phi_{i,j}+0.08\phi_{i-1,j}

Iteración 1 en tiempo:
i=1, j=0

\phi_{1,1} =0.08 \phi_{2,0}+ 0.84\phi_{1,0}+0.08\phi_{0,0} \phi_{1,1} =0.08 (0)+ 0.84(0)+0.08(5)=0.4

i=2,j=0

\phi_{2,1} =0.08 \phi_{3,0}+ 0.84\phi_{2,0}+0.08\phi_{1,0} = 0

Para los proximos valores i>2, todos los resultados son 0

Iteración 2 en tiempo
i=1, j=1

\phi_{1,2} =0.08 \phi_{2,0}+ 0.84\phi_{1,0}+0.08\phi_{0,0}

\phi_{1,2} =0.08 (0)+ 0.84(0.4)+0.08(5)=0.736
i=2, j=1

\phi_{2,2} =0.08 \phi_{3,1}+ 0.84\phi_{2,1}+0.08\phi_{1,1} \phi_{2,2} =0.08(0)+ 0.84(0)+0.08(0.4) = 0.032

i=3, j=1

\phi_{3,2} =0.08\phi_{4,1}+ 0.84\phi_{3,1}+0.08\phi_{2,1}=0

Para los proximos valores i>3, todos los resultados son 0

Tarea: Desarrollar la iteración 3 en el tiempo.

siguiendo las iteraciones se tiene la siguiente tabla:

[[5.0, 0.000, 0.000, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 0.400, 0.000, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 0.736, 0.032, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 1.021, 0.085, 0.00256, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 1.264, 0.153, 0.00901, 0.00020, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
...
]

Con lo que se obtiene la siguiente gráfica.

El resultado se interpreta mejor con una animación: (tarea)

Tarea: Presentar el orden de error de la ecuación basado en las fórmulas de diferenciación

---

Algorirmo en Python

# 3Eva_IT2019_T3 Difusión en sólidos
# EDP parabólicas. método explícito,usando diferencias finitas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 5
Tb = 0
T0 = 0
# longitud en x
a = 0
b = 0.1
# Constante K
K = 1/(1.6e-1)
# Tamaño de paso
dx = 0.02
dt = dx/100
# iteraciones en tiempo
n = 50

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
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
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('phi[i,j]')
plt.title('Solución EDP parabólica')
plt.show()

La animación se complementa con lo mostrado en la sección de Unidades.

Ejercicio: 3Eva_IT2019_T2 Integral con interpolación

El ejercicio considera dos partes: interpolación e integración

a. Interpolación

Se requiere aproximar la función usando tres puntos. Para comprender la razón del método solicitado, se compara la función con dos interpolaciones:

a.1 Lagrange
a.2 Trazador cúbico sujeto

Observando la gráfica se aclara que en éste caso, una mejor aproximación se obtiene con el método  de trazador cúbico sujeto. Motivo por lo que el tema tiene un peso de 40/100 puntos

Los valores a considerar para la evaluación son:

puntos referencia xi,yi: 
[0.         0.78539816 1.57079633]
[ 0.          0.62426595 -0.97536797]
derivadas en los extremos:  
    3.141592653589793 
    0.6929852019184021
Polinomio de Lagrange
-1.80262534301178*x**2 + 2.21061873102778*x
Trazadores cúbicos sujetos
[0.         0.78539816]
-0.548171611756137*x**3 - 2.55744517923506*x**2 + 3.14159265358979*x

[0.78539816 1.57079633]
4.66299098804068*x**3 - 14.8359577843727*x**2 + 12.7851139029174*x - 2.52466795930204

------------------
Valores calculados para Trazadores cúbicos sujetos:
Matriz A: 
[[-0.26179939 -0.13089969  0.        ]
 [ 0.78539816  3.14159265  0.78539816]
 [ 0.          0.13089969  0.26179939]]
Vector B: 
[  2.34675256 -16.9893436    2.72970237]
coeficientes S: 
[-5.11489036 -7.69808822 14.27573913]
coeficientes a,b,c,d
[-0.54817161  4.66299099]
[-2.55744518 -3.84904411]
[ 3.14159265 -1.89005227]
[0.         0.62426595]

b. Integración

Como forma de comparacíon de resultados, se requiere integrar con varios métodos para comparar resultados y errores.

b.1 Integración con Cuadratura de Gauss, usando el resultado de trazador cúbico.

Se integra en cada tramo de cada polinomio:

Trazadores cúbicos sujetos
[0.         0.78539816]
-0.548171611756137*x**3 - 2.55744517923506*x**2 + 3.14159265358979*x

Se obtienen los puntos del método de cuadratura desplazados en el rango:

xa:  0.16597416116944688
xb:  0.6194240022280014
area:  0.5037962958529855

Para el segundo tramo:

[0.78539816 1.57079633]
4.66299098804068*x**3 - 14.8359577843727*x**2 + 12.7851139029174*x - 2.52466795930204
xa:  0.9513723245668951
xb:  1.4048221656254496
area:  -0.2706563884589365

Con lo que el integral total es:

Integral total:  0.23313990739404894

b.2 Integración analítica

\int_0^{\pi /2}sin(\pi x) dx

u = πx
du/dx = π
dx = du/π

se convierte en:

\frac{1}{\pi}\int sin(u) du \frac{1}{\pi}(-cos(u))

volviendo a la variable x:

\frac{1}{\pi}(-cos(\pi x)) \Big\rvert_{0}^{\frac{\pi}{2}} -\frac{1}{\pi}(cos(\pi \frac{\pi}{2})-cos(\pi(0))) = 0.24809580527879377

c. Estimación del error

Se restan los resultados de las secciones b.1 y b.2

error = |0.24809580527879377 - 0.23313990739404894 |

error = 0.014955897884744829

---

Algoritmo en Python

separado por literales

# 3Eva I T 2019 Interpola e Integra
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if (j!=i):
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

def traza3sujeto(xi,yi,u,v):
    '''
    Trazador cúbico sujeto, splines
    resultado: polinomio en forma simbólica
    '''
    n = len(xi)
    # Valores h
    h = np.zeros(n-1, dtype=float)
    # Sistema de ecuaciones
    A = np.zeros(shape=(n,n), dtype=float)
    B = np.zeros(n, dtype=float)
    S = np.zeros(n-1, dtype=float)
    # coeficientes
    a = np.zeros(n-1, dtype=float)
    b = np.zeros(n-1, dtype=float)
    c = np.zeros(n-1, dtype=float)
    d = np.zeros(n-1, dtype=float)
    
    polinomios=[]
    
    if (n>=3):
        for i in range(0,n-1,1):
            h[i]=xi[i+1]-xi[i]
        A[0,0] = -h[0]/3
        A[0,1] = -h[0]/6
        B[0] = u-(yi[1]-yi[0])/h[0]
        for i in range(1,n-1,1):
            A[i,i-1] = h[i-1]
            A[i,i] = 2*(h[i-1]+h[i])
            A[i,i+1] = h[i]
            B[i] = 6*((yi[i+1]-yi[i])/h[i] - (yi[i]-yi[i-1])/h[i-1])
        A[n-1,n-2] = h[n-2]/6
        A[n-1,n-1] = h[n-2]/3
        B[n-1] = v-(yi[n-1]-yi[n-2])/h[n-2]

        # Resolver sistema de ecuaciones
        S = np.linalg.solve(A,B)

        # Coeficientes
        for i in range(0,n-1,1):
            a[i]=(S[i+1]-S[i])/(6*h[i])
            b[i]=S[i]/2
            c[i]=(yi[i+1]-yi[i])/h[i]-(2*h[i]*S[i]+h[i]*S[i+1])/6
            d[i]=yi[i]
      
        # polinomio en forma simbólica
        x=sym.Symbol('x')
        polinomios=[]
        for i in range(0,n-1,1):
            ptramo = a[i]*(x-xi[i])**3 + b[i]*(x-xi[i])**2 + c[i]*(x-xi[i])+ d[i]
            ptramo = ptramo.expand()
            polinomios.append(ptramo)
        parametros = [A,B,S,a,b,c,d]                                                           
    return(polinomios, parametros)

# INGRESO
f = lambda x: np.sin(np.pi*x)
muestrasf = 20
a = 0
b = np.pi/2
# Derivadas en los extremos
u = np.pi*np.cos(np.pi*a)
v = np.pi*np.cos(np.pi*b)
muestras = 3

# literal a
# PROCEDIMIENTO
xif = np.linspace(a,b,muestrasf)
yif = f(xif)

xi = np.linspace(a,b,muestras)
yi = f(xi)

# Usando Lagrange
x = sym.Symbol('x')
pL = interpola_lagrange(xi,yi)
pxL = sym.lambdify(x,pL)
pxiL =  pxL(xif)

# Trazador cúbico sujeto
pS, parametros = traza3sujeto(xi,yi,u,v)
pxiS = np.zeros(muestrasf,dtype=float)

# Evalua trazadores cúbicos sujetos
i=0
ap = xi[i]
bp = xi[i+1]
poli = sym.lambdify(x, pS[i])
for j in range(0,muestrasf,1):
    punto = xif[j]
    if (punto>bp):
        i = i+1
        ap = xi[i]
        bp = xi[i+1]
        poli = sym.lambdify(x,pS[i])
    pxiS[j] = poli(punto)

# SALIDA
print('puntos referencia xi,yi: ')
print(xi)
print(yi)
print('derivadas en los extremos: ',u,v)
print('Polinomio de Lagrange')
print(pL)
print('Trazadores cúbicos sujetos')
n = len(xi)
for i in range(0,n-1,1):
    print(xi[i:i+2])
    print(pS[i])
# Parametros de Trazadores cúbicos sujetos
print('Matriz A: ')
print(parametros[0])
print('Vector B: ')
print(parametros[1])
print('coeficientes S: ')
print(parametros[2])
print('coeficienetes a,b,c,d')
print(parametros[3])
print(parametros[4])
print(parametros[5])
print(parametros[6])

# Gráficas
plt.plot(xif,yif, label='funcion')
plt.plot(xi,yi,'o', label='muestras')
plt.plot(xif,pxiL, label='p(x)_Lagrange')
plt.plot(xif,pxiS, label='p(x)_Traza3Sujeto')
plt.legend()
plt.xlabel('x')
plt.show()

# literal b
# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    print('xa: ',xa)
    print('xb: ',xb)
    print('area: ', area)
    return(area)

# INGRESO
f0 = sym.lambdify(x,pS[0])
f1 = sym.lambdify(x,pS[1])
# Procedimiento
I0 = integraCuadGauss2p(f0,xi[0],xi[1])
I1 = integraCuadGauss2p(f1,xi[1],xi[2])
It = I0+I1

# SALIDA
print('Integral 1: ', I0)
print('Integral 2: ', I1)
print('Integral total: ',It)