para la evaluación numérica de 1 se usa un valor muy cercano desplazado con la tolerancia aplicada.
Igx≅83(0.3333)[0+3(0.2703)+3(0.3662)+0]=0.2387
d. Indique el resultado obtenido para el área requerida y la cota de error
Area = I_{fx} – I_{gx} = 1.3333 – 0.2387 = 1.0945
cota de error = O(0.03125) + O(0.00411) = 0.03536
e. Encuentre el valor del tamaño de paso si se requiere una cota de error de 0.00032
Si el factor de mayor error es de Simpson 1/3, se considera como primera aproximación que:
cota de error O(h5) = O(0.00032), h = (0.00032)(1/5) = 0.2
es decir el número de tramos es de al menos (b-a)/tramos = 0.2 , tramos = 5.
El número de tramos debe ser par en Simpson de 1/3, por lo que se toma el entero mayor tramos=6 y el tamaño de paso recomendado es al menos 1/6. EL error al aplicar 3 veces la formula es 3(O((1/6)5)) = 0.0003858.
Lo que podría indicar que es necesario al menos dos tramos adicionales con h=1/8 y error O(0,00012) que cumple con el requerimiento.
Se puede aplicar el mismo criterio para Simpson 3/8 y se combinan los errores para verificar que cumplen con el requerimiento.
# 2Eva_2023PAOI_T1 Material para medalla de academiaimport numpy as np
import matplotlib.pyplot as plt
defintegrasimpson13_fi(xi,fi,tolera = 1e-10):
''' sobre muestras de fi para cada xi
integral con método de Simpson 1/3
respuesta es np.nan para tramos desiguales,
no hay suficientes puntos.
'''
n = len(xi)
i = 0
suma = 0
whilenot(i>=(n-2)):
h = xi[i+1]-xi[i]
dh = abs(h - (xi[i+2]-xi[i+1]))
if dh<tolera:# tramos iguales
unS13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
suma = suma + unS13
else: # tramos desiguales
suma = 'tramos desiguales'
i = i + 2
if i<(n-1): # incompleto, faltan tramos por calcular
suma = 'tramos incompletos, faltan '
suma = suma ++str((n-1)-i)+' tramos'return(suma)
defintegrasimpson38_fi(xi,fi,tolera = 1e-10):
''' sobre muestras de fi para cada xi
integral con método de Simpson 3/8
respuesta es np.nan para tramos desiguales,
no hay suficientes puntos.
'''
n = len(xi)
i = 0
suma = 0
whilenot(i>=(n-3)):
h = xi[i+1]-xi[i]
h1 = (xi[i+2]-xi[i+1])
h2 = (xi[i+3]-xi[i+2])
dh = abs(h-h1)+abs(h-h2)
if dh<tolera:# tramos iguales
unS38 = fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]
unS38 = (3/8)*h*unS38
suma = suma + unS38
else: # tramos desiguales
suma = 'tramos desiguales'
i = i + 3
if (i+1)<n: # incompleto, tramos por calcular
suma = 'tramos incompletos, faltan '
suma = suma +str(n-(i+1))+' tramos'return(suma)
# INGRESO
fx = lambda x: 2-8*(0.5-x)**2
gx = lambda x: -(1-x)*np.log(1-x)
a = 0
b = 1-1e-4
muestras1 = 2+1
muestras2 = 3+1
# PROCEDIMIENTO
xi1 = np.linspace(a,b,muestras1)
xi2 = np.linspace(a,b,muestras2)
fi = fx(xi1)
gi = gx(xi2)
Ifx = integrasimpson13_fi(xi1,fi)
Igx = integrasimpson38_fi(xi2,gi)
Area = Ifx - Igx
# SALIDAprint('Ifx: ', Ifx)
print('Igx: ', Igx)
print('Area: ', Area)
plt.plot(xi1,fi,'ob',label='f(x)')
plt.plot(xi2,gi,'or', label='g(x)')
plt.grid()
plt.legend()
plt.xlabel('xi')
# curvas suave con mas muestras (no en evaluación)
xi = np.linspace(a,b,51)
fxi = fx(xi)
gxi = gx(xi)
plt.fill_between(xi,fxi,gxi,color='navajowhite')
plt.plot(xi,fxi,color='blue',linestyle='dotted')
plt.plot(xi,gxi,color='red',linestyle='dotted')
plt.show()
Tema 1 (30 puntos)
Una academia encarga a un joyero un modelo de medalla cuyo costo unitario se determina por el área descrita entre las funciones f(x) y g(x) presentadas.
Se considera que el grosor de la medalla es único e independiente de la forma de la medalla.
f(x)=2−8(21−x)20≤x<1g(x)=−(1−x)ln(1−x)
Para el desarrollo numérico, use diferentes métodos de Simpson para cada función.
a. Realice el planteamiento de las ecuaciones para el ejercicio.
b. Describa el criterio usado para determinar el número de tramos usado en cada caso.
c. Desarrolle las expresiones completas del ejercicio para cada función.
d. Indique el resultado obtenido para el área requerida y la cota de error.
e. Encuentre el valor del tamaño de paso si se requiere una cota de error de 0.00032
Nota: en Python ln(x) se escribe np.log(x).
Rúbrica: literal a (5 puntos), literal b (5 puntos), literal c (10 puntos), literal d (5 puntos), literal e (5 puntos)
Referencia: Star Trek https://intl.startrek.com/
¿A quien se le ocurrió crear la moneda? | Discovery en Español Youtube.com 8 nov 2016.
El desarrollo analítico inicia convirtiendo la ecuación diferencial parcial elíptica de la forma contínua a su representación discreta usando las expresiones de derivadas en diferencias divididas. El cambio de la expresión se realiza usando Sympy.
La EDP se escribe en formato Sympy en el bloque de ingreso. La ecuación EDP se compone del lado izquierdo (LHS) y lado derecho (RHS), indicando u como una función de las variables x,y.
las diferencias divididas se ingresan como un diccionario:
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
sym.diff(u,y,2): (U.subs(j,j-1)-2*U+U.subs(j,j+1))/(Dy**2),
sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}
las condiciones iniciales en los bordes, pueden ser funciones matemáticas, por lo que se usa el formato lambda. En el ejercicio básico presentado en la parte teórica, los bordes tienen temperaturas constantes, por lo que en ésta sección se escriben las condiciones como la constante mas cero por la variable independiente, facilitando la evaluación de un vector xi por ejemplo.
los demás parámetros del ejercicio se ingresan más adelante de forma semejante al ejercicio con Numpy.
Todas las expresiones se escriben al lado izquierdo (LHS) de la igualdad, la parte del lado derecho(RHS) se la prefiere mantener en cero. La expresión se organiza manteniendo el coeficiente en 1 para el termino de Dx de mayor orden. Se sustituyen las derivadas por diferencias divididas, para luego simplificar la expresión usando λ como referencia.
El resultado de las operaciones en la expresión usando el algoritmo es:
3. EDP Elípticas – Método iterativo con Sympy-Python
El desarrollo del método iterativo para una ecuación diferencial parcial elíptica usando Sympy, sigue a continuación del anterior.
Se añaden las instrucciones para iterar la expresión discreta para cada i,j dentro de la matriz U creada para el ejercicio. Se aplican las condiciones iniciales a los bordes, para luego proceder a iterar.
Como en las actividades del curso se requiere realizar al menos 3 iteraciones para las expresiones del algoritmo con papel y lápiz, solo se presentarán los resultados para la primera fila en j=0
Los resultados del algoritmo luego del resultado del algoritmo anterior se presentan como:
Se pueden sustituir los valores de las expresiones, evaluando cada u[i,j] y completando la matriz y para la gráfica. Esta parte queda como tarea.
Instrucciones en Python
Las instrucciones adicionales al algoritmo de EDP elíptica contínua a discreta empiezan con la creación de la matriz u_xy para los valores, los vectores xi,yj y una matriz u_mask que indica si se dispone de un valor calculado o conocido en ese punto (x,y) . Para el método iterativo, la matriz se rellena con el promedio de los valores máximos de las funciones dadas para los bordes de la placa.
Tarea: El algoritmo desarrolla el cálculo para j_k=1, solo la primera fila. Se deja como tarea desarrollar para toda la matriz.
# ITERAR para cada i,j dentro de U ------------# x[i] , y[j] valor en posición en cada eje
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)
# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u# llena u con valores en fronteras
u_xy[0,:] = fya(yj) # izquierda Ta
u_xy[n-1,:] = fyb(yj) # derecha Tb
u_xy[:,0] = fxc(xi) # inferior Tc
u_xy[:,m-1] = fxd(xi) # superior Td
u0 = np.copy(u_xy) # matriz u inicial# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:] = True# izquierda
u_mask[-1,:] = True# derecha
u_mask[:,0] = True# inferior
u_mask[:,-1] = True# superior# valor inicial de iteración dentro de u# promedio = (Ta+Tb+Tc+Td)/4
promedio = (np.max(fya(yj))+np.max(fyb(yj))+\
np.max(fxc(xi))+np.max(fxd(xi)))/4
u_xy[1:n-1,1:m-1] = promedio
u_mask[1:n-1,1:m-1] = True
u0 = np.copy(u_xy) # copia para revisióndefedp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
'''Sustituye en edp discreta los valores conocidos de U[i,j]
tomados desde u_xy, marcados con u_mask
u_mask indica si el valor se ha calculado con edp.
'''
LHS = discreta.lhs # lado izquierdo de ecuacion
RHS = discreta.rhs # lado derecho# sustituye U[i,j] con valores conocidos
A_diagonal = [] # lista de i,j para matriz de coeficientes A# Separa términos suma
term_suma = sym.Add.make_args(LHS)
for term_k in term_suma:
# busca U[i,j] y cambia por valor uxt[i,j]
cambiar = 0 ; cambiar_valor = 0 ; cambiar_factor = 0
# separa cada factor de término
factor_Mul = sym.Mul.make_args(term_k)
for factor_k in factor_Mul:
# busca U[i,j] en matriz uxt[i,j]if factor_k.is_Function:
[i_k,j_k] = factor_k.args
ifnot(np.isnan(u_xy[i_k,j_k])):
cambiar = u_mask[i_k,j_k]
cambiar_factor = factor_k
cambiar_valor = u_xy[i_k,j_k]
else:
A_diagonal.append([i_k,j_k,term_k/factor_k])
# encontró valor U[i,j],term_k va a la derecha de ecuaciónif cambiar:
LHS = LHS - term_k
term_ki = term_k.subs(cambiar_factor,cambiar_valor)
RHS = RHS - term_ki
discreta = sym.Eq(LHS,RHS)
B_diagonal = RHS
resultado = [discreta,A_diagonal,B_diagonal]
return (resultado)
# Desarrollo de iteraciones en cada nodo i,j, para la fila j_k# genera el valor en cada punto# Nota: Solo la primera fila, Tarea desarrollar para la matrizprint('\n iterar i,j:')
j_k = 1
for i_k inrange(1,n-1,1):
print('j: ',j_k,' ; ','i: ',i_k)
discreta_ij = discreta.subs({i:i_k,j:j_k,
x:xi[i_k],y:yj[j_k]})
sym.pprint(discreta_ij)
RHS = discreta_ij.rhs
discreta_ij = sym.Eq(RHS,0)
# usa valores de frontera segun u_mask con True
discreta_k = edp_sustituyeValorU(discreta_ij,
xi,yj,u_xy,u_mask)
u_xy[i_k,j_k] = sym.Float(-discreta_k[2])
print(buscar,'=',u_xy[i_k,j_k])
# SALIDA
np.set_printoptions(precision=verdigitos)
print('\nMétodo iterativo EDP Elíptica')
print('continuar el desarrollo con: ')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k inrange(2,-1,-1):
print(j_k, (u_xy[:,j_k]))
4. EDP Elípticas – Método Implícito con Sympy-Python
Desarrollo analítico del método implicito para una ecuación diferencial parcial elíptica usando Sympy. El algoritmo reutiliza el algoritmo para la EDP Elíptica de contínua a discreta, la creación de la matriz de valores u_xy, y la función edp_sustituyeValorU() para buscar los valores conocidos de la u(x,y).
El algoritmo usa la ecuación discreta para en cada iteración i,j reemplazar los valores de U conocidos. Los valores de U se escriben en una matriz u_xy, para diferenciar si el valor existe se usa la matriz de estados u_mask.
Con las ecuaciones de cada iteración se llena la matriz A de coeficientes y el vector B de las constantes. Al resolver el sistema de ecuaciones se obtienen todos los valores de la matriz U, completando el ejercicio. Se desarrola la solución al sistema de ecuaciones usando Sympy como alternativa a usar numpy con np.linalg.solve(A.B) que se encuentra como comentario entre las instrucciones.
# Ecuaciones Diferenciales Parciales Elipticas# EDP Elípticas contínua a discreta con Sympyimport numpy as np
import sympy as sym
# funciones continuas y variables simbólicas usadas
y = sym.Symbol('y',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,y)
f = sym.Function('f')(x,y) # funcion complemento# funciones discretas y variables simbólicas usadas
i = sym.Symbol('i',integer=True,positive=True)
j = sym.Symbol('j',integer=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L = sym.Symbol('L',real=True)
U = sym.Function('U')(i,j)
# INGRESO
fxy = 0*x+0*y # f(x,y) = 0 , ecuacion de Poisson# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) + sym.diff(u,y,2)
RHS = fxy
edp = sym.Eq(LHS-RHS,0)
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
sym.diff(u,y,2): (U.subs(j,j-1)-2*U+U.subs(j,j+1))/(Dy**2),
sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}
# Condiciones iniciales en los bordes
fya = lambda y: 60 +0*y # izquierda
fyb = lambda y: 25 +0*y # derecha
fxc = lambda x: 50 +0*x # inferior, función inicial
fxd = lambda x: 70 +0*x # superior, función inicial # dimensiones de la placa
x0 = 0 # longitud en x
xn = 2
y0 = 0 # longitud en y
yn = 1.5
# discretiza, supone dx=dy
dx = 0.25 # Tamaño de paso
dy = dx
iteramax = 100
tolera = 0.0001
verdigitos = 2 # para mostrar en tabla de resultados# PROCEDIMIENTOdefedp_discreta(edp,dif_dividida,x,y,u):
''' EDP contínua a discreta, usa diferencias divididas
proporcionadas en parámetros, indica las variables x,y
con función u de (x,y)
'''
resultado={}
# expresión todo a la izquierda LHS (Left Hand side)
LHS = edp.lhs
RHS = edp.rhs
ifnot(edp.rhs==0):
LHS = LHS-RHS
RHS = 0
edp = sym.Eq(LHS,RHS)
# para orden de derivada por x, y
edp_x = edp.subs(x,0)
edp_y = edp.subs(y,0)
ordenDx = sym.ode_order(edp_x,u)
ordenDy = sym.ode_order(edp_y,u)
resultado['ordenDx'] = ordenDx
resultado['ordenDy'] = ordenDy
# convierte coeficiente d2u/dx2 a 1
coeff_x = edp_coef_Dx(edp,x,ordenDx)
ifnot(coeff_x==1):
LHS = LHS/coeff_x
RHS = RHS/coeff_x
edp = sym.Eq(LHS,RHS)
K_ = edp_coef_Dx(edp,y,ordenDy)
resultado['edp=0'] = edp
resultado['K_'] = K_
# edp discreta
discreta = edp.lhs
for derivada in dif_dividida:
discreta = discreta.subs(derivada,dif_dividida[derivada])
resultado['discreta=0'] = discreta
return (resultado)
defedp_coef_Dx(edp,x,ordenx):
''' Extrae el coeficiente de la derivada Dx de ordenx
'''
coeff_x = 1.0
# separa cada término de suma
term_suma = sym.Add.make_args(edp.lhs)
for term_k in term_suma:
if term_k.is_Mul: # término de mas de un factor
factor_Mul = sym.Mul.make_args(term_k)
coeff_temp = 1; coeff_usar=0
# separa cada factor de término for factor_k in factor_Mul:
ifnot(factor_k.is_Derivative):
coeff_temp = coeff_temp*factor_k
else: # factor con derivada de ordenx
partes = factor_k.args
if partes[1]==(x,ordenx):
coeff_usar = 1
if coeff_usar==1:
coeff_x = coeff_x*coeff_temp
return(coeff_x)
defedp_simplificaLamba(resultado,dx,dy):
'''simplifica ecuacion con valores de lambda, dx y dy
entregando la edp discreta simplificada
'''
discreta = resultado['discreta=0']
ordenDy = resultado['ordenDy']
ordenDx = resultado['ordenDx']
K_ = resultado['K_']
lamb = (Dy**ordenDy)/(Dx**ordenDx)
if ordenDy==1 and ordenDx==2:
lamb = lamb/K_
resultado['Lambda_L'] = lamb
# valor de Lambda en ecuacion edp
L_k = lamb.subs([(Dx,dx),(Dy,dy),(1.0,1)])
resultado['Lambda L_k'] = L_k
# simplifica con lambda L
discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
resultado['(discreta=0)*Dy**ordeny'] = discreta_L
discreta_L = edp_sustituye_L(resultado)
discreta_L = discreta_L.subs(lamb,L)
discreta_L = sym.collect(discreta_L,U)
discreta_L = sym.Eq(discreta_L,0)
resultado['discreta_L=0'] = discreta_L
# sustituye constantes en ecuación a iterar
discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),
(L,L_k),(1.0,1)])
resultado['discreta'] = discreta_L
return(resultado)
defedp_sustituye_L(resultado):
''' sustituye lambda con Dy**ordeny/Dx**x/K_
por L, al simplificar Lambda
'''
discreta = resultado['(discreta=0)*Dy**ordeny']
ordenDy = resultado['ordenDy']
ordenDx = resultado['ordenDx']
discreta_L = 0
# separa cada término de suma
term_suma = sym.Add.make_args(discreta)
for term_k in term_suma:
# busca partes de L y cambia por valor L
cambiar = 0 # por orden de derivadaif term_k.has(Dx) and term_k.has(Dy):
partes = term_k.args
ordeny=1
ordenx=1
for unaparte in partes:
if unaparte.has(Dy):
if unaparte.is_Pow:
partey = unaparte.args
ordeny = partey[1]
if unaparte.has(Dx):
if unaparte.is_Pow:
partey = unaparte.args
ordenx = partey[1]
if (ordeny<=ordenDy and ordenx<=-ordenDx):
cambiar=1
if cambiar:
term_k = term_k*L/resultado['Lambda_L']
discreta_L = discreta_L + term_k
# simplifica unos con decimal a entero
discreta_L = discreta_L.subs(1.0,1)
return(discreta_L)
defedp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
'''Sustituye en edp discreta los valores conocidos de U[i,j]
tomados desde u_xy, marcados con u_mask
u_mask indica si el valor se ha calculado con edp.
'''
LHS = discreta.lhs # lado izquierdo de ecuacion
RHS = discreta.rhs # lado derecho# sustituye U[i,j] con valores conocidos
A_diagonal = [] # lista de i,j para matriz de coeficientes A# Separa términos suma
term_suma = sym.Add.make_args(LHS)
for term_k in term_suma:
# busca U[i,j] y cambia por valor uxt[i,j]
cambiar = 0 ; cambiar_valor = 0 ; cambiar_factor = 0
# separa cada factor de término
factor_Mul = sym.Mul.make_args(term_k)
for factor_k in factor_Mul:
# busca U[i,j] en matriz uxt[i,j]if factor_k.is_Function:
[i_k,j_k] = factor_k.args
ifnot(np.isnan(u_xy[i_k,j_k])):
cambiar = u_mask[i_k,j_k]
cambiar_factor = factor_k
cambiar_valor = u_xy[i_k,j_k]
else:
A_diagonal.append([i_k,j_k,term_k/factor_k])
# encontró valor U[i,j],term_k va a la derecha de ecuaciónif cambiar:
LHS = LHS - term_k
term_ki = term_k.subs(cambiar_factor,cambiar_valor)
RHS = RHS - term_ki
discreta = sym.Eq(LHS,RHS)
B_diagonal = RHS
resultado = [discreta,A_diagonal,B_diagonal]
return (resultado)
# PROCEDIMIENTO# transforma edp continua a discreta
resultado = edp_discreta(edp,dif_dividida,x,y,u)
resultado = edp_simplificaLamba(resultado,dx,dy)
discreta = resultado['discreta']
# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Elíptica contínua a discreta')
for entrada in resultado:
tipo = type(resultado[entrada])
if tipo in algun_numero or entrada in algun_numero:
print('',entrada,':',resultado[entrada])
else:
print('\n',entrada,':')
sym.pprint(resultado[entrada])
# ITERAR para cada i,j dentro de U ------------# x[i] , y[j] valor en posición en cada eje
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)
# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u# llena u con valores en fronteras
u_xy[0,:] = fya(yj) # izquierda Ta
u_xy[n-1,:] = fyb(yj) # derecha Tb
u_xy[:,0] = fxc(xi) # inferior Tc
u_xy[:,m-1] = fxd(xi) # superior Td
u0 = np.copy(u_xy) # matriz u inicial# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:] = True# izquierda
u_mask[-1,:] = True# derecha
u_mask[:,0] = True# inferior
u_mask[:,-1] = True# superior# Método implícito para EDP Elíptica# ITERAR para plantear las ecuaciones en [i,j]
resultado = {}
eq_itera = [] ; tamano = (n-2)*(m-2)
A = np.zeros(shape=(tamano,tamano),dtype=float)
B = np.zeros(tamano,dtype=float)
for j_k inrange(1,m-1,1): # no usar valores en bordesfor i_k inrange(1,n-1,1):
eq_conteo = (j_k - 1)*(n-2)+(i_k-1)
discreta_ij = discreta.subs({i:i_k,j:j_k,
x:xi[i_k],y:yj[j_k]})
resultado[eq_conteo]= {'j':j_k, 'i':i_k,
'discreta_ij': discreta_ij}
# usa valores de frontera segun u_mask con True
discreta_k = edp_sustituyeValorU(discreta_ij,
xi,yj,u_xy,u_mask)
discreta_ij = discreta_k[0]
A_diagonal = discreta_k[1] # lista de (i,j,coeficiente)
B_diagonal = discreta_k[2]
resultado[eq_conteo]['discreta_k'] = discreta_k[0]
# añade ecuacion a resolver
eq_itera.append(discreta_ij)
# Aplica coeficientes de ecuacion en A y B:# A_diagonal tiene lista de (i,j,coeficiente) for uncoeff in A_diagonal:
columna = (uncoeff[1]-1)*(n-2)+(uncoeff[0]-1)
fila = (j_k - 1)*(n-2)+(i_k-1)
A[fila,columna] = uncoeff[2]
B[eq_conteo] = float(B_diagonal) # discreta_ij.rhs
resultado['A'] = np.copy(A)
resultado['B'] = np.copy(B)
# resuelve el sistema de ecuaciones en eq_itera en Sympy
X_k = sym.solve(eq_itera)[0]
# actualiza uxt[i,j] , u_mask segun X_k en Sympyfor nodo_Uij in X_k:
[i_k,j_k] = nodo_Uij.args
u_xy[i_k,j_k] = X_k[nodo_Uij]
u_mask[i_k,j_k] = True# resuelve el sistema A.X=B en Numpy#X = np.linalg.solve(A,B)# tarea: llevar valores X a u_xy# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('\nMétodo Implícito - EDP Elíptica')
for entrada in resultado:
tipo = type(resultado[entrada])
if tipo in algun_numero or entrada in algun_numero:
print('',entrada,':',resultado[entrada])
elif (tipo==dict):
print('j:',resultado[entrada]['j'],'; ''i:',resultado[entrada]['i'],'; ''ecuacion:',entrada)
sym.pprint(resultado[entrada]['discreta_ij'])
if resultado[entrada]['discreta_k'].rhs!=0:
sym.pprint(resultado[entrada]['discreta_k'])
else:
print('\n',entrada,': \n',resultado[entrada])
print('Resultados para U(x,y)')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k inrange(m-1,-1,-1):
print(j_k, (u_xy[:,j_k]))
El desarrollo analítico comienza en escribir la ecuación diferencial parcial parabólica de la forma contínua a la representación discreta, siguiendo los criterios indicados en la parte teórica usando derivadas expresadas en diferencias divididas.
La EDP se escribe en formato Sympy en el bloque de ingreso en dos partes: lado izquierdo (LHS) de la ecuación y lado derecho (RHS), definiendo u como una función de variables x,y.
Las condiciones iniciales en los extremos a y b, y la temperatura en la barra, pueden ser funciones matemáticas, por lo que se propone usar el formato lambda de las variables x o y. En el ejercicio básico presentado en la parte teórica, la barra tiene temperatura constante, por lo que escribe como una constante mas cero por la variable independiente, esta forma facilita la evaluación de un vector xi por ejemplo, en caso de disponer una expresión diferente.
los demás parámetros del ejercicio se ingresan más adelante de forma semejante al ejercicio con Numpy.
Todas las expresiones se escriben al lado izquierdo (LHS) de la igualdad, la parte del lado derecho(RHS) se la prefiere mantener en cero. La expresión se organiza manteniendo el coeficiente en 1 para el termino de Dx de mayor orden. Se sustituyen las derivadas por diferencias divididas, para luego simplificar la expresión usando λ como referencia.
Al desarrollar el método explícito se requiere indicar el nodo u(i,j+1) a buscar
buscar = U.subs(j,j+1) # método explícito
El valor de K se puede determinar en la expresión al reordenar para que el término de Dx de mayor grado sea la unidad. El valor de K será el del término Dy con signo negativo, al encontrarse en el lado izquierdo (LHS) de la expresión, en la propuesta teórica se encontraba como positivo del lado derecho.
El término de mayor grado se puede obtener al revisar la expresión por cada término de la suma, buscando el factor que contiene Dx o Dy en cada caso. Por ser un poco ordenado se añade la función edp_coef_Dx().
El resultado del algoritmo para el ejercicio propuesto es:
Para revisar los resultados de algunos pasos para obtener la expresión discreta de EDP, se usa un diccionario a ser usado en el bloque de salida. Para las expresiones de ecuación se usa sym.pprint().
# Ecuaciones Diferenciales Parciales Parabólicas# EDP Parabólicas contínua a discreta con Sympyimport numpy as np
import sympy as sym
# funciones continuas y variables simbólicas usadas
t = sym.Symbol('t',real=True)
y = sym.Symbol('y',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,y)
K = sym.Symbol('K',real=True)
f = sym.Function('f')(x,y) # funcion complemento# funciones discretas y variables simbólicas usadas
i = sym.Symbol('i',integer=True,positive=True)
j = sym.Symbol('j',integer=True,positive=True)
Dt = sym.Symbol('Dt',real=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L = sym.Symbol('L',real=True)
U = sym.Function('U')(i,j)
# INGRESO
fxy = 0*x+0*y # f(x,y) = 0 , ecuacion complementaria# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2)
RHS = 4*sym.diff(u,y,1) + fxy
edp = sym.Eq(LHS-RHS,0)
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
sym.diff(u,y,1): (U.subs(j,j+1)-U)/Dy}
buscar = U.subs(j,j+1) # método explícito# Valores de frontera
fya = lambda y: 60 +0*y # izquierda
fyb = lambda y: 40 +0*y # derecha
fxc = lambda x: 25 +0*x # inferior, función inicial# dimensiones de la barra
a = 0 # longitud en x
b = 1
y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones# discretiza, supone dx=dy
x_tramos = 10
dx = (b-a)/x_tramos
dy = dx/10
verdigitos = 2 # para mostrar en tabla de resultados
n_iteraciones = 4 # iteraciones en tiempo# PROCEDIMIENTOdefedp_coef_Dx(edp,x,ordenx):
''' Extrae el coeficiente de la derivada Dx de ordenx
'''
coeff_x = 1.0
# separa cada término de suma
term_suma = sym.Add.make_args(edp.lhs)
for term_k in term_suma:
if term_k.is_Mul: # mas de un factor
factor_Mul = sym.Mul.make_args(term_k)
coeff_temp = 1; coeff_usar=0
# separa cada factor de término for factor_k in factor_Mul:
ifnot(factor_k.is_Derivative):
coeff_temp = coeff_temp*factor_k
else: # factor con derivada de ordenx
partes = factor_k.args
if partes[1]==(x,ordenx):
coeff_usar = 1
if coeff_usar==1:
coeff_x = coeff_x*coeff_temp
return(coeff_x)
resultado={}
# expresión todo a la izquierda LHS (Left Hand side)ifnot(edp.rhs==0):
LHS = EDP.lhs
RHS = EDP.rhs
edp = sym.Eq(LHS,RHS)
# para orden de derivada por x, y
edp_x = edp.subs(x,0)
edp_y = edp.subs(y,0)
ordenDx = sym.ode_order(edp_x,u)
ordenDy = sym.ode_order(edp_y,u)
resultado['ordenDx'] = ordenDx
resultado['ordenDy'] = ordenDy
# convierte coeficiente d2u/dx2 a 1
coeff_x = edp_coef_Dx(edp,x,ordenDx)
ifnot(coeff_x==1):
LHS = LHS/coeff_x
RHS = RHS/coeff_x
edp = sym.Eq(LHS-RHS,0)
K_ = -edp_coef_Dx(edp,y,ordenDy)
resultado['edp=0'] = edp
resultado['K_'] = float(K_)
lamb = (Dy**ordenDy)/(Dx**ordenDx)
if ordenDy==1 and ordenDx==2:
lamb = lamb/K_
resultado['Lambda_L'] = lamb
# valor de Lambda en ecuacion edp
L_k = lamb.subs([(Dx,dx),(Dy,dy),(1.0,1)])
resultado['Lambda L_k'] = float(L_k)
discreta = edp.lhs # EDP discretafor derivada in dif_dividida:
discreta = discreta.subs(derivada,dif_dividida[derivada])
# simplifica con Lambda L
discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
resultado['(discreta=0)*Dy**ordeny'] = discreta_L
discreta_L = discreta_L.subs(lamb,L)
discreta_L = sym.collect(discreta_L,U)
discreta_L = sym.Eq(discreta_L,0)
resultado['discreta_L=0'] = discreta_L
# sustituye constantes en ecuación a iterar
discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),
(L,L_k),(1.0+0*x,1)])
discreta_L = discreta_L.subs(1.0+0*x,1)
resultado['discreta'] = discreta_L
# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Parabólica contínua a discreta')
for entrada in resultado:
tipo = type(resultado[entrada])
if tipo in algun_numero or entrada in algun_numero:
print('',entrada,':',resultado[entrada])
else:
print('\n',entrada,':')
sym.pprint(resultado[entrada])
El desarrollo del método explícito para una ecuación diferencial parcial parabólica usando Sympy, requiere añadir instrucciones para procesar la expresión obtenida en el algoritmo anterior con los valores de i,j.
Se añaden las instrucciones para iterar la expresión discreta para cada i,j dentro de la matriz U creada para el ejercicio. Se aplican las condiciones iniciales a los bordes de la matriz, para luego proceder a iterar.
Como en las actividades del curso se requiere realizar al menos 3 iteraciones para las expresiones del algoritmo con papel y lápiz, solo se presentarán los resultados de las expresiones para la primera fila en j=0.
El algoritmo para la EDP del ejemplo muestra el siguiente resultado:
Para el ejercicio presentado en el numeral 1, se desarrolla el método implícito.
Las diferencias finitas centradas y hacia atras se definen en el diccionario dif_dividida.
Para desarrollar la expresión discreta, se incorpora la función edp_sustituye_L(), que busca el en cada término de la suma los componentes de λ y los sustituye por la variable L.
Para cada valor de j, se crean las ecuaciones desde la forma discreta de la EDP, moviendo a la derecha los valores conocidos para generar un sistema de ecuaciones A.X=B. Se resuelve el sistema de ecuaciones y se actualiza la matriz u[i,j] de resultados.
Por la extensión de los pasos a mostrar, para las iteraciones en cada fila, solo se presentan para la primera fila. En caso de requerir observar más filas, puede cambiar el parámetro en el bloque de salida, en el condicional:
Para el método implícito se añade la sección donde se reemplazan los valores en los bordes y en la barra en la matriz U(x,y) en el algoritmo identificada como u_xy . Para diferenciar si el valor existe se usa la matriz de estados u_mask con estados True/False.
También se incorpora la función edp_sustituyeValorU() que para cada término suma de la ecuación busca en la matriz u_xy si el valor U(i,j) existe y lo cambia. Luego pasa el valor al lado derecho de la ecuación par formar el sistema de ecuaciones.
Con los valores existentes antes de la iteración, se obtienen las ecuaciones por cada punto i,j de una fila para el sistema de ecuaciones en la forma A.x=B. El sistema se resuelve con las instrucciones de Numpy. El resultado se actualiza en la matriz u_xy y la matriz u_mask para indicar que el valor ya está disponible.
Compare los resultados con el método numérico y mida los tiempos de ejecución de cada algoritmo para escribir sus observaciones.
# Ecuaciones Diferenciales Parciales Parabólicas# EDP Parabólicas contínua a discreta con Sympyimport numpy as np
import sympy as sym
# funciones continuas y variables simbólicas usadas
t = sym.Symbol('t',real=True)
y = sym.Symbol('y',real=True)
x = sym.Symbol('x',real=True)
u = sym.Function('u')(x,y)
K = sym.Symbol('K',real=True)
f = sym.Function('f')(x,y) # funcion complemento# funciones discretas y variables simbólicas usadas
i = sym.Symbol('i',integer=True,positive=True)
j = sym.Symbol('j',integer=True,positive=True)
Dt = sym.Symbol('Dt',real=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L = sym.Symbol('L',real=True)
U = sym.Function('U')(i,j)
# INGRESO
fxy = 0*x+0*y # f(x,y) = 0 , ecuacion complementaria# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2)
RHS = 4*sym.diff(u,y,1) + fxy
edp = sym.Eq(LHS-RHS,0)
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}
# Valores de frontera
fya = lambda y: 60 +0*y # izquierda
fyb = lambda y: 40 +0*y # derecha
fxc = lambda x: 25 +0*x # inferior, función inicial# dimensiones de la barra
a = 0 # longitud en x
b = 1
y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones# discretiza
x_tramos = 10
dx = (b-a)/x_tramos
dy = dx/10
verdigitos = 2 # para mostrar en tabla de resultados
n_iteraciones = 4 # iteraciones en tiempo# PROCEDIMIENTOdefedp_coef_Dx(edp,x,ordenx):
''' Extrae el coeficiente de la derivada Dx de ordenx
'''
coeff_x = 1.0
# separa cada término de suma
term_suma = sym.Add.make_args(edp.lhs)
for term_k in term_suma:
if term_k.is_Mul: # mas de un factor
factor_Mul = sym.Mul.make_args(term_k)
coeff_temp = 1; coeff_usar=0
# separa cada factor de término for factor_k in factor_Mul:
ifnot(factor_k.is_Derivative):
coeff_temp = coeff_temp*factor_k
else: # factor con derivada de ordenx
partes = factor_k.args
if partes[1]==(x,ordenx):
coeff_usar = 1
if coeff_usar==1:
coeff_x = coeff_x*coeff_temp
return(coeff_x)
defedp_sustituye_L(resultado):
''' sustituye lambda con Dy**ordeny/Dx**x/K_
por L, al simplificar Lambda
'''
discreta = resultado['(discreta=0)*Dy**ordeny']
ordenDy = resultado['ordenDy']
ordenDx = resultado['ordenDx']
discreta_L = 0
# separa cada término de suma
term_suma = sym.Add.make_args(discreta)
for term_k in term_suma:
# busca partes de L y cambia por valor L
cambiar = 0 # por orden de derivadaif term_k.has(Dx) and term_k.has(Dy):
partes = term_k.args
ordeny=1
ordenx=1
for unaparte in partes:
if unaparte.has(Dy):
if unaparte.is_Pow:
partey = unaparte.args
ordeny = partey[1]
if unaparte.has(Dx):
if unaparte.is_Pow:
partey = unaparte.args
ordenx = partey[1]
if (ordeny<=ordenDy and ordenx<=-ordenDx):
cambiar=1
if cambiar:
term_k = term_k*L/resultado['Lambda_L']
discreta_L = discreta_L + term_k
# simplifica unos con decimal a entero
discreta_L = discreta_L.subs(1.0,1)
return(discreta_L)
resultado={}
# expresión todo a la izquierda LHS (Left Hand side)ifnot(edp.rhs==0):
LHS = EDP.lhs
RHS = EDP.rhs
edp = sym.Eq(LHS,RHS)
# para orden de derivada por x, y
edp_x = edp.subs(x,0)
edp_y = edp.subs(y,0)
ordenDx = sym.ode_order(edp_x,u)
ordenDy = sym.ode_order(edp_y,u)
resultado['ordenDx'] = ordenDx
resultado['ordenDy'] = ordenDy
# convierte coeficiente d2u/dx2 a 1
coeff_x = edp_coef_Dx(edp,x,ordenDx)
ifnot(coeff_x==1):
LHS = LHS/coeff_x
RHS = RHS/coeff_x
edp = sym.Eq(LHS-RHS,0)
K_ = -edp_coef_Dx(edp,y,ordenDy)
resultado['edp=0'] = edp
resultado['K_'] = float(K_)
lamb = (Dy**ordenDy)/(Dx**ordenDx)
if ordenDy==1 and ordenDx==2:
lamb = lamb/K_
resultado['Lambda_L'] = lamb
# valor de Lambda en ecuacion edp
L_k = lamb.subs([(Dx,dx),(Dy,dy),(1.0,1)])
resultado['Lambda L_k'] = float(L_k)
discreta = edp.lhs # EDP discretafor derivada in dif_dividida:
discreta = discreta.subs(derivada,dif_dividida[derivada])
# simplifica con Lambda L
discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
resultado['(discreta=0)*Dy**ordeny'] = discreta_L
discreta_L = edp_sustituye_L(resultado)
discreta_L = sym.collect(discreta_L,U)
discreta_L = sym.Eq(discreta_L,0)
resultado['discreta_L=0'] = discreta_L
# sustituye constantes en ecuación a iterar
discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),
(L,L_k),(1.0+0*x,1)])
resultado['discreta'] = discreta_L
# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Parabólica contínua a discreta')
for entrada in resultado:
tipo = type(resultado[entrada])
if tipo in algun_numero or entrada in algun_numero:
print('',entrada,':',resultado[entrada])
else:
print('\n',entrada,':')
sym.pprint(resultado[entrada])
# Método implícito EDP Parabólica ------------# ITERAR para cada i,j dentro de U # x[i] , y[j] valor en posición en cada eje
xi = np.arange(a,b+dx,dx)
yj = np.arange(y0,n_iteraciones*dy+dy,dy)
n = len(xi)
m = len(yj)
# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u# llena u con valores en fronteras
u_xy[:,0] = fxc(xi) # inferior Tc
u_xy[0,:] = fya(yj) # izquierda Ta
u_xy[n-1,:] = fyb(yj) # derecha Tb
u0 = np.copy(u_xy) # matriz u inicial# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:] = True# izquierda Ta
u_mask[-1,:] = True# derecha Tb
u_mask[:,0] = True# inferior Tc# Método implícito para EDP Parabólica# a usar en iteraciones
discreta = resultado['discreta']
defedp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
'''Sustituye en edp discreta los valores conocidos de U[i,j]
tomados desde u_xy, marcados con u_mask
u_mask indica si el valor se ha calculado con edp.
'''
LHS = discreta.lhs # lado izquierdo de ecuacion
RHS = discreta.rhs # lado derecho# sustituye U[i,j] con valores conocidos
A_diagonal = [] # lista de i,j para matriz de coeficientes A# Separa términos suma
term_suma = sym.Add.make_args(LHS)
for term_k in term_suma:
# busca U[i,j] y cambia por valor uxt[i,j]
cambiar = 0 ; cambiar_valor = 0 ; cambiar_factor = 0
# separa cada factor de término
factor_Mul = sym.Mul.make_args(term_k)
for factor_k in factor_Mul:
# busca U[i,j] en matriz uxt[i,j]if factor_k.is_Function:
[i_k,j_k] = factor_k.args
ifnot(np.isnan(u_xy[i_k,j_k])):
cambiar = u_mask[i_k,j_k]
cambiar_factor = factor_k
cambiar_valor = u_xy[i_k,j_k]
else:
A_diagonal.append([i_k,j_k,term_k/factor_k])
# encontró valor U[i,j],term_k va a la derecha de ecuaciónif cambiar:
LHS = LHS - term_k
term_ki = term_k.subs(cambiar_factor,cambiar_valor)
RHS = RHS - term_ki
discreta = sym.Eq(LHS,RHS)
B_diagonal = RHS
resultado = [discreta,A_diagonal,B_diagonal]
return (resultado)
# ITERAR para plantear las ecuaciones en [i,j]
resultado = {}
eq_itera = [] ; tamano = (n-2)
A = np.zeros(shape=(tamano,tamano),dtype=float)
B = np.zeros(tamano,dtype=float)
for j_k inrange(1,m,1): # no usar valores en bordes
resultado[j_k] ={}
for i_k inrange(1,n-1,1):
eq_conteo = i_k-1
discreta_ij = discreta.subs({i:i_k,j:j_k,
x:xi[i_k],y:yj[j_k]})
resultado[j_k][eq_conteo]= {'j':j_k, 'i':i_k,
'discreta_ij': discreta_ij}
# usa valores de frontera segun u_mask con True
discreta_k = edp_sustituyeValorU(discreta_ij,
xi,yj,u_xy,u_mask)
discreta_ij = discreta_k[0]
A_diagonal = discreta_k[1] # lista de (i,j,coeficiente)
B_diagonal = discreta_k[2] # constante para B
resultado[j_k][eq_conteo]['discreta_k'] = discreta_k[0]
# añade ecuacion a resolver
eq_itera.append(discreta_ij)
# Aplica coeficientes de ecuacion en A y B:# A_diagonal tiene lista de (i,j,coeficiente) for uncoeff in A_diagonal:
columna = uncoeff[0]-1
fila = eq_conteo
A[fila,columna] = uncoeff[2]
B[eq_conteo] = float(B_diagonal) # discreta_ij.rhs# resuelve el sistema A.X=B en Numpy
X_k = np.linalg.solve(A,B)
# actualiza uxt[i,j] , u_mask segun X_k
u_xy[1:n-1,j_k] = X_k
u_mask[1:n-1,j_k] = True# almacena resultados
resultado[j_k]['A'] = np.copy(A)
resultado[j_k]['B'] = np.copy(B)
resultado[j_k]['X'] = np.copy(X_k)
# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('\nMétodo implícito EDP Parabólica')
for entrada in resultado:
tipo = type(resultado[entrada])
if tipo in algun_numero or entrada in algun_numero:
print('',entrada,':',resultado[entrada])
elif (tipo==dict):
if entrada<2:
eq_conteo = resultado[entrada].keys()
for una_eq in eq_conteo:
if una_eq=='A'or una_eq=='B'or una_eq=='X':
print(una_eq,':')
print(resultado[entrada][una_eq])
else:
print('j:',resultado[entrada][una_eq]['j'],'; ''i:',resultado[entrada][una_eq]['i'])
sym.pprint(resultado[entrada][una_eq]['discreta_ij'])
sym.pprint(resultado[entrada][una_eq]['discreta_k'])
else:
print('\n',entrada,': \n',resultado[entrada])
print('Resultados para U(x,y)')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k inrange(m-1,-1,-1):
print(j_k, (u_xy[:,j_k]))
Solo para fines didácticos, y como complemento para los ejercicios presentados en la unidad para Ecuaciones Diferenciales Parciales, se presentan las instrucciones para las animaciones usadas en la presentación de los conceptos y ejercicios. Los algoritmos para animación NO son necesarios para realizar los ejercicios, que requieren una parte analítica con al menos tres iteraciones en papel y lápiz. Se lo adjunta como una herramienta didáctica de asistencia para las clases.
# EDP parabólicas d2u/dx2 = K du/dt# método explícito,usando diferencias divididas# Referencia: Chapra 30.2 p.888 pdf.912, Rodriguez 10.2 p.406import numpy as np
# INGRESO# Valores de frontera
Ta = 60
Tb = 40
#T0 = 25 # estado inicial de barra
fx = lambda x: 25 + 0*x # función inicial para T0
a = 0 # longitud en x
b = 1
K = 4 # Constante K
dx = 0.1 # Tamaño de paso
dt = dx/10
n = 100 # iteraciones en tiempo
verdigitos = 2 # decimales a mostrar en tabla de resultados# PROCEDIMIENTO# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
fi = fx(xi)
m = len(xi)
ultimox = m-1
ultimot = n-1
# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n+1), dtype=float)
# valores iniciales de u[:,j]
j = 0
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb
# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb
# Calcula U para cada tiempo + dt
tj = np.arange(0,(n+1)*dt,dt)
j = 0
whilenot(j>ultimot): # igual con lazo forfor i inrange(1,ultimox,1):
u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
j=j+1
# SALIDA
np.set_printoptions(precision=verdigitos)
print('Método explícito EDP Parabólica')
print('lambda: ',np.around(lamb,verdigitos))
print('x:',xi)
print('t:',tj[0:len(xi)],'...',tj[-1])
print('Tabla de resultados en malla EDP Parabólica')
print('j, U[:,:',ultimox,'], primeras iteraciones')
for j inrange(ultimox,-1,-1):
print(j,u[:,j])
# GRAFICA ------------import matplotlib.pyplot as plt
tramos = 10
salto = int(n/tramos) # evita muchas líneasif (salto == 0):
salto = 1
for j inrange(0,n,salto):
vector = u[:,j]
plt.plot(xi,vector)
plt.plot(xi,vector, '.',color='red')
plt.xlabel('x[i]')
plt.ylabel('u[j]')
plt.title('Solución EDP parabólica')
#plt.show()
1.2 Gráficas EDP Parabólica, método explícito en 3D
Esta sección es la base para la generación de las animaciones siguientes, por lo que debe ser incluida antes de las dos siguientes animaciones 2D o 3D. Al menos hasta la parte donde se realiza la transpuesta de U.
# GRAFICA en 3D ------
tj = np.arange(0,n*dt,dt)
tk = np.zeros(tramos,dtype=float)
# Extrae parte de la matriz U,acorde a los tramos
U = np.zeros(shape=(m,tramos),dtype=float)
for k inrange(0,tramos,1):
U[:,k] = u[:,k*salto]
tk[k] = tj[k*salto]
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi,tk)
U = np.transpose(U)
fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U, color ='blue')
graf_3D.plot(xi[0],tk[1],U[1,0],'o',color ='orange')
graf_3D.plot(xi[1],tk[1],U[1,1],'o',color ='green')
graf_3D.plot(xi[2],tk[1],U[1,2],'o',color ='green')
graf_3D.plot(xi[1],tk[2],U[2,1],'o',color ='salmon',
label='U[i,j+1]')
graf_3D.set_title('EDP Parabólica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('t')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()
1.3 Gráficas EDP Parabólica, método explícito en 2D con animación
La animación en 2D para la EDP parabólica método explícito se añade, luego de comentar el #plt.show() anterior
# GRAFICA CON ANIMACION 2D--------# import matplotlib.pyplot as pltimport matplotlib.animation as animation
unmetodo = 'EDP Parabólica - Método explícito'
narchivo = 'EdpParabolica'# nombre archivo GIF# Parametros de trama/foto
retardo = 500 # milisegundos entre tramas# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()
# Lineas y puntos base
linea_f, = graf_ani.plot(xi,u[:,0])
graf_ani.axhline(0, color='k') # Eje en x=0
maximoy = np.max(u)
txt_y = maximoy*(1-0.1)
txt_f = graf_ani.text(xi[0],txt_y,'tiempo:',
horizontalalignment='left')
# Configura gráfica
graf_ani.set_xlim([xi[0],xi[ultimox]])
graf_ani.set_ylim([ np.min(u),maximoy])
graf_ani.set_title(unmetodo)
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('u')
#graf_ani.legend()
graf_ani.grid()
# Nueva Tramadefunatrama(i,xi,u):
# actualiza cada linea
linea_f.set_ydata(u[:,i])
linea_f.set_xdata(xi)
linea_f.set_label('tiempo linea: '+str(i))
if (i<=9): # color de la línea
lineacolor = 'C'+str(i)
else:
numcolor = i%10
lineacolor = 'C'+str(numcolor)
linea_f.set_color(lineacolor)
txt_f.set_text('tiempo['+str(i)+'] = ' + str(i*dt))
return linea_f, txt_f,
# Limpia Trama anteriordeflimpiatrama():
linea_f.set_ydata(np.ma.array(xi, mask=True))
linea_f.set_label('')
txt_f.set_text('')
return linea_f, txt_f,
# contador de tramas
i = np.arange(0,len(tk),1)
ani = animation.FuncAnimation(fig_ani, unatrama,
i ,
fargs=(xi,u),
init_func=limpiatrama,
interval=retardo,
blit=True)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'.mp4')
plt.draw()
plt.show()
1.4 Gráficas EDP Parabólica, método explícito en 3D con animación
Para el caso de animación 3D, se cambia la sección de animación del algoritmo anterior por:
# GRAFICA CON ANIMACION 3D--------# import matplotlib.pyplot as pltimport matplotlib.animation as animation
unmetodo = 'EDP Parabólica - Método explícito'
narchivo = 'EdpParabolica2'# nombre archivo GIF# Parametros de trama/foto
retardo = 500 # milisegundos entre tramas# GRAFICA animada en fig_ani
fig_ani3D = plt.figure()
graf_ani3 = fig_ani3D.add_subplot(projection='3d')
# Lineas y puntos base##graf_ani3.plot_wireframe(Xi,Yi,U,## color ='blue',label='EDP Parabólica')
punto_a, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',color ='green')
punto_b, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',
color ='salmon',label='U[i,j+1]')
linea_c, = graf_ani3.plot(xi[0],tk[0],U[1,0],'-',color ='blue')
# Configura gráfica
graf_ani3.set_xlim(min(xi),max(xi))
graf_ani3.set_ylim(min(tk),max(tk))
graf_ani3.set_title(unmetodo)
graf_ani3.set_xlabel('x')
graf_ani3.set_ylabel('t')
graf_ani3.set_zlabel('U')
graf_ani3.legend()
graf_ani3.grid()
graf_ani3.view_init(35, -45)
mallaEDP = None# Nueva Tramadefunatrama(i,xi,tk,U,k):
f = i//k
c = i%k
# actualiza cada punto
punto_a.set_data([[xi[c],xi[c+1],xi[c+2]],
[tk[f],tk[f],tk[f]]])
punto_a.set_3d_properties([U[f,c],U[f,c+1],U[f,c+2]])
punto_b.set_data([xi[c+1]],[tk[f+1]])
punto_b.set_3d_properties([U[f+1,c+1]])
linea_c.set_data(xi[0:c+2],[tk[f+1]]*(c+2))
linea_c.set_3d_properties([U[f+1,0:c+2]])
global mallaEDP
if mallaEDP:
graf_ani3.collections.remove(mallaEDP)
mallaEDP = graf_ani3.plot_wireframe(Xi[0:f+1,:],
Yi[0:f+1,:],
U[0:f+1,:],
color ='blue',
label='EDP Parabólica')
return (punto_a,punto_b,linea_c)
# contador de tramas
k = len(xi)-2
k2 = k*(len(tk)-1)
i = np.arange(0,k2,1)
ani = animation.FuncAnimation(fig_ani3D, unatrama,
i ,
fargs=(xi,tk,U,k),
interval=retardo,
blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'.mp4')#plt.draw()
plt.show()
Puesto que la solución de la tabla de resultados en U es la misma que el los dos métodos explícito e implícito, solo se adjunta la parte de la animación.
Para el caso del la gráfica 3D del método implícito se cambia la sección de gráfica animada por:
# GRAFICA CON ANIMACION 3D--------# import matplotlib.pyplot as pltimport matplotlib.animation as animation
unmetodo = 'EDP Parabólica - Método implícito'
narchivo = 'EdpParabImpli'# nombre archivo GIF# Parametros de trama/foto
retardo = 500 # milisegundos entre tramas# GRAFICA animada en fig_ani
fig_ani3D = plt.figure()
graf_ani3 = fig_ani3D.add_subplot(projection='3d')
# Lineas y puntos base##graf_ani3.plot_wireframe(Xi,Yi,U,## color ='blue',label='EDP Parabólica')
punto_a, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',color ='salmon')
punto_b, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',
color ='green',label='U[i,j-1]')
linea_c, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',color ='blue')
# Configura gráfica
graf_ani3.set_xlim(min(xi),max(xi))
graf_ani3.set_ylim(min(tk),max(tk))
graf_ani3.set_title(unmetodo)
graf_ani3.set_xlabel('x')
graf_ani3.set_ylabel('t')
graf_ani3.set_zlabel('U')
graf_ani3.legend()
graf_ani3.grid()
graf_ani3.view_init(35, -45)
mallaEDP = None# Nueva Tramadefunatrama(i,xi,tk,U,k):
f = i//k
c = i%k
# actualiza cada punto
punto_a.set_data([[xi[c],xi[c+1],xi[c+2]],
[tk[f+1],tk[f+1],tk[f+1]]])
punto_a.set_3d_properties([U[f+1,c],U[f+1,c+1],U[f+1,c+2]])
punto_b.set_data([xi[c+1]],[tk[f]])
punto_b.set_3d_properties([U[f,c+1]])
linea_c.set_data(xi[0:c],[tk[f+1]]*(c))
linea_c.set_3d_properties([U[f+1,0:c]])
global mallaEDP
if mallaEDP:
graf_ani3.collections.remove(mallaEDP)
mallaEDP = graf_ani3.plot_wireframe(Xi[0:f+1,:],
Yi[0:f+1,:],
U[0:f+1,:],
color ='blue',
label='EDP Parabólica')
return (punto_a,punto_b,linea_c)
# contador de tramas
k = len(xi)-2
k2 = k*(len(tk)-1)
i = np.arange(0,k2,1)
ani = animation.FuncAnimation(fig_ani3D, unatrama,
i ,
fargs=(xi,tk,U,k),
interval=retardo,
blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'.mp4')#plt.draw()
plt.show()
Al algoritmo desarrollado en EDP Elípticas método iterativo, se crean una lista para almacenar los valores de la matriz U por cada iteración, se usan para graficar cada iteración.
# Ecuaciones Diferenciales Parciales# Elipticas. Método iterativoimport numpy as np
# INGRESO# Condiciones iniciales en los bordes
Ta = 60 # izquierda
Tb = 25 # derecha#Tc = 50 # inferior
fxc = lambda x: 50+0*x # función inicial inferior# Td = 70 # superior
fxd = lambda x: 70+0*x # función inicial superior# dimensiones de la placa
x0 = 0 # longitud en x
xn = 2
y0 = 0 # longitud en y
yn = 1.5
# discretiza, supone dx=dy
dx = 0.25 # Tamaño de paso
dy = dx
iteramax = 100
tolera = 0.0001
verdigitos = 2 # decimales a mostrar en tabla de resultados# 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)
# llena u con valores en fronteras
u[0,:] = Ta # izquierda
u[n-1,:] = Tb # derecha
fic = fxc(xi)
u[:,0] = fic # Tc inferior
fid = fxd(xi)
u[:,m-1] = fid # Td superior# valor inicial de iteración dentro de u# promedio = (Ta+Tb+Tc+Td)/4
promedio = (Ta+Tb+np.max(fic)+np.max(fid))/4
u[1:n-1,1:m-1] = promedio
np.set_printoptions(precision=verdigitos)
print('u inicial:')
print(np.transpose(u))
# iterar puntos interiores
U_3D = [np.copy(u)]
U_error = ['--']
itera = 0
converge = Falsewhile (itera<=iteramax and converge==False):
itera = itera +1
Uantes = 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])/4
diferencia = u - Uantes
errado_u = np.max(np.abs(diferencia))
if (errado_u<tolera):
converge=True
U_3D.append(np.copy(u))
U_error.append(np.around(errado_u,verdigitos))
if itera<=3:
np.set_printoptions(precision=verdigitos)
print('itera:',itera-1,' ; ','u:')
print(np.transpose(u))
print('errado_u: ', errado_u)
if itera==4:
print('... continua')
# SALIDAprint('Método Iterativo EDP Elíptica')
print('iteraciones:',itera,' converge = ', converge)
print('xi:', xi)
print('yj:', yj)
print('Tabla de resultados en malla EDP Elíptica')
print('j, U[i,j]')
for j inrange(m-1,-1,-1):
print(j,u[:,j])
# GRAFICA en 3D ------import matplotlib.pyplot as plt
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U,
color ='blue')
graf_3D.plot(Xi[1,0],Yi[1,0],U[1,0],'o',color ='orange')
graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='red',
label='U[i,j]')
graf_3D.plot(Xi[1,2],Yi[1,2],U[1,2],'o',color ='salmon')
graf_3D.plot(Xi[0,1],Yi[0,1],U[0,1],'o',color ='green')
graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')
graf_3D.set_title('EDP elíptica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('y')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()
# GRAFICA CON ANIMACION 3D--------# import matplotlib.pyplot as pltimport matplotlib.animation as animation
unmetodo = 'EDP Eliptica - Método iterativo'
narchivo = 'EdpElipticaItera'# nombre archivo GIF# Parametros de trama/foto
retardo = 500 # milisegundos entre tramas# GRAFICA animada en fig_ani
fig_ani3D = plt.figure()
graf_ani3 = fig_ani3D.add_subplot(projection='3d')
# Lineas y puntos base
punto_a, = graf_ani3.plot(xi[0],yj[0],U[0,0],'.',color ='blue')
# Configura gráfica
graf_ani3.set_xlim(min(xi),max(xi))
graf_ani3.set_ylim(min(yj),max(yj))
graf_ani3.set_title(unmetodo)
graf_ani3.set_xlabel('x')
graf_ani3.set_ylabel('t')
graf_ani3.set_zlabel('U')
# graf_ani3.legend()
graf_ani3.grid()
graf_ani3.view_init(35, -45)
etiqueta_i ='itera: '+str(0)
etiqueta_err = 'errado['+str(0)+']: '+str(U_error[0])
txt_i = graf_ani3.text2D(0.05, 0.95, etiqueta_i,
transform=graf_ani3.transAxes)
txt_errado = graf_ani3.text2D(0.05, 0.90, etiqueta_err,
transform=graf_ani3.transAxes)
mallaEDP = None# Nueva Tramadefunatrama(i,xi,yj,U_3D,U_error):
etiqueta_i ='itera: '+str(i)
txt_i.set_text(etiqueta_i)
etiqueta_err = 'errado['+str(i)+']: '+str(U_error[i])
txt_errado.set_text(etiqueta_err)
global mallaEDP
if mallaEDP:
graf_ani3.collections.remove(mallaEDP)
mallaEDP = graf_ani3.plot_wireframe(Xi,Yi,
np.transpose(U_3D[i]),
color ='blue')
return (txt_i,txt_errado,)
# contador de tramas
k = 11 #len(U_3D)
i = np.arange(0,k,1)
ani = animation.FuncAnimation(fig_ani3D, unatrama,
i ,
fargs=(xi,yj,U_3D,U_error),
interval=retardo,
blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='imagemagick')
#ani.save(narchivo+'.mp4')#plt.draw()
plt.show()
Continuando con el ejercicio de la sección anterior de una ecuación diferencial ordinaria, lineal y de coeficientes constantes, se obtuvo la solución homogénea o ecuación complementaria con dsolve(). Para el desarrollo analítico de la solución se requieren describir los pasos tales como la ecuación auxiliar, como se describe a continuación como un ejercicio de análisis de expresiones con Sympy.
Ejemplo 1. Ecuación diferencial de un circuito RLC, ecuación complementaria.
El circuito RLC con entrada x(t) de la figura tiene una corriente y(t) o salida del sistema que se representa por medio de una ecuación diferencial.
dt2d2y(t)+3dtdy(t)+2y(t)=dtdx(t)
Las condiciones iniciales del sistema a tiempo t=0 son y(0)=0 , y'(0)=-5
1. Ecuación diferencial y condiciones de frontera o iniciales
La ecuación diferencial del ejercicio con Sympy se escribe con el operador sym.diff(y,t,1). Indicando la variable independiente ‘t‘, que ‘y‘ es una función y el orden de la derivada con 1. esquema que se mantiene para la descripción de las condiciones iniciales.
# ecuacion: lado izquierdo = lado derecho# Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)
# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
sym.diff(y(t),t,1).subs(t,0) : -5}
2. Clasificar la ecuación diferencial ordinaria
Sympy permite revisar el tipo de EDO a resolver usando la instrucción:
Para el análisis de la respuesta del circuito, se inicia haciendo la entrada x(t)=0, también conocida como ecuación para «respuesta a entrada cero» o ecuación homogenea.
dt2d2y(t)+3dtdy(t)+2y(t)=0
En el algoritmo, La ecuación homogénea se obtiene al substituir x(t)=0 en cada lado de la ecuación, que también es un paso para encontrar la respuesta a entrada cero. La ecuacion homogenea se escribe de la forma f(t)=0.
En la ecuación homogenea , se procede a sustituir en la ecuación el operador dy/dt de la derivada con una variable r elevada al orden de la derivada de cada término . El resultado es una ecuación algebraica que se analiza encontrando las raíces de r.
r2+3r+2=0
Los valores de ‘r’ que resuelven la ecuación permiten estimar la forma de la solución para y(t) conocida como la ecuación general
Se usan diferentes formas para mostrar la ecuación auxiliar, pues en algunos casos se requiere usar la forma de factores, en otros casos los valores de las raíces y las veces que se producen. Formas usadas para generar diagramas de polos y ceros, o expresiones de transformada de Laplace. Motivo por el que los resultados se los presenta en un diccionario, y asi usar la respuesta que sea de interés en cada caso.
# Ecuación Diferencial Ordinaria EDO# ecuación auxiliar o característica # Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dximport sympy as sym
# INGRESO
t = sym.Symbol('t', real=True)
r = sym.Symbol('r')
y = sym.Function('y')
x = sym.Function('x')
# ecuacion: lado izquierdo = lado derecho# Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)
# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
sym.diff(y(t),t,1).subs(t,0) : -5}
# PROCEDIMIENTOdefedo_lineal_auxiliar(ecuacion,
t = sym.Symbol('t'),r = sym.Symbol('r'),
y = sym.Function('y'),x = sym.Function('x')):
''' ecuacion auxiliar o caracteristica de EDO
t independiente
'''# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0
homogenea = sym.expand(homogenea,t)
# ecuación auxiliar o característica
Q = 0*r
term_suma = sym.Add.make_args(homogenea)
for term_k in term_suma:
orden_k = sym.ode_order(term_k,y)
coef = 1 # coefientes del término suma
factor_mul = sym.Mul.make_args(term_k)
for factor_k in factor_mul:
cond = factor_k.has(sym.Derivative)
cond = cond or factor_k.has(y(t))
ifnot(cond):
coef = coef*factor_k
Q = Q + coef*(r**orden_k)
# Q factores y raices
Q_factor = sym.factor(Q,r)
Q_poly = sym.poly(Q,r)
Q_raiz = sym.roots(Q_poly)
auxiliar = {'homogenea' : sym.Eq(homogenea,0),
'auxiliar' : Q,
'Q' : Q,
'Q_factor' : Q_factor,
'Q_raiz' : Q_raiz }
return(auxiliar)
# analiza la ecuación diferencial
edo_tipo = sym.classify_ode(ecuacion, y(t))
auxiliar = edo_lineal_auxiliar(ecuacion,t,r)
# SALIDAprint('clasifica EDO:')
for untipo in edo_tipo:
print(' ',untipo)
print('ecuacion auxiliar:')
for entrada in auxiliar:
print(' ',entrada,':',auxiliar[entrada])
Otra forma de mostrar el resultado es usando un procedimiento creado para mostrar elementos del diccionario según corresponda a cada tipo. el procedimiento se adjunta al final.
Las funciones se incorporan a los algoritmos del curso en matg1052.py
6. Ecuación diferencial, ecuación general y complementaria con Sympy-Python
La ecuación general se encuentra con la instrucción sym.dsolve(), sin condiciones iniciales, por que el resultado presenta constantes Ci por determinar.
# solucion general de ecuación homogénea
general = sym.dsolve(homogenea, y(t))
general = general.expand()
el resultado para el ejercicio es
general :
-t -2*t
y(t) = C1*e + C2*e
Para encontrar los valores de las constantes, se aplica cada una de las condiciones iniciales a la ecuación general, obteniendo un sistema de ecuaciones.
# Aplica condiciones iniciales o de frontera
eq_condicion = []
for cond_k in y_cond: # cada condición
valor_k = y_cond[cond_k]
orden_k = sym.ode_order(cond_k,y)
if orden_k==0: # condicion frontera
t_k = cond_k.args[0] # f(t_k)
expr_k = general.rhs.subs(t,t_k)
else: # orden_k>0# f.diff(t,orden_k).subs(t,t_k)
subs_param = cond_k.args[2] # en valores
t_k = subs_param.args[0] # primer valor
dyk = general.rhs.diff(t,orden_k)
expr_k = dyk.subs(t,t_k)
eq_condicion.append(sym.Eq(valor_k,expr_k))
con el siguiente resultado:
eq_condicion :
0 = C1 + C2
-5 = -C1 - 2*C2
El sistema de ecuaciones se resuelve con la instrucción
constante = sym.solve(eq_condicion)
que entrega un diccionario con cada valor de la constante
constante : {C1: -5, C2: 5}
La ecuación complementaria se obtiene al sustituir en la ecuación general los valores de las constantes.
Las instrucciones detalladas se presentan en el algoritmo integrado.
# Solución complementaria a una Ecuación Diferencial Ordinaria EDO# Lathi 2.1.a pdf 155, (D^2+ 3D + 2)y = Dximport sympy as sym
equivalentes = [{'DiracDelta': lambda x: 1*(x==0)},
{'Heaviside': lambda x,y: np.heaviside(x, 1)},
'numpy',]
import matg1052 as fcnm
# INGRESO
t = sym.Symbol('t', real=True)
r = sym.Symbol('r')
y = sym.Function('y')
x = sym.Function('x')
# ecuacion: lado izquierdo = lado derecho# Left Hand Side = Right Hand Side
LHS = sym.diff(y(t),t,2) + 3*sym.diff(y(t),t,1) + 2*y(t)
RHS = sym.diff(x(t),t,1)
ecuacion = sym.Eq(LHS,RHS)
# condiciones de frontera o iniciales
y_cond = {y(0) : 0,
sym.diff(y(t),t,1).subs(t,0) : -5}
# PROCEDIMIENTOdefedo_lineal_complemento(ecuacion,y_cond):
# ecuación homogénea x(t)=0, entrada cero
RHSx0 = ecuacion.rhs.subs(x(t),0).doit()
LHSx0 = ecuacion.lhs.subs(x(t),0).doit()
homogenea = LHSx0 - RHSx0
# solucion general de ecuación homogénea
general = sym.dsolve(homogenea, y(t))
general = general.expand()
# Aplica condiciones iniciales o de frontera
eq_condicion = []
for cond_k in y_cond: # cada condición
valor_k = y_cond[cond_k]
orden_k = sym.ode_order(cond_k,y)
if orden_k==0: # condicion frontera
t_k = cond_k.args[0] # f(t_k)
expr_k = general.rhs.subs(t,t_k)
else: # orden_k>0
subs_param = cond_k.args[2] # en valores
t_k = subs_param.args[0] # primer valor
dyk = sym.diff(general.rhs,t,orden_k)
expr_k = dyk.subs(t,t_k)
eq_condicion.append(sym.Eq(valor_k,expr_k))
constante = sym.solve(eq_condicion)
# ecuacion complementaria# reemplaza las constantes en general
yc = general
for Ci in constante:
yc = yc.subs(Ci, constante[Ci])
complemento = {'homogenea' : sym.Eq(homogenea,0),
'general' : general,
'eq_condicion' : eq_condicion,
'constante' : constante,
'complementaria' : yc}
return(complemento)
# analiza ecuacion diferencial
edo_tipo = sym.classify_ode(ecuacion, y(t))
auxiliar = fcnm.edo_lineal_auxiliar(ecuacion,t,r)
complemento = edo_lineal_complemento(ecuacion,y_cond)
yc = complemento['complementaria'].rhs
# SALIDAprint('clasifica EDO:')
for elemento in edo_tipo:
print(' ',elemento)
fcnm.print_resultado_dict(auxiliar)
fcnm.print_resultado_dict(complemento)
7. Gráfica de la ecuación complementaria yc
Para la gráfica se procede a convertir la ecuación yc en la forma numérica con sym.lambdify(). Se usa un intervalo de tiempo entre [0,5] y 51 muestras con el siguiente resultado:
8. Mostrar el resultado del diccionario según el tipo de entrada
Para mostrar los resultados de una forma más fácil de leer, se usará sym.pprint() para las ecuaciones, y print() simple para los demás elementos del resultado
defprint_resultado_dict(resultado):
''' print de diccionario resultado
formato de pantalla
'''for entrada in resultado:
tipo = type(resultado[entrada])
if tipo == sym.core.relational.Equality:
print(entrada,':')
sym.pprint(resultado[entrada])
elif tipo==listor tipo==tuple:
tipoelem = type(resultado[entrada][0])
if tipoelem==sym.core.relational.Equality:
print(entrada,':')
for fila in resultado[entrada]:
sym.pprint(fila)
else:
print(entrada,':')
for fila in resultado[entrada]:
print(' ',fila)
else:
print(entrada,':',resultado[entrada])
return()
Para presentar la parte analítica puede seleccionar menos puntos, se busca mostrar la aplicación del algoritmo, no necesariamente cuan largo puede ser.
Tema 3 (35 puntos) Suponga que se requiere planificar la operación de recolección de restos del accidente del sumergible ocurrido en junio del 2023 usando un robot autónomo sumergible.
Existen limitaciones debidas a la profundidad del lecho submarino, las fuertes corrientes que fluyen desde el oeste, limitada energía en las baterías del robot, etc.
El recorrido del robot se puede programar usando un polinomio para el intervalo de las coordenadas presentadas.
Se estima poder realizar tan solo dos recorridos desde el oeste hacia el este.
a. Plantee el ejercicio usando interpolación polinómica para realizar un recorrido para la recolección de restos, describiendo los criterios usados para la selección de los puntos.
b. Desarrolle un primer polinomio usando el método de diferencias divididas de Newton
c. Excluyendo los puntos usados en el literal anterior, desarrolle otro polinomio usando el método de Lagrange.
d. Grafique los resultados obtenidos para cada caso.
e. Determine el error de cada polinomio como la distancia entre cada punto por el que no pasa el polinomio y la coordenada equivalente en el polinomio. Realice observaciones a los resultados.
Nota: Los valores de la tabla se establecen para el ejercicio y no corresponden a una referencia publicada.
Rúbrica: literal a (5 puntos), literal b (10 puntos), literal c (10 puntos), literal d (5 puntos), literal e (5 puntos)
Referencia: [1] La Guardia Costera de EE.UU. anuncia una investigación oficial sobre la implosión del sumergible Titán. RTVE.es / EFE. 26.06.2023 https://www.rtve.es/noticias/20230626/eeuu-anuncia-investigacion-oficial-sumergible-titan/2450440.shtml
[2] Así será la recuperación de los restos del sumergible Titán, según un experto militar. CNN en Español. 23 jun 2023.
Para el ejercicio solo es posible plantear tres ecuaciones, se muestran datos solo para tres semanas. Se tienen 4 incógnitas, por lo que tendría infinitas soluciones y no se podría resolver.
Sin embargo desde el punto de vista práctico se puede usar una variable libre, considerando algunos criterios como:
– que en la nota se da el precio para la docena de huevos a 2.5 USD
– que la variable para huevos es la de menor cantidad se puede incurrir en un menor error si el precio ha variado en los últimos días respecto a lo que se podría haber tenido como referencia.
luego, observando la segunda columna desde el elemento de la diagonal hacia abajo, se observa que el mayor valor en magnitud ya se encuentra en la diagonal. No se requieren intercambios de fila para la tercera columna.
literal d
Para el método iterativo de Gauss-Seidel se requiere el vector inicial, tomando las observaciones de la nota. Se da un precio de 50 USD por un quintal métrico de 100Kg, por lo que se estima que el precio por kilogramo ronda 50/100= 0.5 USD. Se indica además que los demás valores se encuentran en el mismo orden de magnitud, por lo que se tiene como vector inicial
X0 = [0.5, 0.5, 0.5]
Para un método iterativo se despejan las ecuaciones como:
Para el cálculo del error se considera que se busca un precio, lo regular es considerar el centavo, es decir una tolerancia 10-2. Para el desarrollo se considera usar cuatro decimales por si se considera truncar o redondear.
Tema 2 (35 puntos) Un restaurante usa materia prima como: arroz, cebolla, papa, huevos, leche, etc.
Sin embargo los precios han presentado variaciones en el mercado presentando como causas: las lluvias fuertes, el fenómeno del niño, daños en las cosechas y variaciones en insumos por la guerra Ucrania.
Un proveedor mayorista ofrece “combo” de productos por un valor total registrado en varias semanas como:
arroz [Kg]
cebolla [Kg]
papa [Kg]
huevo [docena]
Total Pagado [USD]
Semana 1
500
600
400
90
1660
Semana 2
800
450
300
100
1825
Semana 3
400
300
600
80
1430
Nota: Los valores de la tabla se establecen para el ejercicio y no corresponden a una referencia publicada.
Se requiere conocer el precio unitario de los insumos para estimar el presupuesto operativo semanal del restaurante.
a. Plantee el ejercicio usando los datos de la tabla. Considere justificar el uso de una variable libre.
b. Establezca la forma matricial del sistema de ecuaciones y como matriz aumentada
c. De ser necesario realice el pivoteo parcial por filas
d. Use el método iterativo de Gauss-Seidel, realizando al menos 3 iteraciones. Estime la tolerancia, y justifique.
e. Comente sobre la convergencia del método y justifique sus observaciones usando los errores entre iteraciones.
Nota: considere el precio de 50 USD para el quintal métrico (100Kg) de arroz, y precios menores del mismo orden de magnitud para los demás productos. La docena de huevo a 2,5 USD.
Rúbrica: literal a (2 puntos), literal b (5 puntos), literal c (3puntos), literal d (15 puntos), literal e (5 puntos)
Referencia: [1] Papa, arroz, huevos, cebolla y gas escasean en las islas Galápagos por problemas de abastecimiento. Eluniverso.com. 9 Junio 2023. https://www.eluniverso.com/noticias/ecuador/escasez-de-productos-en-islas-galapagos-persiste-nota/
[2] El Gobierno importará arroz para controlar que el precio no suba por especulación o escasez. Ecuavisa. 16 junio 2023. https://www.ecuavisa.com/noticias/ecuador/el-gobierno-importara-arroz-para-controlar-que-el-precio-no-suba-por-especulacion-o-escasez-CK5396253
[3] Estos son los alimentos de la canasta básica que se han encarecido por las lluvias. Ecuavisa. 05 junio 2023. https://www.ecuavisa.com/la-noticia-a-fondo/estos-son-los-alimentos-de-la-canasta-basica-que-se-han-encarecido-por-las-lluvias-FX5303236
[4] Cebolla, arroz, piña y tomate son los productos que han subido su precio en mercados de Quito. La Hora. Junio 23, 2023. https://www.lahora.com.ec/pais/cebolla-arroz-pina-y-tomate-son-los-productos-que-han-subido-su-precio-en-mercados-de-quito/
[5] Escasez de arroz ya se refleja en los mercados de Ecuador | Televistazo | Ecuavisa. 14 jun 2023