1. Ejercicio
Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694
\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0Valores de frontera: Ta = 60, Tb = 25, Tc = 50, Td = 70
Longitud en x0 = 0, xn = 2, y0 = 0, yn = 1.5
Tamaño de paso dx = 0.25, dy = dx
iteramax=100, tolera = 0.0001
(ecuación de Laplace, Ecuación de Poisson con f(x,y)=0)
2. EDP Elípticas contínua a discreta
El resultado de las operaciones en la expresión usando el algoritmo es:
EDP Elíptica contínua a discreta
ordenDx : 2
ordenDy : 2
edp=0 :
2 2
∂ ∂
───(u(x, y)) + ───(u(x, y)) = 0
2 2
∂x ∂y
K_ : 1
discreta=0 :
-2⋅U(i, j) + U(i, j - 1) + U(i, j + 1) -2⋅U(i, j) + U(i - 1, j) + U(i + 1, j)
────────────────────────────────────── + ──────────────────────────────────────
2 2
Dy Dx
Lambda_L :
2
Dy
───
2
Dx
Lambda L_k : 1
(discreta=0)*Dy**ordeny :
2 2 2
2⋅Dy ⋅U(i, j) Dy ⋅U(i - 1, j) Dy ⋅ U(i + 1, j)
-2⋅U(i, j) + U(i, j - 1) + U(i, j + 1) - ───────────── + ─────────────── + ───────────────
2 2 2
Dx Dx Dx
discreta_L = 0 :
L⋅U(i - 1, j) + L⋅U(i + 1, j) + (-2⋅L - 2)⋅U(i, j) + U(i, j - 1) + U(i, j + 1) = 0
discreta :
-4⋅U(i, j) + U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j) = 0
discreta_iterativa
U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j)
U(i, j) = ─────────────────────────────────────────────────────
4
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.
# INGRESO
fxy = 0*x+0*y # f(x,y) = 0 , ecuacion de Poisson
# ecuacion: LHS=RHS
LHS = sym.diff(u,x,2) + sym.diff(u,y,2)
RHS = fxy
EDP = sym.Eq(LHS-RHS,0)
las diferencias divididas se ingresan como un diccionario:
# centrada, centrada, atras
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.
# 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
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.
3. Algoritmo contínua a discreta en Python
Los resultados al aplicar una operación a la expresión se guardan en un diccionario, para revisar cada paso intermedio al final
# Ecuaciones Diferenciales Parciales Elipticas
# EDP Elípticas contínua a discreta con Sympy
import numpy as np
import sympy as sym
# u(x,y) funciones continuas y variables simbólicas usadas
x = sym.Symbol('x',real=True)
y = sym.Symbol('y',real=True)
u = sym.Function('u')(x,y) # funcion
f = sym.Function('f')(x,y) # funcion complemento
# U[i,j] funciones discretas y variables simbólicas usadas
i = sym.Symbol('i',integer=True,positive=True) # indices
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)
# centrada, centrada, atras
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
fxd = lambda x: 70 +0*x # superior
# dimensiones de la placa
x0 = 0 # longitud en x
xn = 2
y0 = 0 # longitud en y
yn = 1.5
# muestreo en ejes, discreto, supone dx=dy
dx = 0.25 # Tamaño de paso
dy = dx # supone dx=dy
iteramax = 100 # revisa convergencia
tolera = 0.0001
verdigitos = 2 # para mostrar en tabla de resultados
casicero = 1e-15 # para redondeo de términos en ecuacion
# PROCEDIMIENTO
def edp_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
if not(edp.rhs==0):
LHS = LHS-RHS
RHS = 0
edp = sym.Eq(LHS,RHS)
# 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 # guarda en resultados
resultado['ordenDy'] = ordenDy
# coeficiente derivada orden mayor a 1 (d2u/dx2)
coeff_x = edp_coef_Dx(edp,x,ordenDx)
if not(coeff_x==1):
LHS = LHS/coeff_x
RHS = RHS/coeff_x
edp = sym.Eq(LHS,RHS)
K_ = edp_coef_Dx(edp,y,ordenDy)
if abs(K_)%1<casicero: # si es entero
K_ = int(K_)
resultado['edp=0'] = edp
resultado['K_'] = K_
discreta = edp.lhs # EDP discreta
for derivada in dif_dividida: # reemplaza diferencia dividida
discreta = discreta.subs(derivada,dif_dividida[derivada])
resultado['discreta=0'] = discreta
return (resultado)
def edp_coef_Dx(edp,x,ordenx):
''' Extrae el coeficiente de la derivada Dx de ordenx
edp es la ecuación como lhs=rhs
'''
coeff_x = 1.0 # valor inicial
# 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=False
# separa cada factor de término
for factor_k in factor_Mul:
if not(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==True:
coeff_x = coeff_x*coeff_temp
return(coeff_x)
def redondea_coef(ecuacion, precision=6,casicero = 1e-15):
''' redondea coeficientes de términos suma de una ecuacion
ecuación como lhs=rhs
'''
tipo = type(ecuacion)
tipo_eq = False
if tipo == sym.core.relational.Equality:
RHS = ecuacion.rhs
ecuacion = ecuacion.lhs
tipo = type(ecuacion)
tipo_eq = True
if tipo == sym.core.add.Add: # términos suma de ecuacion
term_sum = sym.Add.make_args(ecuacion)
ecuacion = sym.S.Zero # vacia
for term_k in term_sum:
# factor multiplicativo de termino suma
term_mul = sym.Mul.make_args(term_k)
producto = sym.S.One
for factor in term_mul:
if not(factor.has(sym.Symbol)): # es numerico
factor = np.around(float(factor),precision)
if (abs(factor)%1)<casicero: # si es entero
factor = int(factor)
producto = producto*factor
ecuacion = ecuacion + producto
if tipo == sym.core.mul.Mul: # termino único, busca factores
term_mul = sym.Mul.make_args(ecuacion)
producto = sym.S.One
for factor in term_mul:
if not(factor.has(sym.Symbol)): # es numerico
factor = np.around(float(factor),precision)
if (abs(factor)%1)<casicero: # si es entero
factor = int(factor)
producto = producto*factor
ecuacion = producto
if tipo == float: # # solo un numero
if (abs(ecuacion)%1)<casicero:
ecuacion = int(ecuacion)
if tipo_eq==True: # era igualdad, integra lhs=rhs
ecuacion = sym.Eq(ecuacion,RHS)
return(ecuacion)
def edp_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)])
if abs(L_k)%1<casicero: # si es entero
L_k = int(L_k)
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)])
discreta_L = redondea_coef(discreta_L)
resultado['discreta'] = discreta_L
return(resultado)
def edp_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 = False # por orden de derivada
if 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=True
if cambiar==True:
term_k = term_k*L/resultado['Lambda_L']
discreta_L = discreta_L + term_k
# simplifica unos con decimal a entero
discreta_L = redondea_coef(discreta_L)
return(discreta_L)
# 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])
4. 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:
discreta_iterativa
U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j)
U(i, j) = ─────────────────────────────────────────────────────
4
iterar i,j:
j: 1 ; i: 1
U(0, 1) U(1, 0) U(1, 2) U(2, 1)
U(1, 1) = ─────── + ─────── + ─────── + ───────
4 4 4 4
U(1, 1) = 53.125
j: 1 ; i: 2
U(1, 1) U(2, 0) U(2, 2) U(3, 1)
U(2, 1) = ─────── + ─────── + ─────── + ───────
4 4 4 4
U(2, 1) = 51.40625
j: 1 ; i: 3
U(2, 1) U(3, 0) U(3, 2) U(4, 1)
U(3, 1) = ─────── + ─────── + ─────── + ───────
4 4 4 4
U(3, 1) = 50.9765625
j: 1 ; i: 4
U(3, 1) U(4, 0) U(4, 2) U(5, 1)
U(4, 1) = ─────── + ─────── + ─────── + ───────
4 4 4 4
U(4, 1) = 50.869140625
j: 1 ; i: 5
U(4, 1) U(5, 0) U(5, 2) U(6, 1)
U(5, 1) = ─────── + ─────── + ─────── + ───────
4 4 4 4
U(5, 1) = 50.84228515625
j: 1 ; i: 6
U(5, 1) U(6, 0) U(6, 2) U(7, 1)
U(6, 1) = ─────── + ─────── + ─────── + ───────
4 4 4 4
U(6, 1) = 50.8355712890625
j: 1 ; i: 7
U(6, 1) U(7, 0) U(7, 2) U(8, 1)
U(7, 1) = ─────── + ─────── + ─────── + ───────
4 4 4 4
U(7, 1) = 44.271392822265625
Método iterativo EDP Elíptica
continuar el desarrollo con:
xi: [0. 0.25 0.5 0.75 1. 1.25 1.5 1.75 2. ]
yj: [0. 0.25 0.5 0.75 1. 1.25 1.5 ]
j, U[i,j]
2 [60. 51.25 51.25 51.25 51.25 51.25 51.25 51.25 25. ]
1 [60. 53.12 51.41 50.98 50.87 50.84 50.84 44.27 25. ]
0 [50. 50. 50. 50. 50. 50. 50. 50. 50.]
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/2,dx)
yj = np.arange(y0,yn+dy/2,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ón
def edp_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 = False ; 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
if not(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ón
if cambiar==True:
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)
# Método iterativo para EDP Elíptica
# separa término U[j,j] del centro,con valor no conocido
buscar = U.subs(j,j)
discreta = sym.solve(discreta,buscar)
discreta = sym.factor_terms(discreta[0])
discreta = sym.Eq(buscar,discreta)
resultado['discreta_iterativa'] = discreta
print('\ndiscreta_iterativa')
sym.pprint(discreta)
# 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 matriz
print('\n iterar i,j:')
j_k = 1
for i_k in range(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.subs({i:i_k,j:j_k}),
'=',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 in range(2,-1,-1):
print(j_k, (u_xy[:,j_k]))
5. Grafica en 3D para EDP Elípticas
A partir de la solución en matriz U_xy, Xi, Yi
# GRAFICA en 3D ------
import matplotlib.pyplot as plt
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi, yj)
U_xy = np.transpose(u_xy) # 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_xy,
color ='blue')
graf_3D.plot(Xi[1,0],Yi[1,0],U_xy[1,0],'o',color ='orange')
graf_3D.plot(Xi[1,1],Yi[1,1],U_xy[1,1],'o',color ='red',
label='U[i,j]')
graf_3D.plot(Xi[1,2],Yi[1,2],U_xy[1,2],'o',color ='salmon')
graf_3D.plot(Xi[0,1],Yi[0,1],U_xy[0,1],'o',color ='green')
graf_3D.plot(Xi[2,1],Yi[2,1],U_xy[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()