s2Eva_IT2012_T1_MN Longitud de teleférico

Ejercicio: 2Eva_IT2012_T1_MN Longitud de teleférico

Los datos tomados para el problema son:

x    = [0.00, 0.25, 0.50, 0.75, 1.00]
f(x) = [25.0, 22.0, 32.0, 51.0, 75.0]

Se debe considerar que los datos tienen tamaño de paso (h) del mismo valor.

Literal a)

Fórmulas de orden 2, a partir de:

http://blog.espol.edu.ec/analisisnumerico/formulas-de-diferenciacion-por-diferencias-divididas/

considere que el término del Error O(h2), perdió una unidad del exponente en el proceso, por lo que las fórmulas de orden 2 tienen un error del mismo orden.

Se puede usar por ejemplo:

Para los términos x en el intervalo [0,0.50] hacia adelante

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

Para el término x con 0.75, centradas:

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

y para el término x con 1.0, hacia atras:

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

Luego se aplica el resultado en la fórmula:

g(x) = \sqrt{1+[f'(x)]^2}

L = \int_a^b g(x) \delta x .

Literal b)

Use las fórmulas de integración numérica acorde a los intervalos. Evite repetir intervalos, que es el error más común.

Por ejemplo, se puede calcular el integral de g(x) aplicando dos veces Simpson de 1/3, que sería la más fácil de aplicar dado los h iguales.

Otra opción es Simpson de 3/8 añadido un trapecio, otra forma es solo con trapecios en todo el intervalo.

Como ejemplo de cálculo usando un algoritmo en Python se muestra que:

f'(x): [-38.  22.  66.  86. 106.]
 g(x): [ 38.0131  22.0227  66.0075  86.0058 106.0047]
L con S13:  59.01226169578733
L con trapecios:  61.511260218050175

los cálculos fueron realizados a partir de la funciones desarrolladas durante la clase. Por lo que se muestran 3 de ellas en el algoritmo.

import numpy as np
import matplotlib.pyplot as plt

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

def isimpson13(f,h):
    n = len(f)
    integral = 0
    for i in range(0,n-2,2):
        area = (h/3)*(f[i]+4*f[i+1]+f[i+2])
        integral = integral + area
    return(integral)

def isimpson38 (f,h):
    n=len(f)
    integral=0
    for i in range(0,n-3,3):
        area=(3*h/8)*(f[i]+3*f[i+1]+3*f[i+2] +f[i+3] )
        integral=integral + area
    return(integral)

# INGRESO
x = np.array( [0.0,0.25,0.50,0.75,1.00])
fx= np.array([ 25.0, 22.0, 32.0, 51.0, 75.0])

# PROCEDIMIENTO
n = len(fx)
dx = x[1]-x[0]

# Diferenciales
dy = np.zeros(n)

for i in range(0,n-2,1):
    dy[i] = (-fx[i+2]+4*fx[i+1]-3*fx[i])/(2*dx)
# Diferenciales penúltimo
i = n-2
dy[i] = (fx[i+1]-fx[i-1])/(2*dx)
# Diferenciales último
i = n-1
dy[i] = (3*fx[i]-4*fx[i-1]+fx[i-2])/(2*dx)

# Función gx
gx = np.sqrt(1+(dy**2))

# Integral
integral = isimpson13(gx,dx)
integrartr = itrapecio(gx,dx)

# SALIDA 
print('f\'(x):', dy)
print(' g(x):', gx)
print("L con S13: ", integral )
print("L con trapecios: ", integrartr )

plt.plot(x,fx)
plt.show()

La gráfica del cable es:

s2Eva_IT2012_T2 Modelo de clima

Ejercicio: 2Eva_IT2012_T2 Modelo de clima

Se deja de tarea realizar las tres primeras iteraciones en papel.

Se presenta el resultado usando el algoritmo de segundo grado como una variante a la respuesta usada como ejemplo.

 [ ti, xi, yi, zi]
[[ 0.0000e+00  1.0000e+01  7.0000e+00  7.0000e+00]
 [ 2.5000e-03  9.9323e+00  7.5033e+00  7.1335e+00]
 [ 5.0000e-03  9.8786e+00  7.9988e+00  7.2774e+00]
 ...
 [ 2.4995e+01 -8.4276e+00 -2.7491e+00  3.3021e+01]
 [ 2.4998e+01 -8.2860e+00 -2.6392e+00  3.2858e+01]
 [ 2.5000e+01 -8.1453e+00 -2.5346e+00  3.2692e+01]]

Algoritmo en Python

# 2Eva_IT2012_T2 Modelo de clima
# MATG1013 Análisis Numérico
import numpy as np
import matplotlib.pyplot as plt

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


#INGRESO
to = 0
xo = 10
yo = 7
zo = 7
f = lambda t,x,y,z: 10*(y-x)
g = lambda t,x,y,z: x*(28-z) - y
j = lambda t,x,y,z: -(8/3)*z + (x*y)
h = 0.0025
muestras = 10000

#PROCEDIMIENTO
#Rugen-Kutta 2_orden
tabla = rungekutta2_fg(f,g,j,to,xo,yo,zo,h,muestras)

#SALIDA
np.set_printoptions(precision=4)
print(' [ ti, xi, yi, zi]')
print(tabla)

# Gráfica
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(tabla[:,1], tabla[:,2], tabla[:,3])
plt.show()

s2Eva_IIT2011_T3 EDP Parabólica, explícito

Ejercicio: 2Eva_IIT2011_T3 EDP Parabólica, explícito

La ecuación a resolver es:

\frac{\delta u}{\delta t} - \frac{\delta^2 u}{\delta x^2} =2

Como el método requerido es explícito se usan las diferencias divididas:

\frac{d^2 u}{dx^2} = \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2} \frac{du}{dt} = \frac{u_{i,j+1} - u_{i,j} }{\Delta t}

Reordenando la ecuación al modelo realizado en clase:

\frac{\delta^2 u}{\delta x^2} = \frac{\delta u}{\delta t} - 2 \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2} = \frac{u_{i,j+1} - u_{i,j} }{\Delta t} - 2

multiplicando cada lado por Δt

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

Se establece el valor de

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

obteniendo la ecuación general, ordenada por índices de izquierda a derecha:

\lambda u_{i-1,j} +(1- 2 \lambda)u_{i,j} + \lambda u_{i+1,j} + 2\Delta t = u_{i,j+1}

con valores de:

λ = Δt/(Δx)2 = (0.04)/(0.25)2 = 0.64
P = λ = 0.64
Q = (1-2λ) = -0.28
R = λ = 0.64
0.64 u_{i-1,j} - 0.28 u_{i,j} + 0.64 u_{i+1,j} + 0.08 = u_{i,j+1}

Se realizan 3 iteraciones:

i= 1, j=0

u_{1,1} = 0.64 u_{0,0} -0.28u_{1,0}+0.64 u_{2,0}+0.08 u_{1,1} = 0.64 [\sin(\pi*0)+0*(1-0)]- 0.28[\sin(\pi*0.25)+0.25*(1-0.25)] +0.64[\sin(\pi*0.5)+ 0.5*(1-0.5)]+0.08

u[1,1] =0.89

i = 2, j=0

0.64 u_{1,0} - 0.28 u_{2,0} + 0.64 u_{3,0} + 0.08 = u_{2,1}

u[1,0] = 1.25

i = 3, j=0

0.64 u_{2,0} - 0.28 u_{3,0} + 0.64 u_{4,0} + 0.08 = u_{3,1}

u[3,1] = 0.89


Algoritmo en Python

con los valores y ecuación del problema

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito, usando diferencias finitas
# Referencia: Chapra 30.2 p.888 pdf.912
#       Rodriguez 10.2 p.406
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 0
Tb = 0
T0 = lambda x: np.sin(np.pi*x)+x*(1-x)
# longitud en x
a = 0
b = 1
# Constante K
K = 1
# Tamaño de paso
dx = 0.1
dt = dx/10/2
# 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(xi[1:ultimox])
u[ultimox,j] = Tb

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

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

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

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

s2Eva_IT2010_T3 EDP elíptica, Placa no rectangular

Ejercicio: 2Eva_IT2010_T3 EDP elíptica, Placa no rectangular

la fórmula a resolver, siendo f(x,y)=20

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

Siendo las condiciones de frontera en los bordes marcados:

Se convierte a forma discreta la ecuación

\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} = 20 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) = 20 \Delta x^2

siendo λ = 1, al tener la malla en cuadrícula Δx=Δy

u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1}+u_{i,j-1} = 20 \Delta x^2

despejando u[i,j]

u_{i,j} =\frac{1}{4}\Big(u_{i+1,j}+u_{i-1,j} + u_{i,j+1}+u_{i,j-1} - 20 \Delta x^2 \Big)

Para el ejercicio, se debe usar las fronteras para los valores de cálculo que se puede hacer de dos formas:
1. Por posición del valor en la malla [i,j]
2. Por valor xi,yj que es la forma usada cuando la función es contínua.

1. Por posición del valor en la malla

Se busca la posición de la esquina derecha ubicada en xi = 0 y yj=1.5 y se la ubica como variable ‘desde’ como referencia para ubicar los valores de las fronteras usando los índices i, j.

# valores en fronteras
# busca poscición de esquina izquierda
desde = -1 
for j in range(0,m,1):
    if ((yj[j]-1.5)<tolera):
        desde = j
# llena los datos de bordes de placa
for j  in range(0,m,1):
    if (yj[j]>=1):
        u[j-desde,j] = Ta
    u[n-1,j] = Tb(xi[n-1],yj[j])
for i in range(0,n,1):
    if (xi[i]>=1):
        u[i,m-1]=Td
    u[i,desde-i] = Tc(xi[i],yj[i])
# valor inicial de iteración

2. Por valor xi, yj

Se establecen las ecuaciones que forman la frontera:

ymin = lambda x,y: 1.5-(1.5/1.5)*x
ymax = lambda x,y: 1.5+(1/1)*x

que se usan como condición para hacer el cálculo en cada nodo

if ((yj[j]-yjmin)>tolera and (yj[j]-yjmax)<-tolera):
                u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-20*(dx**2))/4

Para minimizar errores por redondeo o truncamiento al seleccionar los puntos, la referencia hacia cero se toma como «tolerancia»; en lugar de más que cero o menos que cero.

Con lo que se obtiene el resultado mostrado en la gráfica aumentando la resolución con Δx=Δy=0.05:

converge =  1
xi=
[0.   0.05 0.1  0.15 0.2  0.25 0.3  0.35 0.4  0.45 0.5  0.55 0.6  0.65
 0.7  0.75 0.8  0.85 0.9  0.95 1.   1.05 1.1  1.15 1.2  1.25 1.3  1.35
 1.4  1.45 1.5 ]
yj=
[0.   0.05 0.1  0.15 0.2  0.25 0.3  0.35 0.4  0.45 0.5  0.55 0.6  0.65
 0.7  0.75 0.8  0.85 0.9  0.95 1.   1.05 1.1  1.15 1.2  1.25 1.3  1.35
 1.4  1.45 1.5  1.55 1.6  1.65 1.7  1.75 1.8  1.85 1.9  1.95 2.   2.05
 2.1  2.15 2.2  2.25 2.3  2.35 2.4  2.45 2.5 ]
matriz u
[[  0.     0.     0.   ...   0.     0.     0.  ]
 [  0.     0.     0.   ...   0.     0.     0.  ]
 [  0.     0.     0.   ...   0.     0.     0.  ]
 ...
 [  0.     0.     6.67 ...  96.1   98.03 100.  ]
 [  0.     3.33   5.31 ...  96.03  98.   100.  ]
 [  0.     2.     4.   ...  96.    98.   100.  ]]
>>>

Algoritmo en Python

# Ecuaciones Diferenciales Parciales
# Elipticas. Método iterativo para placa NO rectangular
import numpy as np

# INGRESO
# Condiciones iniciales en los bordes
Ta = 100
Tb = lambda x,y:(100/2.5)*y
Tc = lambda x,y: 100-(100/1.5)*x
Td = 100
# dimensiones de la placa no rectangular
x0 = 0
xn = 1.5
y0 = 0
yn = 2.5
ymin = lambda x,y: 1.5-(1.5/1.5)*x
ymax = lambda x,y: 1.5+(1/1)*x
# discretiza, supone dx=dy
dx = 0.05 
dy = 0.05 
maxitera = 100
tolera = 0.0001

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)
# Matriz u
u = np.zeros(shape=(n,m),dtype = float)

# valores en fronteras
# busca posición de esquina izquierda
desde = -1 
for j in range(0,m,1):
    if ((yj[j]-1.5)<tolera): desde = j 

# llena los datos de bordes de placa 

for j in range(0,m,1): if (yj[j]>=1):
        u[j-desde,j] = Ta
    u[n-1,j] = Tb(xi[n-1],yj[j])
for i in range(0,n,1):
    if (xi[i]>=1):
        u[i,m-1]=Td
    u[i,desde-i] = Tc(xi[i],yj[i])
# valor inicial de iteración
# se usan los valores cero

# iterar puntos interiores
itera = 0
converge = 0
while not(itera>=maxitera and converge==1):
    itera = itera + 10000
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            yjmin = ymin(xi[i],yj[j])
            yjmax = ymax(xi[i],yj[j])
            if ((yj[j]-yjmin)>tolera and (yj[j]-yjmax)<-tolera):
                u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-20*(dx**2))/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))

    if (erroru<tolera):
        converge=1

# SALIDA
np.set_printoptions(precision=2)
print('converge = ', converge)
print('xi=')
print(xi)
print('yj=')
print(yj)
print('matriz u')
print(u)

Para incorporar la gráfica en forma de superficie.

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

X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()
ax = Axes3D(figura)
ax.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()

s3Eva_IIT2010_T4 EDO con Taylor

Ejercicio: 3Eva_IIT2010_T4 EDO con Taylor

\frac{\delta y}{\delta x} = \frac{y^3}{1-2xy^2} y(0) = 1, 0 \leq x \leq 1

escrita en forma simplificada

y' = \frac{y^3}{1-2xy^2}

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

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

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

y_{i+1} = y_i +\frac{h}{1!}

Con lo que se puede realizar el algoritmo

estimado[xi,yi]
[[0.  1. ]
 [0.2 1.2 ]
 [0.4 2.01509434]
 [0.6 1.28727044]
 [0.8 0.85567954]
 [1.  0.12504631]]

Algoritmo en Python

# 3Eva_IIT2010_T4 EDO con Taylor
# estima la solución para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor2t(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        y = y + h*d1y(x,y)
        x = x+h
        estimado[i] = [x,y]
    return(estimado)

# PROGRAMA PRUEBA
# 3Eva_IIT2010_T4 EDO con Taylor

# INGRESO.
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: (y**3)/(1-2*x*(y**2))
x0 = 0
y0 = 1
h = 0.2
muestras = 5

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

# SALIDA
print('estimado[xi,yi]')
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 2 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

Nota: Revisar los resultados no lineales con los valores de h=0.02

s3Eva_IT2010_T2 EDO problema con valor inicial dy/dx

Ejercicio: 3Eva_IT2010_T2 EDO problema con valor inicial dy/dx

Ecuación del problema:

(2y^2 + 4x^2)\delta x -xy \delta y =0

se despeja dy/dx:

(2y^2 + 4x^2)\delta x = xy \delta y \frac{\delta y}{\delta x} =\frac{2y^2 + 4x^2}{xy}

con valores iniciales de x0 = 1, y0=-2 , h=0.2 e intervalo [1,2]

Usando Runge-Kutta de 2do Orden

Iteración 1

K_1= h\frac{\delta y}{\delta x}(1,-2) =(0.2)\frac{2(-2)^2 + 4(1)^2}{(1)(-2)}= -1.2 K_2= h\frac{\delta y}{\delta x}(1+0.2,-2+(-1.2)) =(0.2)\frac{2(-2-1.2)^2 + 4(1+0.2)^2}{(1+0.2)(-2-1.2)}= -1.3667 y_1 = -2 + \frac{-1.2-1.3667}{2} = -3.2833 x_1 = x_0 + h = 1 + 0.2 = 1.2 error = O(h^3) = O(0.2^3) = 0.008

Iteración 2

K_1= h\frac{\delta y}{\delta x}(1.2,-3.2833) =(0.2)\frac{2(-3.2833)^2 + 4(1.2)^2}{(1.2)(-3.2833)}= -1.3868 K_2= h\frac{\delta y}{\delta x}(1.2+0.2,-3.2833+(-1.3868)) =(0.2)\frac{2(-3.2833+(-1.3868))^2 + 4(1.2+0.2)^2}{(1.2+0.2)(-3.2833+(-1.3868))}= -1.5742 y_2 = -3.2833 + \frac{-1.3868-1.5742}{2} = -4.7638 x_2 = x_1 + h = 1.2 + 0.2 = 1.4

Iteración 3

K_1= h\frac{\delta y}{\delta x}(1.4,-4.7638) =(0.2)\frac{2(-4.76383)^2 + 4(1.4)^2}{(1.4)(-4.7638)}= -1.5962 K_2= h\frac{\delta y}{\delta x}(1.4+0.2,-4.7638+(-1.5962)) =(0.2)\frac{2(-4.7638+(-1.5962))^2 + 4(1.4+0.2)^2}{(1.4+0.2)(-4.7638+(-1.5962))}= -1.7913 y_3 = -4.7638 + \frac{-1.5962-1.7913}{2} = -6.4576 x_3 = x_2 + h = 1.4 + 0.2 = 1.6

con lo que usando el algoritmo se obtiene la tabla y gráfica:

 [xi,       yi,       K1,     K2     ]
[[  1.      -2.       0.       0.    ]
 [  1.2     -3.2833  -1.2     -1.3667]
 [  1.4     -4.7638  -1.3868  -1.5742]
 [  1.6     -6.4576  -1.5962  -1.7913]
 [  1.8     -8.3698  -1.8126  -2.0119]
 [  2.     -10.5029  -2.032   -2.2342]

Algoritmo en Python

# 3Eva_IT2010_T2 EDO
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
d1y = lambda x,y : (2*(y**2)+4*(x**2))/(x*y)
x0 = 1
y0 = -2
h = 0.2
a = 1
b = 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_IT2010_T1 Envase cilíndrico

Ejercicio: 3Eva_IT2010_T1 Envase cilíndrico

Se conoce que el volumen del cilindro es 1000 cm3
que se calcula como:

Volumen = area_{base} . altura = \pi r^2 h \pi r^2 h = 1000 h = \frac{1000}{\pi r^2}

con lo que la altura queda en función del radio, pues el volumen es una constante.

Para conocer el área de la hoja de material para el envase se conoce que las tapas cilíndricas deben ser 0.25 mas que el radio para sellar el recipiente

Area de una tapa circular:

Area_{tapa} = \pi (r+ 0.25)^2

para la hoja lateral:

Area_{lateral} = base * altura = (2 \pi r + 0.25) h

que en función del radio, usando la fórmula para h(r) se convierte en:

Area_{lateral} = (2 \pi r + 0.25) \frac{1000}{\pi r^2}

por lo que el total de material de hoja a usar corresponde a dos tapas circulares y una hoja lateral.

Area total = 2 \pi (r+ 0.25)^2 + (2 \pi r + 0.25) \frac{1000}{\pi r^2}

la gráfica permite observar el comportamiento entre el área de la hoja necesaria para fabricar el envase y el radio de las tapas circulares.

El objetivo es encontrar el área mínima para consumir menos material. Con la fórmula mostrada se puede aplicar uno de los métodos para encontrar raíces a partir de la derivada de la fórmula.

Tarea: Encontrar el valor requerido para r con la tolerancia indicada para el examen.


Instrucciones para la gráfica:

# 3Eva_IT2010_T1 Envase cilíndrico
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
a = 0.5
b = 20
muestras = 51

Area = lambda r: (2*np.pi)*(r+0.25)**2 +(2*np.pi*r+0.25)*1000/(np.pi*(r**2))

# PROCEDIMIENTO
ri = np.linspace(a,b,muestras)
Areai = Area(ri)

# SALIDA
plt.plot(ri,Areai)
plt.xlabel('r (cm)')
plt.ylabel('Area (cm2)')
plt.title('Area hoja para material cilindro')
plt.show()

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