La gráfica usando una mayor resolución para tener una idea de la solución:
Los resultados se obtienen usando las siguientes instrucciones:
# 2da Evaluación I Término 2018# Tema 3. EDP Elipticaimport numpy as np
# INGRESO# ejes x,y
x0 = 0 ; xn = 1 ; hx = (1/3)# (1/3)/10
y0 = 0 ; yn = 1 ; hy = (1/3) # (1/3)/10# Fronteras
fux0 = lambda x: x+1
fu0y = lambda y: y+1
fux1 = lambda x: x**2 + x + 2
fu1y = lambda y: y**2 + y + 2
fxy = lambda x,y: 2*(x**2+y**2)
# PROCEDIMIENTO
xi = np.arange(x0,xn+hx,hx)
yj = np.arange(y0,yn+hy,hy)
n = len(xi)
m = len(yj)
# funcion f[xi,yi]
fij = np.zeros(shape=(n,m), dtype = float)
for i inrange(0,n,1):
for j inrange(0,m,1):
fij[i,j]=fxy(xi[i],yj[j])
# matriz inicial u[i,j]
u = np.zeros(shape=(n,m), dtype = float)
u[:,0] = fux0(xi)
u[0,:] = fu0y(yj)
u[:,m-1] = fux1(xi)
u[n-1,:] = fu1y(yj)
uinicial = u.copy()
# Calcular de forma iterativa
maxitera = 100
tolera = 0.0001
# valor inicial de iteración
promedio = (np.max(u)+np.min(u))/2
u[1:n-1,1:m-1] = promedio
# iterar puntos interiores
itera = 0
converge = 0
erroru = 2*tolera # para calcular al menos una matrizwhilenot(erroru=maxitera):
itera = itera +1
nueva = np.copy(u)
for i inrange(1,n-1):
for j inrange(1,m-1):
u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-(hy**2)*fij[i,j])/4
diferencia = nueva-u
erroru = np.linalg.norm(np.abs(diferencia))
if (erroru<tolera):
converge=1
# SALIDA
np.set_printoptions(precision=2)
print('x[i]:')
print(xi)
print('y[j]:')
print(yj)
print('f[i,j]:')
print(fij)
print('matriz inicial u[i,j]:')
print(uinicial)
print('resultado para u, iterando: ')
print('converge = ', converge)
print('iteraciones = ', itera)
print(u)
para obtener la gráfica se debe añadir:
# Gráficaimport matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
X, Y = np.meshgrid(xi, yj)
figura = plt.figure()
ax = Axes3D(figura)
U = np.transpose(u) # ajuste de índices fila es x
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()
Para el problema, se usan los puntos: [a,f(a)], [b,f(b)] y [c,f(c)]
por donde pasa la curva f(x) aproximada a un polinomio de grado 2, f(x) \approx p(x)
# 2da Evaluación I Término 2018# Tema 4. canal acceso a Puertos de Guayaquilimport numpy as np
# INGRESO
profcanal = 11
xi = np.array([ 0., 50., 100., 150., 200., 250.])
yi = np.array([ 0., 100., 200., 300., 400., 500.])
batimetria = [[ -6.79,-12.03,-10.04,-11.60, -7.24,-7.91],
[ -8.85,-10.89, -8.95, -7.23,-11.42,-7.93],
[-11.90, -9.86, -9.35,-12.05, -9.38,-9.65],
[ -7.30,-11.55,-10.41, -8.67,-11.84,-6.77],
[-12.17, -9.62, -7.47, -6.51, -9.02,-9.60],
[-11.90,-10.23,-10.68, -9.94, -6.76,-7.46]]
batimetria = np.array(batimetria)
# PROCEDIMIENTO
[n,m] = np.shape(batimetria)
# Matriz remover sedimentos
remover = batimetria + profcanal
for i inrange(0,n,1):
for j inrange(0,m,1):
if remover[i,j]<0:
remover[i,j]=0
# SALIDAprint('matriz remover sedimentos: ')
print(remover)
b) el volumen se calcula usando el algoritmo de Simpson primero por un eje, y luego con el resultado se continúa con el otro eje,
Considere que existen 6 puntos, o 5 tramos integrar en cada eje.
Al usar Simpson de 1/3 que usan tramos pares, faltaría integrar el último tramo.
En el caso de Simpson de 3/8 se requieren tramos múltiplos de 3, por lo que faltaría un tramo para volver a usar la fórmula.
La solución por filas se implementa usando una combinación de Simpson 3/8 para los puntos entre remover[i, 0:(2+1)] y Simpson 1/3 para los puntos entre remover[i, 3:(4+1)].
Luego se completa el integral del otro eje con el resultado anterior, aplicando el mismo método.
Se obtuvieron los siguientes resultados:
Integral en eje x:
[435.1 346.5 287.44 255.58 590.54 440.52]
Volumen: 199161.11111111115
que se obtiene usando las instrucciones a continuación de las anteriores:
a. Planteamiento con Runge-Kutta 2do Orden para Segunda derivada
La expresión:
\frac{dv}{dt} = g - \frac{cd}{m} v^2
se puede plantear sustituir la variable con v = -\frac{dy}{dt} al considerar el sentido contrario entre la velocidad al caer y la referencia de altura hacia arriba. Ver figura.
\frac{d^2y}{dt^2} = g - \frac{cd}{m} \Big( \frac{dy}{dt} \Big) ^2
con las condiciones iniciales del ejercicio t0 = 0 , y0 = 1000, v0 = 0
la velocidad se inicia con cero, si el paracaidista se deja caer desde el borde el risco, como en el video adjunto al enunciado.
Para las iteraciones, recuerde que
t se corrige con t+h (en el algoritmo era la posición para x)
y se corrige con y+K1y
v se corrige con v+K1v (en el algoritmo era la posición para z)
itera = 0
K1y = h f(t,y,v) = 2(-(0)) = 0 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (0)^2) = 19.6
.. K2y = h f(t+h, y+K1y, v + K1v) = 2(-(0 + 19.6)) = -39.2
K2v = h g(t+h, y+K1y, v + K1v) = 2(9.8 - \frac{0.225}{90} (0+19.6)^2) =17.6792
Resultados
La velocidad máxima, si no hubiese límite en la altura, se encuentra en alrededor de 62.39 m/s.
Sin embargo, luego de 20 segundos se observa que la altura alcanza el valor de cero, es decir se alcanzó tierra con velocidad de 62 m/s, que son algo mas de 223 Km/h, el objeto se estrella...!
# 2Eva_IT2018_T1 Paracaidista wingsuitimport numpy as np
import matplotlib.pyplot as plt
defrungekutta2_fg(f,g,x0,y0,z0,h,muestras,
vertabla=False, precision = 6):
''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
x0,y0 son valores iniciales, h es tamano de paso,
muestras es la cantidad de puntos a calcular.
'''
tamano = muestras + 1
tabla = np.zeros(shape=(tamano,3+4),dtype=float)
# incluye el punto [x0,y0,z0]
tabla[0] = [x0,y0,z0,0,0,0,0]
xi = x0
yi = y0
zi = z0
i=0
if vertabla==True:
print('Runge-Kutta Segunda derivada')
print('i ','[ xi, yi, zi',']')
print(' [ K1y, K1z, K2y, K2z ]')
np.set_printoptions(precision)
print(i,tabla[i,0:3])
print(' ',tabla[i,3:])
for i inrange(1,tamano,1):
K1y = h * f(xi,yi,zi)
K1z = h * g(xi,yi,zi)
K2y = h * f(xi+h, yi + K1y, zi + K1z)
K2z = h * g(xi+h, yi + K1y, zi + K1z)
yi = yi + (K1y+K2y)/2
zi = zi + (K1z+K2z)/2
xi = xi + h
tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
if vertabla==True:
txt = ' 'if i>=10:
txt=' 'print(str(i)+'',tabla[i,0:3])
print(txt,tabla[i,3:])
return(tabla)
# PROGRAMA PRUEBA# INGRESO
f = lambda t,y,v: -v # el signo, revisar diagrama cuerpo libre
g = lambda t,y,v: 9.8 - (0.225/90)*(v**2)
t0 = 0
y0 = 1000
v0 = 0
h = 2
muestras = 15+1
# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,y0,v0,h,muestras, vertabla=True)
ti = tabla[:,0]
yi = tabla[:,1]
vi = tabla[:,2]
# SALIDA# print('tabla de resultados')# print(tabla)# GRAFICA
plt.subplot(211)
plt.plot(ti,vi,label='velocidad v(t)', color='green')
plt.plot(ti,vi,'o')
plt.ylabel('velocidad (m/s)')
plt.title('paracaidista Wingsuit con Runge-Kutta')
plt.legend()
plt.grid()
plt.subplot(212)
plt.plot(ti,yi,label='Altura y(t)',)
plt.plot(ti,yi,'o',)
plt.axhline(0, color='red')
plt.xlabel('tiempo (s)')
plt.ylabel('Altura (m)')
plt.legend()
plt.grid()
plt.show()
b. Usando Runge-Kutta 2do Orden para Primera Derivada o velocidad,
El problema para un tiempo de observación t>0, se puede dividir en dos partes: velocidad y altura.
Determinar velocidad: Se aplica Runge-Kutta a la expresión con primera derivada o velocidad. Use los valores iniciales dados, descarte calcular las alturas.
Determinar las altura: Con los valores de velocidades y la altura inicial de 1 km = 1000 m puede integrar con trapecio para obtener la tabla de alturas. Se integra tramo a tramo.
Observe las unidades de medida y que la velocidad es contraria al eje de altura dy/dt = -v.
La expresión:
\frac{dv}{dt} = g - \frac{cd}{m} v^2 f(t,v) = g - \frac{cd}{m} v^2
las siguientes iteraciones se completan con el algoritmo.
Resultados
La velocidad máxima, si no hubiese límite en la altura, se encuentra en alrededor de 62.39 m/s.
Sin embargo, luego de 20 segundos se observa que la altura alcanza el valor de cero, es decir se alcanzó tierra con velocidad de 62 m/s, que son algo mas de 223 Km/h, el objeto se estrella...!
Los cálculos se realizaron usando las instrucciones en Python:
# 2da Evaluación I Término 2018# Tema 1. Paracaidista wingsuitimport numpy as np
defrungekutta2(d1y,x0,y0,h,muestras):
# Runge Kutta de 2do orden
tamano = muestras + 1
tabla = np.zeros(shape=(tamano,2+2),dtype=float)
# incluye el punto [x0,y0]
tabla[0] = [x0,y0,0,0]
xi = x0
yi = y0
for i inrange(1,tamano,1):
K1 = h * d1y(xi,yi)
K2 = h * d1y(xi+h, yi + K1)
yi = yi + (1/2)*(K1+K2)
xi = xi + h
tabla[i] = [xi,yi,K1,K2]
return(tabla)
defintegratrapecio_fi_tabla(xi,fi,y0):
tamano = len(xi)
yi = np.zeros(tamano,dtype = float)
yi[0] = y0
for i inrange(1,tamano,1):
h = xi[i]-xi[i-1]
trapecio = h*(fi[i]+fi[i-1])/2
yi[i]= yi[i-1] + trapecio
return(yi)
# PROGRAMA -------------------------# INGRESO
g = 9.8
cd = 0.225
m = 90
d1v = lambda t,v: g - (cd/m)*(v**2)
t0 = 0
v0 = 0
h = 2
y0 = 1000
muestras = 11
# PROCEDIMIENTO
velocidad = rungekutta2(d1v,t0,v0,h,muestras)
ti = velocidad[:,0]
vi = velocidad[:,1]
# Altura, velocidad es contraria altura,# integrar en cada tramo por trapecios o Cuadratura de Gauss
altura = integratrapecio_fi_tabla(ti,-vi,y0)
# Tabla de resultados de tiempo, velocidad, altura
altura = np.transpose([altura])
tabla = np.concatenate((velocidad[:,0:2],altura), axis = 1)
# SALIDA
np.set_printoptions(precision=4)
print('velocidad con Runge-Kutta primera derivada')
print(' [tiempo, velocidad, K1,K2]')
print(velocidad)
np.set_printoptions(precision=2)
print('velocidad con RK2 y altura con trapecio')
print(' [tiempo, velocidad, altura]')
print(tabla)
Runge kutta de 2do Orden
f: y' = z
g: z' = .....
K1y = h f(xi, yi, zi)
K1z = h g(xi, y1, zi)
K2y = h f(xi+h, yi+K1y, zi+K1z)
K2z = h g(xi+h, yi+K1y, zi+K1z)
y(i+1) = yi + (1/2)(K1y + K2y)
z(i+1) = zi + (1/2)(K1z + K2z)
x(i+1) = xi + h
E = O(h3)
xi ≤ z ≤ x(i+1)
Instruccciones en python, usando el algoritmo desarrollado en clase
# Runge Kutta de 2do# EDO de 2do orden con condiciones de inicioimport numpy as np
import matplotlib.pyplot as plt
defrungekutta2_fg(f,g,v0,h,m):
tabla = [v0]
xi = v0[0]
yi = v0[1]
zi = v0[2]
for i inrange(0,m,1):
K1y = h * f(xi,yi,zi)
K1z = h * g(xi,yi,zi)
K2y = h * f(xi+h, yi + K1y, zi+K1z)
K2z = h * g(xi+h, yi + K1y, zi+K1z)
yi1 = yi + (1/2)*(K1y+K2y)
zi1 = zi + (1/2)*(K1z+K2z)
xi1 = xi + h
vector = [xi1,yi1,zi1]
tabla.append(vector)
xi = xi1
yi = yi1
zi = zi1
tabla = np.array(tabla)
return(tabla)
# Programa Prueba# Funciones
f = lambda x,y,z : z
g = lambda x,y,z : (-gr/L)*np.sin(y)
gr = 9.8
L = 2
x0 = 0
y0 = np.pi/6
z0 = 0
v0 = [x0,y0,z0]
h = 0.1
xn = 2
m = int((xn-x0)/h)
# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,v0,h,m)
xi = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]
# SALIDA
np.set_printoptions(precision=6)
print('x, y, z')
print(tabla)
plt.plot(xi,yi, label='y')
plt.plot(xi,zi, label='z')
plt.legend()
plt.grid()
plt.show()
-1.05625 U(1,1) + 0.028125 U(2,1) = -1.732050
0.028125 U(1,1) - 1.05625 U(2,1) = -1.732050
A = [[-1.05625 , 0.028125],
[ 0.028125, -1.05625 ]]
B = [-1.732050,-1.732050]
que resuelta con un método numérico:
[ 1.68, 1.68]
Por lo que la solución para una gráfica, con los índices de (fila,columna) como (t,x):
U = [[0, 1.732050, 1.732050, 0],
[0, 1.680000, 1,680000, 0]]
El error del procedimiento, tal como fué planteado es del orden de O(Δt) y O(Δx2), o error de truncamiento E = O(Δx2) + O(Δt). Δt debe ser menor que Δx en aproximadamente un orden de magnitud
Usando algoritmo en python.
Usando lo resuelto en clase y laboratorio, se comprueba la solución con el algoritmo, con hx y ht mas pequeños y más iteraciones:
# EDP parabólicas d2u/dx2 = K du/dt# método implícito# Referencia: Chapra 30.3 p.895 pdf.917# Rodriguez 10.2.5 p.417import numpy as np
import matplotlib.pyplot as plt
# INGRESO# Valores de frontera
Ta = 0
Tb = 0
# longitud en x
a = 0
b = 1
# Constante K
K = 16
# Tamaño de paso
dx = 0.1
dt = 0.01
# temperatura en barra
tempbarra = lambda x: 2*np.sin(np.pi*x)
# iteraciones
n = 100
# PROCEDIMIENTO# Valores de x
xi = np.arange(a,b+dx,dx)
m = len(xi)
# Resultados en tabla de u
u = np.zeros(shape=(m,n), dtype=float)
# valores iniciales de u
j=0
u[0,j] = Ta
for i inrange(1,m-1,1):
u[i,j] = tempbarra(xi[i])
u[m-1,j] = Tb
# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = -1 -2*lamb
R = lamb
vector = np.array([P,Q,R])
tvector = len(vector)
# Calcula U para cada tiempo + dt
j=1
whilenot(j>=n):
u[0,j] = Ta
u[m-1,j] = Tb
# Matriz de ecuaciones
tamano = m-2
A = np.zeros(shape=(tamano,tamano), dtype = float)
B = np.zeros(tamano, dtype = float)
for f inrange(0,tamano,1):
for c inrange(0,tvector,1):
c1 = f+c-1
if(c1>=0 and c1<tamano):
A[f,c1] = vector[c]
B[f] = -u[f+1,j-1]
B[0] = B[0]-P*u[0,j]
B[tamano-1] = B[tamano-1]-R*u[m-1,j]
# Resuelve sistema de ecuaciones
C = np.linalg.solve(A, B)
# copia resultados a u[i,j]for f inrange(0,tamano,1):
u[f+1,j] = C[f]
j=j+1 # siguiente iteración# SALIDAprint('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)
# Gráfica
salto = int(n/10)
if (salto == 0):
salto = 1
for j inrange(0,n,salto):
vector = u[:,j]
plt.plot(xi,vector)
plt.plot(xi,vector, '.m')
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()
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
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 clasedefitrapecio (datos,dt):
n=len(datos)
integral=0
for i inrange(0,n-1,1):
area=dt*(datos[i]+datos[i+1])/2
integral=integral + area
return(integral)
defisimpson13(f,h):
n = len(f)
integral = 0
for i inrange(0,n-2,2):
area = (h/3)*(f[i]+4*f[i+1]+f[i+2])
integral = integral + area
return(integral)
defisimpson38 (f,h):
n=len(f)
integral=0
for i inrange(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 inrange(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()
