s3Eva_IIT2009_T2 Valor inicial Runge-Kutta 4to orden dy/dx

Ejercicio: 3Eva_IIT2009_T2 Valor inicial Runge-Kutta 4to orden dy/dx

La ecuación del problema:

(1-x^2)y' - xy = x (1-x^2)

se despeja:

(1-x^2)y' = x (1-x^2) - xy y' = x - \frac{xy}{(1-x^2)}

con valores iniciales x0 = 0 , y0 = 2, h = 0.1 en el intervalo [0, 0.5]

Usando Runge-Kutta de 2do Orden

Iteración 1

K_1 = h y'(0,2) = (0.1)\Big[0 - \frac{0 (2)}{(1-(0^2))}\Big] = 0 K_2 = h y'(0+0.1, 2 + 0) = (0.1)\Big[0.1 - \frac{0.1 (2+0)}{(1-(0.1^2))}\Big] = 0.0302 y_1= 2 + \frac{0+0.0302}{2} = 2.0151 x_1 = x_0 + h = 0 + 0.1 = 0.1

Iteración 2

K_1 = h y'(0.1,2.0151) = (0.1)\Big[0.1 - \frac{0.1 (2.0151)}{(1-(0.1^2))}\Big] = 0.0304 K_2 = h y'(0.1+0.1, 2.0151 + 0.0304) = (0.1)\Big[0.2 - \frac{0.2 (2.0151 + 0.0304)}{(1-(0.2^2))}\Big] = 0.0626 y_2 = 2.0151 + \frac{0.0304+0.0626}{2} = 2.0616 x_2 = x_1 + h = 0.1 + 0.1 = 0.2

Iteración 3

K_1 = h y'(0.2,2.0616) = (0.1)\Big[0.2 - \frac{0.2 (2.0616)}{(1-(0.2^2))}\Big] = 0.0629 K_2 = h y'(0.2+0.1, 2.0616 + 0.0629) = (0.1)\Big[0.3 - \frac{0.3 (2.0616 + 0.0629)}{(1-(0.3^2))}\Big] = 0.1 y_3 = 2.0151 + \frac{0.0629+0.1}{2} = 2.1431 x_3 = x_2 + h = 0.2 + 0.1 = 0.3

siguiendo el algoritmo se completa la tabla:

 [xi,    yi,    K1,    K2    ]
[[0.     2.     0.     0.    ]
 [0.1    2.0151 0.     0.0302]
 [0.2    2.0616 0.0304 0.0626]
 [0.3    2.1431 0.0629 0.1   ]
 [0.4    2.2668 0.1007 0.1468]
 [0.5    2.4463 0.1479 0.211 ]]

Algoritmo en Python

# 3Eva_IIT2009_T2 Valor inicial Runge-Kutta 4to orden dy/dx
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
d1y = lambda x,y : x + x*y/(1-x**2)
x0 = 0
y0 = 2
h = 0.1
a = 0
b = 1/2

# PROCEDIMIENTO
muestras = int((b -a)/h)+1
tabla = np.zeros(shape=(muestras,4),dtype=float)

i = 0
xi = x0
yi = y0
tabla[i,:] = [xi,yi,0,0]

i = i+1
while not(i>=muestras):
    K1 = h* d1y(xi,yi)
    K2 = h* d1y(xi+h,yi+K1)
    yi = yi + (K1+K2)/2
    xi = xi +h
    tabla[i,:] = [xi,yi,K1,K2]
    i = i+1
# vector para gráfica
xg = tabla[:,0]
yg = tabla[:,1]

# SALIDA
# muestra 4 decimales
np.set_printoptions(precision=4)
print(' [xi, yi, K1, K2]')
print(tabla)

# Gráfica
plt.plot(xg,yg)
plt.xlabel('xi')
plt.ylabel('yi')
plt.grid()
plt.show()

Tarea: Realizar con Runge-Kutta de 4to Orden

s3Eva_IT2009_T2 EDO Taylor Seno(x)

Ejercicio: 3Eva_IT2009_T2 EDO Taylor Seno(x)

La ecuación del problema es:

xy'+ 2y = \sin (x) \frac{\pi}{2} \leq x \leq \frac{3\pi}{2} y\Big(\frac{\pi}{2} \Big) = 1

Para el algoritmo se escribe separando y’:

y' = \frac{\sin (x)}{x} - 2\frac{y}{x}

tomando como referencia Taylor de 3 términos más el término de error O(h3)

y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i + \frac{h^3}{3!}y'''_i

Se usa hasta el tercer término para el algoritmo.

y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i

Se determina que se requiere la segunda derivada para completar la aproximación. A partir de la ecuación del problema se aplica en cada término:

\Big(\frac{u}{v}\Big)' = \frac{u'v-uv' }{v^2} y'' = \frac{\cos (x) x - sin(x)}{x^2} - 2\frac{y'x-y}{x^2} y'' = \frac{\cos (x) x - sin(x)}{x^2} - 2\frac{\Big(\frac{\sin (x)}{x} - 2\frac{y}{x} \Big)x-y}{x^2}

Con lo que se realizan las iteraciones para llenar la tabla

iteración 1

x0 = π/2= 1.57079633

y0= 1

y' = \frac{\sin (π/2)}{π/2} - 2\frac{1}{π/2} = -0.40793719 y'' = \frac{\cos (π/2) π/2 - sin(π/2)}{(π/2)^2}+ - 2\frac{(-0.40793719)(π/2)-1}{(π/2)^2}= 0.48531359 y_{1} = 1 +\frac{\pi/10}{1!}(-0.40793719) + \frac{(\pi/10)^2}{2!}(0.48531359) = 0.86

x1 = x0+h = 1.57079633 + π/10 =1.88495559

iteración 2

x1 = 1.88495559

y1 = 0.86

y' = \frac{\sin (1.88495559)}{1.88495559} - 2\frac{0.86}{1.88495559} = -0.31947719 y'' = \frac{\cos (1.88495559)(1.88495559) - sin(1.88495559)}{(1.88495559)^2} - 2\frac{(-0.31947719) (1.88495559)-y}{(1.88495559)^2} = 0.16854341 y_{2} = 0.86 +\frac{\pi /10}{1!}(-0.31947719) + \frac{(\pi /10)^2}{2!}(0.16854341) = 0.75579202

x2 = x1+ h = 1.88495559 + π/10 = 2.19911486

iteración 3

x2 = 2.19911486

y2 = 0.75579202

y' = \frac{\sin (2.19911486)}{2.19911486} - 2\frac{y}{2.19911486} = -0.29431724 y'' = \frac{\cos (2.19911486)(2.19911486) - sin(2.19911486)}{(2.19911486)^2} + - 2\frac{(-0.29431724)(2.19911486)-0.75579202}{(2.19911486)^2} = 0.0294177 y_{3} = 0.75579202 +\frac{\pi/10}{1!}(-0.29431724) + \frac{(\pi /10)^2}{2!}(0.0294177)

x3 = x3+h = 2.19911486 + π/10 = 2.19911486

Con lo que se puede realizar el algoritmo, obteniendo la siguiente tabla:

estimado
 [xi,          yi,         d1y,        d2y       ]
[[ 1.57079633  1.         -0.40793719  0.48531359]
 [ 1.88495559  0.86       -0.31947719  0.16854341]
 [ 2.19911486  0.75579202 -0.29431724  0.0294177 ]
 [ 2.51327412  0.66374258 -0.29583247 -0.02247944]
 [ 2.82743339  0.5727318  -0.30473968 -0.02730493]
 [ 3.14159265  0.47868397 -0.31027009 -0.00585871]
 [ 3.45575192  0.38159973 -0.30649476  0.02930236]
 [ 3.76991118  0.28383639 -0.29064275  0.06957348]
 [ 4.08407045  0.18899424 -0.26221809  0.10859762]
 [ 4.39822972  0.10111944 -0.22243506  0.14160656]
 [ 4.71238898  0.02410027  0.          0.        ]]

Algoritmo en Python

# 3Eva_IT2009_T2 EDO Taylor Seno(x)
# EDO. Método de Taylor 3 términos 
# estima solucion para muestras separadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

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):
        y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i-1,2:]= [d1y(x,y),d2y(x,y)]
        estimado[i,0:2] = [x,y]
    return(estimado)

# PROGRAMA PRUEBA
# 3Eva_IT2009_T2 EDO Taylor Seno(x).

# INGRESO.
# d1y = y', d2y = y''
d1y = lambda x,y: np.sin(x)/x - 2*y/x
d2y = lambda x,y: (x*np.cos(x)-np.sin(x))/(x**2)-2*(x*(np.sin(x)/x - 2*y/x)-y)/(x**2)
x0 = np.pi/2
y0 = 1
h = np.pi/10
muestras = 10

# PROCEDIMIENTO
puntos = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

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

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

s2Eva_IT2008_T1_MN Producción petroleo

Ejercicio: 2Eva_IT2008_T1_MN Producción petroleo

Literal a

Para el cálculo de las derivadas se hace uso de las fórmulas de diferenciación presentadas en la unidad 6, y basadas en el polinomio de Taylor:

f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} + O(h) f''(x_i) = \frac{f(x_{i+2})-2f(x_{i+1})+f(x_i)}{h^2} + O(h)
[ dia, prod, dprod, d2prod]
[[ 1.000e+00  3.345e+03 -1.000e+02  6.600e+01]
 [ 2.000e+00  3.245e+03 -3.400e+01  1.320e+02]
 [ 3.000e+00  3.211e+03  9.800e+01 -5.600e+01]
 [ 4.000e+00  3.309e+03  4.200e+01  1.900e+01]
 [ 5.000e+00  3.351e+03  6.100e+01 -2.430e+02]
 [ 6.000e+00  3.412e+03 -1.820e+02  8.700e+01]
 [ 7.000e+00  3.230e+03 -9.500e+01  9.200e+01]
 [ 8.000e+00  3.135e+03 -3.000e+00  0.000e+00]
 [ 9.000e+00  3.132e+03 -3.000e+00  0.000e+00]
 [ 1.000e+01  3.129e+03  0.000e+00  0.000e+00]]

representados en las siguientes gráfica:

literal b

Dado que las fórmlas de error usadas tienen error del orden h: O(h), el error de las fórmulas es del orden de:

h= dia[1]-dia[0] = 2-1 = 1

literal c

Para el día dos se observa un decrecimiento en la producción, tal como lo refleja el valor negativo de la primera derivada.
Sin embargo para el día siguiente, la producción no mantiene la tasa de decrecimiento, se observa la segunda derivada positiva, Empieza a «acelerar».


Las instrucciones en Python para la tabla presentada son:

# 2Eva_IT2008_T1_MN Producción petroleo
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
dia = np.array( [1., 2, 3, 4, 5, 6, 7, 8, 9, 10])
produccion = [3345., 3245, 3211, 3309, 3351, 3412, 3230, 3135, 3132, 3129]
produccion = np.array(produccion, dtype = float)

# PROCEDIMIENTO
n = len(dia)

# primera derivada
dp = np.zeros(n,dtype=float)
for i in range(0,n-1,1):
    dp[i] = (produccion[i+1]-produccion[i])/(dia[i+1]-dia[i])

# segunda derivada
d2p = np.zeros(n,dtype=float)
h = dia[1]-dia[0]
for i in range(0,n-2,1):
    d2p[i] = (produccion[i]-2*produccion[i+1]+ produccion[i+2])/(h**2)

tabla = np.concatenate(([dia],[produccion],[dp],[d2p]),axis=0)
tabla = np.transpose(tabla)

# SALIDA
print(" [ dia, prod, dprod, d2prod]")
print(tabla)

# gráfica
plt.subplot(121)
plt.plot(dia,produccion)
plt.xlabel('dia')
plt.ylabel('producción')
plt.grid()
plt.subplot(122)
plt.plot(dia,dp,color='green',label='dp')
plt.xlabel('dia')
plt.plot(dia,d2p,color='orange',label='d2p')
plt.axhline(0)
plt.legend()
plt.grid()
plt.show()

s2Eva_IIT2007_T2_AN Lanzamiento vertical proyectil

Ejercicio: 2Eva_IIT2007_T2_AN Lanzamiento vertical proyectil

la ecuación del problema:

m \frac{\delta v}{\delta t} = -mg - kv|v|

se despeja:

\frac{\delta v}{\delta t} = -g - \frac{k}{m}v|v|

y usando los valores indicados en el enunciado:

\frac{\delta v}{\delta t} = -9,8 - \frac{0.002}{0.11}v|v|

con valores iniciales de:

t0 = 0 , v0 = 8 , h=0.2

Como muestra inicial, se usa Runge-Kutta de 2do Orden

iteración 1

K_1 = h\frac{\delta v}{\delta t}(0, 8) = 0.2[-9,8 - \frac{0.002}{0.11}8|8|] = -2.1927 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 8 -2.1927 ) = 0.2[-9,8 - \frac{0.002}{0.11}(8 -2.1927)|8 -2.1927|] =-2.0826 v_1 = -9,8 +\frac{-2.1927-2.0826 }{2} = 5.8623 t_1 = t_0 + h = 0 + 0.2 = 0.2 error = O(h^3) = O(0.2^3) = O(0.008)

iteración 2

K_1 = h\frac{\delta v}{\delta t}(0.2, 5.8623) = 0.2[-9,8 - \frac{0.002}{0.11}(5.8623)|5.8623|] = -2.085 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 5.8623 -2.085) = 0.2[-9,8 - \frac{0.002}{0.11}(5.8623 -2.085)|5.8623 -2.085|] =-2.0119 v_2 = -9,8 +\frac{-2.085-2.0119}{2} = 3.8139 t_2 = t_1 + h = 0.2 + 0.2 = 0.4

iteración 3

K_1 = h\frac{\delta v}{\delta t}(0.4, 3.8139) = 0.2[-9,8 - \frac{0.002}{0.11}( 3.8139)| 3.8139|] = -2.0129 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 3.8139 -2.0129) = 0.2[-9,8 - \frac{0.002}{0.11}(3.8139 -2.0129)|3.8139 -2.0129|] =-1.9718 v_3 = -9,8 +\frac{-2.0129-1.9718}{2} = 1.8215 t_3 = t_2 + h = 0.4 + 0.2 = 0.6

Tabla y gráfica del ejercicio para todo el intervalo:

 [xi,      yi,     K1,     K2    ]
[[ 0.      8.      0.      0.    ]
 [ 0.2     5.8623 -2.1927 -2.0826]
 [ 0.4     3.8139 -2.085  -2.0119]
 [ 0.6     1.8215 -2.0129 -1.9718]
 [ 0.8    -0.1444 -1.9721 -1.9599]
 [ 1.     -2.0964 -1.9599 -1.9439]]

Algoritmo en Python

# 3Eva_IT2009_T2 EDO Taylor Seno(x)
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
d1y = lambda t,v: -9.8-(0.002/0.11)*v*np.abs(v)
x0 = 0
y0 = 8
h = 0.2
a = 0
b = 1

# PROCEDIMIENTO
muestras = int((b -a)/h)+1
tabla = np.zeros(shape=(muestras,4),dtype=float)

i = 0
xi = x0
yi = y0
tabla[i,:] = [xi,yi,0,0]

i = i+1
while not(i>=muestras):
    K1 = h*d1y(xi,yi)
    K2 = h*d1y(xi+h,yi+K1)
    yi = yi + (K1+K2)/2
    xi = xi +h
    tabla[i,:] = [xi,yi,K1,K2]
    i = i+1
# vector para gráfica
xg = tabla[:,0]
yg = tabla[:,1]

# SALIDA
# muestra 4 decimales
np.set_printoptions(precision=4)
print(' [xi, yi, K1, K2]')
print(tabla)
# Gráfica
plt.plot(xg,yg)
plt.xlabel('ti')
plt.ylabel('yi')
plt.grid()
plt.show()

Tarea: Realizar iteraciones para Runge-Kutta de 4to Orden

s3Eva_IIT2007_T3 EDO Taylor orden 2

Ejercicio: 3Eva_IIT2007_T3 EDO Taylor orden 2

La ecuación del problema es:

y'= 1 +\frac{y}{t} + \Big(\frac{y}{t}\Big) ^2 1\leq t\leq 2 y(1)=0, h=0.2

Tomando como referencia Taylor de 3 términos más el término de error O(h3)

y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i + \frac{h^3}{3!}y'''_i

Se usa hasta el tercer término para el algoritmo.

y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i

Se determina que se requiere la segunda derivada para completar la aproximación. A partir de la ecuación del problema se aplica en cada término:

\Big(\frac{u}{v}\Big)' = \frac{u'v-uv' }{v^2} y'= 1 +\frac{y}{t} + \frac{y^2}{t^2} y''= 0+\frac{y't-y}{t^2} + \frac{2yy't^2-2y^2t}{t^4} y''= \frac{y't-y}{t^2} + \frac{2y\Big(1 +\frac{y}{t} + \frac{y^2}{t^2}\Big) t^2-2y^2t}{t^4}

Con lo que se realizan las iteraciones para llenar la tabla

iteración 1

t0 = 1

y0= 0

y'= 1 +\frac{0}{1} + \Big(\frac{0}{1}\Big)^2 = 1 y''= \frac{(1)(1)-0}{(1)^2} + \frac{2(0)(1)(1)^2-2(0)^2(1)}{(1)^4} = 1 y_{1} = 0 +\frac{0.2}{1!}(1) + \frac{0.2^2}{2!}(1) = 0.22

t1 = t0 + h = 1 + 0.2 = 1.2

iteración 2

t1 = 1.2

y1= 0.22

y'= 1 +\frac{0.22}{1.2} + \Big(\frac{0.22}{1.2}\Big) ^2 = 1.21694444 y''= \frac{y'(1.2)-0.22}{1.2^2} + \frac{2(0.22)y'(1.2)^2-2(0.22)^2(1.2)}{(1.2)^4} = 1.17716821 y_{2} = 0.22 +\frac{0.2}{1!}(1.21694444) + \frac{0.2^2}{2!}(1.17716821) = 0.48693225

t2 = t1 + h = 1.2 + 0.2 = 1.4

iteración 3

t2 = 1.4

y2= 0.48693225

y'= 1 +\frac{0.48693225}{1.4} + \Big(\frac{0.48693225}{1.4}\Big) ^2 = 1.46877968 y''= \frac{(1.46877968)(1.4)-0.48693225}{1.4^2} + \frac{2(0.48693225)(1.46877968)(1.4)^2-2(0.48693225)^2(1.4)}{1.4^4} = 1.35766995 y_{3} = 0.48693225 +\frac{0.2}{1!}(1.46877968) + \frac{0.2^2}{2!}(1.35766995) = 0.80784159

t3 = t2 + h = 1.4 + 0.2 = 1.6

Con lo que se puede realizar el algoritmo, obteniendo la siguiente tabla:

estimado
 [xi,        yi,        d1yi,      d2yi]
[[1.         0.         1.         1.        ]
 [1.2        0.22       1.21694444 1.17716821]
 [1.4        0.48693225 1.46877968 1.35766995]
 [1.6        0.80784159 1.759826   1.57634424]
 [1.8        1.19133367 2.09990017 1.8564435 ]
 [2.         1.64844258 0.         0.        ]]

s3Eva_IIT2007_T1 EDP Eliptica, problema de frontera

Ejercicio: 3Eva_IIT2007_T1 EDP Eliptica, problema de frontera

Con los datos del ejercicio se plantean de la siguiente forma en los bordes:

La ecuación a resolver es:

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

Que en su forma discreta, con diferencias divididas centradas es:

\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} = 4

Al agrupar constantes se convierte en:

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

Siendo Δx= 1/3 y Δy =2/3, se mantendrá la relación λ = (Δx/Δy)2

u[i-1,j]-2u[i,j]+u[i+1,j] + + \lambda u[i,j-1]-2\lambda u[i,j]+\lambda u[i,j+1] = = 4(\Delta x)^2
u[i-1,j]-2(1+\lambda)u[i,j]+u[i+1,j] + + \lambda u[i,j-1]+\lambda u[i,j+1] = 4(\Delta x)^2

Se pueden realizar las iteraciones para los nodos y reemplaza los valores de frontera:

i=1 y j=1

u[0,1]-2(1+\lambda)u[1,1]+u[2,1] + + \lambda u[1,0]+\lambda u[1,2] = 4(\Delta x)^2 \Delta y^2-2(1+\lambda)u[1,1]+u[2,1] + + \Delta x^2+\lambda u[1,2] = 4(\Delta x)^2 -2(1+\lambda)u[1,1]+u[2,1] +\lambda u[1,2] = 4(\Delta x)^2 - \Delta y^2 - \Delta x^2 -2(1+0.25)u[1,1]+u[2,1] +0.25 u[1,2] = 4(1/3)^2 - (2/3)^2- (1/3)^2 -2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111

i=2 , j=1

u[1,1]-2(1+\lambda)u[2,1]+u[3,1] + + \lambda u[2,0]+\lambda u[2,2] = 4(\Delta x)^2 u[1,1]-2(1+\lambda)u[2,1]+(\delta y-1)^2 + + \Delta x^2 +\lambda u[2,2] = 4(\Delta x)^2 u[1,1]-2(1+\lambda)u[2,1] + \lambda u[2,2] = 4(\Delta x)^2 - (\Delta y-1)^2 - \Delta x^2 u[1,1]-2(1+0.25)u[2,1] + 0.25 u[2,2] = 4(1/3)^2 - (2/3-1)^2 - (1/3)^2 u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222

i=1 , j=2

u[0,2]-2(1+\lambda)u[1,2]+u[2,2] + + \lambda u[1,1]+\lambda u[1,3] = 4(\Delta x)^2 (2\Delta y^2)-2(1+\lambda)u[1,2]+u[2,2] + + \lambda u[1,1]+\lambda (\Delta x-1)^2 = 4(\Delta x)^2 -2(1+\lambda)u[1,2]+u[2,2] + \lambda u[1,1] = 4(\Delta x)^2 - (2\Delta y^2) -\lambda (\Delta x-1)^2 -2(1+0.25)u[1,2]+u[2,2] + 0.25 u[1,1] = 4(1/3)^2 - (2(2/3)^2] -0.25 (1/3-1)^2 -2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555

i=2, j=2

u[1,2]-2(1+\lambda)u[2,2]+u[3,2] + + \lambda u[2,1]+\lambda u[2,3] = 4(\Delta x)^2 u[1,2]-2(1+\lambda)u[2,2]+(\Delta y-1)^2 + + \lambda u[2,1]+\lambda (\Delta x-1)^2 = 4(\Delta x)^2 u[1,2]-2(1+\lambda)u[2,2] + \lambda u[2,1] =4(\Delta x)^2- (\Delta y-1)^2 -\lambda (\Delta x-1)^2 u[1,2]-2(1+0.25)u[2,2] + 0.25 u[2,1] =4(1/3)^2- (2/3-1)^2 -0.25(1/3-1)^2 u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222

Obtiene el sistema de ecuaciones a resolver:

-2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111 u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222 -2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555 u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222

Se escribe en forma matricial Ax=B

\begin {bmatrix} -2.5 && 1&&0.25&&0 \\1&&-2.5&&0&&0.25\\0.25&&0&&-2.5&&1\\0&&0.25&&1&&-2.5\end{bmatrix} \begin {bmatrix} u[1,1] \\ u[2,1]\\ u[1,2]\\u[2,2] \end{bmatrix} = \begin {bmatrix}-0.1111\\0.2222\\-0.5555\\-0.2222 \end{bmatrix}

que se resuelve con un algoritmo:

import numpy as np

A = np.array([[-2.5,1,0.25,0],
              [1,-2.5,0,0.25],
              [0.25,0,-2.5,1],
              [0,0.25,1,-2.5]])
B = np.array([-0.1111,0.2222,-0.5555,-0.2222])

x = np.linalg.solve(A,B)

print(x)

y los resultados son:

[ 0.05762549 -0.04492835  0.31156835  0.20901451]

s1Eva_IT2017_T2 Tanque esférico-volumen

Ejercicio: 1Eva_IT2017_T2 Tanque esférico-volumen

a. Planteamiento del problema

V = \frac{\pi h^{2} (3r-h)}{3}

Si r=1 m y V=0.75 m3,

0.75 = \frac{\pi h^{2} (3(1)-h)}{3} 0.75 -\frac{\pi h^{2} (3(1)-h)}{3} = 0 f(h) = 0.75 -\frac{\pi h^{2} (3-h)}{3} = 0.75 -\frac{\pi}{3}(h^{2} (3)-h^{2}h) = 0.75 -\frac{\pi}{3}(3 h^{2} - h^{3}) f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3}

b. Intervalo de búsqueda de raíz

El tanque vacio tiene h=0 y completamente lleno h= 2r = 2(1) = 2, por lo que el intevalo tiene como extremos:

[0,2]

Verificando que exista cambio de signo en la función:

f(0) = 0.75 -\pi (0)^{2} + \frac{\pi}{3}(0)^{3} = 0.75 f(2) = 0.75 -\pi (2)^{2} + \frac{\pi}{3}(2)^{3}= -3.4387

y verificando al usar la gráfica dentro del intervalo:

Tolerancia

Se indica en el enunciado como 0.01 que es la medición mínima a observar con un flexómetro.

tolera = 0.01


c. Método de Newton-Raphson
d. Método de Punto Fijo


c. Método de Newton-Raphson

El método de Newton-Raphson requiere la derivada de la función:

x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)} f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3} f'(h) = -2\pi h + \pi h^{2}

Tomando como punto inicial de búsqueda el extremo izquierdo del intervalo, genera una división para cero. Por lo que se mueve un poco a la derecha, algo más cercano a la raiz, viendo la gráfica por ejemplo 0.1

x0 = 0.1

iteración 1

i =0

f(0.1) = 0.75 -\pi (0.1)^{2} + \frac{\pi}{3}(0.1)^{3} =0.7196 f'(0.1) = -2\pi (0.1) + \pi (0.1)^{2} = -0.5969 x_{1} = x_0 -\frac{0.7496 }{-0.0625} = 1.3056 tramo = |x_{1}-x_{0}| = |0.1-1.3056 | = 1.2056

iteración 2

i =1

f(1.3056) = 0.75 -\pi (1.3056)^{2} + \frac{\pi}{3}(1.3056)^{3} = -2.2746 f'(1.3056) = -2\pi (1.3056) + \pi (1.3056)^{2} =-2.8481 x_{1} = x_0 -\frac{-2.2746}{-2.8481} = 0.5069 tramo = |x_{2}-x_{1}|=|0.5069-1.3056|=0.7987

iteración 3

i =2

f(0.5069) = 0.75 -\pi (0.5069)^{2} + \frac{\pi}{3}(0.5069)^{3} = 0.0789 f'(0.5069) = -2\pi (0.5069) + \pi (0.5069)^{2} =-2.3780 x_{1} = x_0 -\frac{0.0789}{-2.3780} = 0.5401 tramo = |x_{3}-x_{2}| =|0.5401 - 0.5069| = 0.0332

Observe que el tramo disminuye en cada iteración , por lo que el método converge, si se siguen haciendo las operaciones se tiene que:

 [ xi, xnuevo, tramo]
[[1.00000000e-01 1.30560920e+00 1.20560920e+00]
 [1.30560920e+00 5.06991599e-01 7.98617601e-01]
 [5.06991599e-01 5.40192334e-01 3.32007350e-02]
 [5.40192334e-01 5.39518667e-01 6.73667593e-04]]
raiz 0.5395186666699257

Instrucciones en Python

para Método de Newton-Raphson

# 1Eva_IT2017_T2 Tanque esférico-volumen
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda h: 0.75 - np.pi*(h**2)+(np.pi/3)*h**3
dfx = lambda h: -2*np.pi*h+np.pi*(h**2)

# Parámetros de método
a = 0
b = 2
tolera = 0.01
x0 = 0.1
iteramax = 100

# parametros de gráfica
La = a
Lb = b
muestras = 21

# PROCEDIMIENTO
# Newton-Raphson
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)

# calcula para grafica
hi = np.linspace(La,Lb,muestras)
fi = fx(hi)
gi = dfx(hi)

# SALIDA
print(' [ xi, xnuevo, tramo]')
print(tabla)
print('raiz', xnuevo)
plt.plot(hi,fi)
plt.plot(xi,fx(xi),'ro')
plt.axhline(0, color='green')
plt.xlabel('h')
plt.ylabel('V')
plt.title('Método Newton-Raphson')
plt.show()

Planteo con Punto Fijo


d. Método de Punto Fijo

Del planteamiento del problema en el literal a, se tiene que:

0.75 = \frac{\pi h^{2} (3(1)-h)}{3}

de donde se despeja una h:

\frac{3(0.75)}{\pi (3(1)-h) } = h^{2} h = \sqrt{\frac{3*0.75}{\pi (3-h) }} h = \sqrt{\frac{2.25}{\pi (3-h) }}

con lo que se obtienen las expresiones a usar en el método

identidad = h g(h) = \sqrt{\frac{2.25}{\pi (3-h) }}

El punto inicial de búsqueda debe encontrarse en el intervalo, se toma el mismo valor que x0 en el método de Newton-Raphson

x0 = 0.10

Iteracion 1

x_0= 0.10 g(0.10) = \sqrt{\frac{2.25}{\pi (3-(0.10) }}= 0.4969 tramo= |0.4969-0.10| = 0.3869

Iteracion 2

x_1= 0.4969 g(0.4969) = \sqrt{\frac{2.25}{\pi (3-(0.4969 ) }}= 0.5349 tramo= |0.5349- 0.4969| = 0.038

Iteracion 3

x_2 =0.5349 g(0.5349) = \sqrt{\frac{2.25}{\pi (3-(0.5349) }}= 0.5390 tramo= |0.5390 - 0.5349| = 0.0041

con lo que se cumple el criterio de tolerancia, y se obtiene la raiz de:

raiz = 0.5390

Tabla de resultados, donde se observa que el tramo o error en cada iteración disminuye, por lo que el método converge.

 [i,xi,xi+1,tramo]
[[1.         0.1        0.4969553  0.3969553 ]
 [2.         0.4969553  0.5349116  0.03795631]
 [3.         0.5349116  0.53901404 0.00410243]]
raiz 0.5390140355891347
>>> 

Instrucciones en Python

para Método de Punto-Fijo, recordamos que el método puede diverger, por lo que se añade el parámetro iteramax

# 1Eva_IT2017_T2 Tanque esférico-volumen
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda h: h
gx = lambda h: np.sqrt(2.25/(np.pi*(3-h)))

a = 0.1
b = 2
tolera = 0.01
iteramax = 100

La = a
Lb = b
muestras = 21

# PROCEDIMIENTO
# Punto Fijo
tabla = []
i = 1 # iteración
b = gx(a)
tramo = abs(b-a)
while not(tramo<=tolera or i>=iteramax):
    tabla.append([i,a,b,tramo])
    a = b
    b = gx(a)
    tramo = abs(b-a)
    i = i+1
tabla.append([i,a,b,tramo])
respuesta = b

# Validar respuesta
if (i>=iteramax):
    respuesta = np.nan
tabla = np.array(tabla)

# calcula para grafica
hi = np.linspace(La,Lb,muestras)
fi = fx(hi)
gi = gx(hi)

# SALIDA
print(' [i, xi, xi+1, tramo]')
print(tabla)
print('raiz', respuesta)
plt.plot(hi,fi, label = 'identidad', color='green')
plt.plot(hi,gi, label = 'g(h)', color = 'orange')
plt.plot(b,gx(b),'ro')
plt.axhline(0, color='green')
plt.xlabel('h')
plt.ylabel('V')
plt.title('Método Punto Fijo')
plt.legend()
plt.axhline(0, color='green')
plt.show()

Otra forma de probar la convergencia es que |g'(x)|<1 que se observa en la una gráfica adicional, lo que limita aún más el intervalo de búsqueda.

Desarrollo en la siguiente clase.


s1Eva_IT2016_T3_MN Tasa interés anual

Ejercicio: 1Eva_IT2016_T3_MN Tasa interés anual

Propuesta de Solución  empieza con el planteamiento del problema, luego se desarrolla con el método de Bisección y método del Punto Fijo solo con el objetivo de comparar resultados.


Planteamiento del problema

La fórmula del enunciado para el problema es:

A = P \frac{i(1+i)^{n}}{(1+i)^{n} -1}

que con los datos dados se convierte a:

5800 = 35000 \frac{i(1+i)^8}{(1+i)^8 -1} 35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0

que es la forma de f(x) = 0

f(i)=35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800

Intervalo de búsqueda

Como el problema plantea la búsqueda de una tasa de interés, consideramos que:

  • Las tasas de interés no son negativas. ⌉(i<0)
  • Las tasas de interés no son cero en las instituciones bancarias (i≠0)
  • Las tasas muy grandes 1 = 100/100 = 100% tampoco tienen mucho sentido

permite acotar la búsqueda a un intervalo (0,1].
Sin embargo tasas demasiado altas tampoco se consideran en el problema pues el asunto es regulado (superintendencia de bancos), por lo que se podría intentar entre un 1% = 0.01l y la mitad del intervalo 50%= 0.5 quedando

[0.01,0.5]

Tolerancia

si la tolerancia es de tan solo menor a un orden de magnitud que el valor buscado, se tiene que las tasas de interés se representan por dos dígitos después del punto decimal, por lo que la tolerancia debe ser menor a eso.

Por ejemplo: tolerancia < 0.001 o aumentando la precisión

tolera = 0.0001


Método de la Bisección

itera = 1

a = 0.01, b = 0.5

c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255 f(0.01) = 35000 \frac{0.01(1+0.01)^8}{(1+0.01)^8-1}-5800 = -1225.83 f(0.5) = 35000 \frac{0.5(1+0.5i)^8}{(1+0.5)^8-1}-5800 = 12410.54

con lo que se verifica que existe cambio de signo al evaluar f(x) en el intervalo y puede existir una raíz.

f(0.255) = 35000 \frac{0.255(1+0.255)^8}{(1+0.255)^8-1}-5800 = 4856.70

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.255]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.255-0.01 = 0.245

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

itera = 2

a = 0.01, b = 0.255

c = \frac{0.01+0.255}{2} = 0.1325 f(0.01) = -1225.83 f(0.255) = 4856.70 f(0.1325 ) = 35000 \frac{0.1325(1+0.1325)^8}{(1+0.1325)^8-1}-5800 = 1556.06

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.1325]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.1325-0.01 = 0.1225

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

itera = 3

a = 0.01, b = 0.1325

c = \frac{0.01+0.1225}{2} = 0.071 f(0.01) = -1225.83 f(0.1325) = 1556.06 f(0.071 ) = 35000 \frac{0.071(1+0.071)^8}{(1+0.071)^8-1}-5800 = 89.79

lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.071]

El tramo permite estimar el error, reduce el intervalo a:

tramo = b-a = 0.071-0.01 = 0.061

valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

Para una evaluación del tema en forma escrita es suficiente para mostrar el objetivo de aprendizaje, el valor final se lo encuentra usando el algoritmo.

       raiz en:  0.06724243164062502
error en tramo:  5.981445312500111e-05
iteraciones:  13
>>>

Algoritmo en Python

# 1Eva_IT2016_T3_MN Tasa interés anual
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: 35000*(x*(1+x)**8)/((1+x)**8 -1) -5800
a = 0.01
b = 0.5
tolera = 0.0001

# PROCEDIMIENTO
cuenta = 0
np.set_printoptions(precision=3)
tramo = b-a
while not(tramo<tolera):
    c = (a+b)/2
    fa = fx(a)
    fb = fx(b)
    fc = fx(c)
    cambia = np.sign(fa)*np.sign(fc)
    if cambia < 0: a = a b = c if cambia > 0:
        a = c
        b = b
    tramo = b-a
    cuenta = cuenta+1

# SALIDA
print('       raiz en: ', c)
print('error en tramo: ', tramo)
print('iteraciones: ',cuenta)


Método del Punto Fijo

El planteamiendo del punto fijo se realiza con x= g(x), por lo que se reordena la ecuación a:

35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0 35000 \frac{i(1+i)^8}{(1+i)^8 -1}=5800 \frac{i(1+i)^8}{(1+i)^8 -1}=\frac{5800}{35000} i=\frac{58}{350} \frac{(1+i)^8 -1} {i(1+i)^8}

con lo que g(x) es:

g(i)=\frac{58}{350} \frac{(1+i)^8 -1} {(1+i)^8}

valor inicial de búsqueda

Para el punto inicial i0 se puede usar uno de los extremos del intervalo propuesto en la sección de planteamiento. Para reducir aun más la búsqueda se pude seleccionar el punto intermedio

 i0 = 0.255

Itera = 1

 i0 = 0.255

g(0.255i)=\frac{58}{350} \frac{(1+0.255)^8 -1} {(1+0.255)^8} = 0.1387

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.1387 - 0.255| = 0.1163

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Itera = 2

i1 = g(i0 ) = 0.1387

g(0.1387)=\frac{58}{350} \frac{(1+0.1387)^8 -1} {(1+0.1387)^8} = 0.1071

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.1071 - 0.1387| = 0,0316

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Itera = 3

i2 = g(i1 ) = 0.1071

g(0.1071)=\frac{58}{350} \frac{(1+0.1071)^8 -1} {(1+0.1071)^8} = 0.0922

el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

tramo = | nuevo-antes| = |0.0922 - 0.1071| = 0,0149

como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

Observación: Como el error disminuye entre cada iteración, se considera que el método converge, si se realizan suficientes iteraciones se cumplierá con el valor de tolerancia y se habrá llegado a la precisión requerida.

Usando el algoritmo se tiene:

iteraciones: 17
    raiz en: 0.06754389199556779
>>>

Algoritmo en Python

# Algoritmo de Punto Fijo
# x0 valor inicial de búsqueda
# error = tolera

import numpy as np
def puntofijo(gx,antes,tolera, iteramax=50):
    itera = 1 # iteración
    nuevo = gx(antes)
    tramo = abs(nuevo-antes)
    while(tramo>=tolera and itera<=iteramax ):
        antes = nuevo
        nuevo = gx(antes)
        tramo = abs(nuevo-antes)
        itera = itera+1
    respuesta = nuevo
    # Validar respuesta
    if (itera>=iteramax ):
        respuesta = np.nan
    print('iteraciones:',itera)
    return(respuesta)

# PROGRAMA ---------

# INGRESO
gx = lambda i: (58/350)*((1+i)**8-1)/((1+i)**8)

a = 0       # intervalo
b = 0.5

x0 = 0.255
tolera = 0.0001
iteramax  = 50      # itera máximo

# PROCEDIMIENTO
respuesta = puntofijo(gx,0.255,tolera,iteramax)

# SALIDA
print('    raiz en:',respuesta)

s2Eva_IT2015_T2 Deflexión de mástil

Ejercicio: 2Eva_IT2015_T2 Deflexión de mástil

\frac{\delta ^2 y}{\delta x^2} = \frac{f}{2EI} (L-x)^2 x= 0, y = 0, \frac{\delta y}{\delta x} = 0

Ecuación en forma estandarizada:

y'' = \frac{f}{2EI} (L-x)^2

Sistema de ecuaciones

y' = z = f(x,y,z) z' = \frac{f}{2EI} (L-x)^2 = g(x,y,z)

siendo:

\frac{f}{2EI} =\frac{60}{2(1.25 \times 10 ^{-8}) (0.05)} = 4.8 \times 10^{-6}

El mástil mide L=30 m y se requieren 30 intervalos, entonces h = 30/30 = 1.

Usando muestras = tramos +1 = 31

Se plantea una solución usando Runge Kutta de 2do Orden

f(x,y,z) = z g(x,y,z) = (4.8 \times 10^{-6})(30-x)^2

Desarrollo analítico

Las iteraciones se realizan para llenar los datos de la tabla,

tabla de resultados
i xi yi zi
0 0 0 0
1 1 0.00216 0.00417
2 2 0.00835 0.00807
3 3 tarea

iteración 1

i = 0 ; x0= 0; y0 = 0; z0 = 0

K_{1y} = h f(x_0,y_0, z_0)= 1(0) = 0

 

K_{1z} = h g(x_0,y_0, z_0) = = 1(4.8 \times 10^{-6})(30-0)^2 = 0.00432 K_{2y} = h f(x_0+h,y_0+K_{1y}, z_0+K_{1z}) = = 1(0+0.00432) = 0.00432 K_{2z} = h g(x_0+h,y_0+K_{1y}, z_0+K_{1z}) = = 1(4.8 \times 10^{-6})(30-(0+1))^2 = 0.00403 y_1 = y_0 + \frac{K_{1y}+K_{2y}}{2} = = 0 + \frac{0+0.00432}{2} = 0.00216 z_1 = z_0 + \frac{K_{1z}+K_{2z}}{2} = = 0 + \frac{0.00432+0.00403}{2} = 0.00417

iteración 2

i = 2 ; x1= 1; y1 = 0.00216; z1 = 0.00417

K_{1y} = h f(x_1,y_1, z_1)= 1(0.00417) = 0.00417 K_{1z} = h g(x_1,y_1, z_1) = = 1(4.8 \times 10^{-6})(30-1)^2 = 0.00403 K_{2y} = h f(x_1+h,y_1+K_{1y}, z_1+K_{1z}) = = 1(0.00417+0.00403) = 0.00821 K_{2z} = h g(x_1+h,y_1+K_{1y}, z_1+K_{1z}) = = 1(4.8 \times 10^{-6})(30-(1+1))^2 = 0.00376 y_2 = y_1 + \frac{K_{1y}+K_{2y}}{2} = = 0.00216 + \frac{0.00417+0.00821}{2} = 0.00835 z_2 = z_2 + \frac{K_{1z}+K_{2z}}{2} = = 0.00417 + \frac{0.00403+0.00376}{2} = 0.00807

iteración 3
i = 2; x2= 2; y2 = 0.00835; z2 = 0.00807
tarea

Para facilitar los cálculos se propone usa el algoritmo en Python, como se describe en la siguiente sección.


Algoritmo en Python

Al usar el algoritmo se puede comparar los resultados entre Runge-Kutta de 2do orden y de 4to Orden.

De los resultados se presenta la siguiente gráfica

Observe en la gráfica la diferencia de escalas entre los ejes.

Runge Kutta 2do Orden
 [x, y, z]
[[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  1.00000000e+00   2.16000000e-03   4.17840000e-03]
 [  2.00000000e+00   8.35680000e-03   8.07840000e-03]
 [  3.00000000e+00   1.83168000e-02   1.17096000e-02]
 [  4.00000000e+00   3.17760000e-02   1.50816000e-02]
 [  5.00000000e+00   4.84800000e-02   1.82040000e-02]
 ...
 [  2.90000000e+01   9.29856000e-01   4.32216000e-02]
 [  3.00000000e+01   9.73080000e-01   4.32240000e-02]
 [  3.10000000e+01   1.01630400e+00   4.32264000e-02]]
Runge Kutta 4do Orden
 [x, y, z]
[[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  1.00000000e+00   2.11240000e-03   4.17760000e-03]
 [  2.00000000e+00   8.26240000e-03   8.07680000e-03]
 [  3.00000000e+00   1.81764000e-02   1.17072000e-02]
 [  4.00000000e+00   3.15904000e-02   1.50784000e-02]
 [  5.00000000e+00   4.82500000e-02   1.82000000e-02]
 ...
 [  2.90000000e+01   9.28800400e-01   4.31984000e-02]
 [  3.00000000e+01   9.72000000e-01   4.32000000e-02]
 [  3.10000000e+01   1.01520040e+00   4.32016000e-02]]
>>> 

Las instrucciones en Python obtener los resultados:

# 2Eva_IT2015_T2 Deflexión de mástil
# solución propuesta: edelros@espol.edu.ec
import numpy as np

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

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

def rungekutta4fg(fx,gx,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,3),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * fx(xi,yi,zi)
        K1z = h * gx(xi,yi,zi)
        
        K2y = h * fx(xi+h/2, yi + K1y/2, zi + K1z/2)
        K2z = h * gx(xi+h/2, yi + K1y/2, zi + K1z/2)
        
        K3y = h * fx(xi+h/2, yi + K2y/2, zi + K2z/2)
        K3z = h * gx(xi+h/2, yi + K2y/2, zi + K2z/2)

        K4y = h * fx(xi+h, yi + K3y, zi + K3z)
        K4z = h * gx(xi+h, yi + K3y, zi + K3z)

        yi = yi + (K1y+2*K2y+2*K3y+K4y)/6
        zi = zi + (K1z+2*K2z+2*K3z+K4z)/6
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# INGRESO
f = 60
L = 30
E = 1.25e8
I = 0.05
x0 = 0
y0 = 0
z0 = 0
tramos = 30

fx = lambda x,y,z: z
gx = lambda x,y,z: (f/(2*E*I))*(L-x)**2

# PROCEDIMIENTO
muestras = tramos + 1
h = L/tramos
tabla2 = rungekutta2fg(fx,gx,x0,y0,z0,h,muestras)
xi2 = tabla2[:,0]
yi2 = tabla2[:,1]
zi2 = tabla2[:,2]

tabla4 = rungekutta4fg(fx,gx,x0,y0,z0,h,muestras)
xi4 = tabla4[:,0]
yi4 = tabla4[:,1]
zi4 = tabla4[:,2]

# SALIDA
print('Runge Kutta 2do Orden')
print(' [x, y, z]')
print(tabla2)
print('Runge Kutta 4do Orden')
print(' [x, y, z]')
print(tabla4)

# GRAFICA
import matplotlib.pyplot as plt
plt.plot(xi2,yi2, label='Runge-Kutta 2do Orden')
plt.plot(xi4,yi4, label='Runge-Kutta 4do Orden')
plt.title('Deflexión de mástil')
plt.xlabel('x')
plt.ylabel('y: deflexión')
plt.legend()
plt.grid()
plt.show()

s1Eva_IT2015_T4 Lingotes metales

Ejercicio: 1Eva_IT2015_T4 Lingotes metales

Se plantea que cada lingote debe aportar una proporción xi al lingote nuevo a ser fundido.

Se dispone de los compuestos de cada lingote por filas:

compuesto = np.array([[ 20, 50, 20, 10],
                      [ 30, 40, 10, 20],
                      [ 20, 40, 10, 30],
                      [ 50, 20, 20, 10]])
proporcion = np.array([ 27, 39.5, 14, 19.5])

El contenido final de cada componente basado en los aportes xi de cada lingote para cada componente.

Ejemplo para los 27 gramos de oro

20x_1 + 30x_2+ 20x_3 + 50x_4 = 27

se realiza el mismo procedimiento para los otros tipos de metal.

50x_1 + 40x_2+ 40x_3 + 20x_4 = 39.5 20x_1 + 10x_2+ 10x_3 + 20x_4 = 14 10x_1 + 20x_2+ 30x_3 + 10x_4 = 19.5

Las ecuaciones se escriben en la forma matricial Ax=B

\begin{bmatrix} 20 && 30&& 20 &&50 \\ 50 && 40 && 40 && 20 \\ 20 && 10 && 10 && 20 \\ 10 && 20 && 30 && 10 \end{bmatrix} = \begin{bmatrix} x_1 \\ x_2 \\x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} 27 \\ 39.5 \\ 14 \\ 19.5 \end{bmatrix}

Para resolver se plantea la matriz aumentada

\begin{bmatrix} 20 && 30&& 20 &&50 && 27\\ 50 && 40 && 40 && 20 && 39.5\\ 20 && 10 && 10 && 20 && 14 \\ 10 && 20 && 30 && 10 && 19.5 \end{bmatrix}

se pivotea por filas la matriz:

\begin{bmatrix} 50 && 40 && 40 && 20 && 39.5\\ 20 && 30&& 20 &&50 && 27\\ 10 && 20 && 30 && 10 && 19.5 \\ 20 && 10 && 10 && 20 && 14 \end{bmatrix}

Para eliminar hacia adelante:

\begin{bmatrix} 50 && 40 && 40 && 20 && 39.5 \\ 20 - 50\frac{20}{50} && 30-40\frac{20}{50} && 20-40\frac{20}{50} && 50-20\frac{20}{50} && 27-39.5\frac{20}{50}\\ 10 - 50\frac{10}{50} && 20-40\frac{10}{50} && 30-40\frac{10}{50} && 10-20\frac{10}{50} && 19.5-39.5\frac{10}{50} \\ 20 - 50\frac{20}{50} && 10-40\frac{20}{50} && 10-40\frac{20}{50} && 20-20\frac{20}{50} && 14-39.5\frac{20}{50} \end{bmatrix}

continuando con el desarrollo:

Elimina hacia adelante
[[50.  40.  40.  20.  39.5]
 [ 0.  14.   4.  42.  11.2]
 [ 0.  12.  22.   6.  11.6]
 [ 0.  -6.  -6.  12.  -1.8]]
Elimina hacia adelante
[[ 50.  40.  40.      20.  39.5]
 [  0.  14.   4.      42.  11.2]
 [  0.   0.  18.5714 -30.   2. ]
 [  0.   0.  -4.2857  30.   3. ]]
Elimina hacia adelante
[[ 50.  40.  40.      20.     39.5   ]
 [  0.  14.   4.      42.     11.2   ]
 [  0.   0.  18.5714 -30.      2.    ]
 [  0.   0.   0.      23.0769  3.4615]]
Elimina hacia adelante
[[ 50.  40.   40.      20.     39.5  ]
 [  0.  14.    4.      42.     11.2  ]
 [  0.   0.   18.5714 -30.      2.   ]
 [  0.   0.    0.      23.0769  3.4615]]
Elimina hacia atras
[[ 1.    0.    0.    0.    0.25]
 [ 0.    1.    0.    0.    0.25]
 [ 0.    0.    1.   -0.    0.35]
 [ 0.    0.    0.    1.    0.15]]
el vector solución X es:
[[0.25]
 [0.25]
 [0.35]
 [0.15]]
verificar que A.X = B
[[39.5]
 [27. ]
 [19.5]
 [14. ]]

Las proporciones de cada lingote a usar para el nuevo lingote que cumple con lo solicitado son:

[0.25, 0.25, 0.35, 0.15]


Algoritmo en python

usado para la solución es:

# Método de Gauss-Jordan ,
# Recibe las matrices A,B
# presenta solucion X que cumple: A.X=B
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
A = np.array([[ 50., 40, 40, 20],
              [ 20., 30, 20, 50],
              [ 10., 20, 30, 10],
              [ 20., 10, 10, 20]])
B1 = np.array([ 39.5, 27, 19.5, 14])

B = np.transpose([B1])

# PROCEDIMIENTO
casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

print('matriz aumentada: ')
print(AB)
# Gauss elimina hacia adelante
# tarea: verificar términos cero
for i in range(0,n,1):
    pivote = AB[i,i]
    adelante = i+1 
    for k in range(adelante,n,1):
        if (np.abs(pivote)>=casicero):
            factor = AB[k,i]/pivote
            AB[k,:] = AB[k,:] - factor*AB[i,:]
    print('Elimina hacia adelante')
    print(AB)

# Gauss-Jordan elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    # Normaliza a 1 elemento diagonal
    AB[i,:] = AB[i,:]/AB[i,i]
    pivote = AB[i,i] # uno
    # arriba de la fila i
    atras = i-1 
    for k in range(atras,0-1,-1):
        if (np.abs(AB[k,i])>=casicero):
            factor = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*factor - AB[i,:]
        else:
            factor= 'division para cero'
print('Elimina hacia atras')
print(AB)

X = AB[:,ultcolumna]

# Verifica resultado
verifica = np.dot(A,X)

# SALIDA
print('el vector solución X es:')
print(np.transpose([X]))

print('verificar que A.X = B')
print(np.transpose([verifica]))

Tarea: Revisar sobre la última pregunta.