EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]
..
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)
EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]
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.
Instrucciones 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])
EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]
..
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:
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]))
EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]
..
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()
EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]