EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]
..
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
EDP Parabólica [ contínua a discreta ] || [ sympy_explícito ] [sympy_ implícito ]
..
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í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().
Instrucciones 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])
EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]
3. EDP Parabólica – 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]))
EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]