1. Ejercicio
Referencia: Chapra 30.2 p888 pdf912, Burden 9Ed 12.2 p725, Rodríguez 10.2 p406
\frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}Valores de frontera: Ta = 60, Tb = 40, T0 = 25
Longitud en x a = 0, b = 1, Constante K= 4
Tamaño de paso dx = 0.1, dt = dx/10
2. EDP Parabólica contínua a discreta
El resultado del algoritmo para el ejercicio propuesto es:
EDP Parabólica contínua a discreta
ordenDx : 2
ordenDy : 1
edp=0 :
2
∂ ∂
───(u(x, y)) - 4⋅──(u(x, y)) = 0
2 ∂y
∂x
K_ : 4
Lambda_L :
Dy
─────
2
4⋅Dx
Lambda L_k : 0.250000000000000
(discreta=0)*Dy**ordeny :
2⋅Dy⋅U(i, j) Dy⋅U(i - 1, j) Dy⋅U(i + 1, j)
4⋅U(i, j) - 4⋅U(i, j + 1) - ──────────── + ────────────── + ──────────────
2 2 2
Dx Dx Dx
discreta_L=0 :
⎛ 2⋅Dy⎞ Dy⋅U(i - 1, j) Dy⋅U(i + 1, j)
⎜4 - ────⎟⋅U(i, j) - 4⋅U(i, j + 1) + ────────────── + ────────────── = 0
⎜ 2 ⎟ 2 2
⎝ Dx ⎠ Dx Dx
discreta :
2⋅U(i, j) - 4⋅U(i, j + 1) + U(i - 1, j) + U(i + 1, j) = 0
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.
# 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)
La ecuación f(x,y) se añade para ser usada con otros ejercicios que puede revisar en la sección de Evaluaciones.
Las expresiones para diferencias divididas seleccionadas se ingresan como un diccionario:
# centrada, adelante
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}
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.
# 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
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) # U[i,j+1] método explícitoEl 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().
3. Algoritmo contínua a discreta en Python
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 Sympy
import numpy as np
import sympy as sym
# u(x,y) funciones continuas y variables simbólicas usadas
t = sym.Symbol('t',real=True) # independiente
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
K = sym.Symbol('K',real=True)
# 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)
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)
# centrada, adelante
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) # U[i,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
# [a,b] dimensiones de la barra
a = 0 # longitud en x
b = 1
y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones
# muestreo en ejes, discreto
x_tramos = 10
dx = (b-a)/x_tramos
dy = dx/10
n_iteraciones = 4 # iteraciones en tiempo
verdigitos = 2 # para mostrar en tabla de resultados
casicero = 1e-15 # para redondeo de términos en ecuacion
# PROCEDIMIENTO --------------------------------------
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 = True
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 # separa lado derecho
ecuacion = ecuacion.lhs # analiza lado izquierdo
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)
resultado = {} # resultados en diccionario
# 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)
LHS = edp.lhs # lado izquierdo de edp
RHS = edp.rhs # lado derecho de edp
if not(coeff_x==1):
LHS = LHS/coeff_x
RHS = RHS/coeff_x
# toda la expresión a la izquierda
edp = sym.Eq(LHS-RHS,0)
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_
# lambda en ecuacion edp
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:
L_k = int(L_k)
resultado['Lambda L_k'] = L_k
discreta = edp.lhs # EDP discreta
for derivada in dif_dividida: # reemplaza diferencia 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) #agrupa 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
# 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])
4. Algoritmo método explícito con Sympy
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:
Método explícito EDP Parabólica
discreta_itera:
U(i, j) U(i - 1, j) U(i + 1, j)
U(i, j + 1) = ─────── + ─────────── + ───────────
2 4 4
iterar i,j:
j: 0 ; i: 1 ; ecuacion: 0
U(0, 0) U(1, 0) U(2, 0)
U(1, 1) = ─────── + ─────── + ───────
4 2 4
U(1, 1) = 33.75
j: 0 ; i: 2 ; ecuacion: 1
U(1, 0) U(2, 0) U(3, 0)
U(2, 1) = ─────── + ─────── + ───────
4 2 4
U(2, 1) = 25.0
j: 0 ; i: 3 ; ecuacion: 2
U(2, 0) U(3, 0) U(4, 0)
U(3, 1) = ─────── + ─────── + ───────
4 2 4
U(3, 1) = 25.0
j: 0 ; i: 4 ; ecuacion: 3
U(3, 0) U(4, 0) U(5, 0)
U(4, 1) = ─────── + ─────── + ───────
4 2 4
U(4, 1) = 25.0
j: 0 ; i: 5 ; ecuacion: 4
U(4, 0) U(5, 0) U(6, 0)
U(5, 1) = ─────── + ─────── + ───────
4 2 4
U(5, 1) = 25.0
j: 0 ; i: 6 ; ecuacion: 5
U(5, 0) U(6, 0) U(7, 0)
U(6, 1) = ─────── + ─────── + ───────
4 2 4
U(6, 1) = 25.0
j: 0 ; i: 7 ; ecuacion: 6
U(6, 0) U(7, 0) U(8, 0)
U(7, 1) = ─────── + ─────── + ───────
4 2 4
U(7, 1) = 25.0
j: 0 ; i: 8 ; ecuacion: 7
U(7, 0) U(8, 0) U(9, 0)
U(8, 1) = ─────── + ─────── + ───────
4 2 4
U(8, 1) = 25.0
j: 0 ; i: 9 ; ecuacion: 8
U(8, 0) U(9, 0) U(10, 0)
U(9, 1) = ─────── + ─────── + ────────
4 2 4
U(9, 1) = 28.75
... continua calculando
x: [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0. 0.01 0.02 0.03 0.04]
j, U[i,j]
4 [60. 42.77 31.29 26.37 25.14 25. 25.06 25.59 27.7 32.62 40. ]
3 [60. 40.86 29.38 25.55 25. 25. 25. 25.23 26.88 31.8 40. ]
2 [60. 38.12 27.19 25. 25. 25. 25. 25. 25.94 30.62 40. ]
1 [60. 33.75 25. 25. 25. 25. 25. 25. 25. 28.75 40. ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
Instrucciones en Python
# Método explí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/2,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
# busca el término no conocido
u_explicito = sym.solve(discreta_L,buscar)
u_explicito = sym.Eq(buscar,u_explicito[0])
resultado['discreta_itera'] = u_explicito
print('\nMétodo explícito EDP Parabólica')
print('discreta_itera:')
sym.pprint(u_explicito)
# iterando
print('\n iterar i,j:')
j_k = 0
for j_k in range(0,m-1,1):
for i_k in range(1,n-1,1):
eq_conteo = j_k*(n-2)+(i_k-1)
discreta_ij = u_explicito.subs({i:i_k,j:j_k})
valorij = 0 # calcula valor de nodo ij
# separa cada término de suma
term_suma = sym.Add.make_args(discreta_ij.rhs)
for term_k in term_suma:
term_ki = 1
if term_k.is_Function:
[i_f,j_f] = term_k.args
term_ki = u_xy[i_f,j_f]
elif term_k.is_Mul: # mas de un factor
factor_Mul = sym.Mul.make_args(term_k)
# separa cada factor de término
for factor_k in factor_Mul:
if factor_k.is_Function:
[i_f,j_f] = factor_k.args
term_ki = term_ki*u_xy[i_f,j_f]
else:
term_ki = term_ki*factor_k
else:
term_ki = term_k
valorij = valorij + term_ki
[i_f,j_f] = discreta_ij.lhs.args
u_xy[i_f,j_f] = float(valorij) # actualiza matriz u
# muestra las primeras m iteraciones
if eq_conteo<(n-2):
print('j: ',j_k,' ; ','i: ',i_k,
' ; ','ecuacion:',eq_conteo)
sym.pprint(discreta_ij)
print(discreta_ij.lhs,'=',float(valorij))
print('... continua calculando')
# SALIDA
np.set_printoptions(precision=verdigitos)
print('\nx:',xi)
print('t:',yj)
print(' j, U[i,j]')
for j_k in range(m-1,-1,-1):
print(j_k, (u_xy[:,j_k]))