EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]
..
4. EDP Parabólica – Método implícito con Sympy
Para el ejercicio presentado en el numeral 1, se desarrolla el método implícito.
Las diferencias finitas centradas y hacia atrás 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:
elif (tipo==dict):
if entrada<2
Los resultados para el algoritmo son:
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 :
4⋅L⋅U(i - 1, j) + 4⋅L⋅U(i + 1, j) + (-8⋅L - 4)⋅U(i, j) + 4⋅U(i, j - 1) = 0
discreta :
-6⋅U(i, j) + 4⋅U(i, j - 1) + U(i - 1, j) + U(i + 1, j) = 0
Método implícito EDP Parabólica
j: 1 ; i: 1
U(0, 1) + 4⋅U(1, 0) - 6⋅U(1, 1) + U(2, 1) = 0
-6⋅U(1, 1) + U(2, 1) = -160.0
j: 1 ; i: 2
U(1, 1) + 4⋅U(2, 0) - 6⋅U(2, 1) + U(3, 1) = 0
U(1, 1) - 6⋅U(2, 1) + U(3, 1) = -100.0
j: 1 ; i: 3
U(2, 1) + 4⋅U(3, 0) - 6⋅U(3, 1) + U(4, 1) = 0
U(2, 1) - 6⋅U(3, 1) + U(4, 1) = -100.0
j: 1 ; i: 4
U(3, 1) + 4⋅U(4, 0) - 6⋅U(4, 1) + U(5, 1) = 0
U(3, 1) - 6⋅U(4, 1) + U(5, 1) = -100.0
j: 1 ; i: 5
U(4, 1) + 4⋅U(5, 0) - 6⋅U(5, 1) + U(6, 1) = 0
U(4, 1) - 6⋅U(5, 1) + U(6, 1) = -100.0
j: 1 ; i: 6
U(5, 1) + 4⋅U(6, 0) - 6⋅U(6, 1) + U(7, 1) = 0
U(5, 1) - 6⋅U(6, 1) + U(7, 1) = -100.0
j: 1 ; i: 7
U(6, 1) + 4⋅U(7, 0) - 6⋅U(7, 1) + U(8, 1) = 0
U(6, 1) - 6⋅U(7, 1) + U(8, 1) = -100.0
j: 1 ; i: 8
U(7, 1) + 4⋅U(8, 0) - 6⋅U(8, 1) + U(9, 1) = 0
U(7, 1) - 6⋅U(8, 1) + U(9, 1) = -100.0
j: 1 ; i: 9
U(8, 1) + 4⋅U(9, 0) - 6⋅U(9, 1) + U(10, 1) = 0
U(8, 1) - 6⋅U(9, 1) = -140.0
A :
[[-6. 1. 0. 0. 0. 0. 0. 0. 0.]
[ 1. -6. 1. 0. 0. 0. 0. 0. 0.]
[ 0. 1. -6. 1. 0. 0. 0. 0. 0.]
[ 0. 0. 1. -6. 1. 0. 0. 0. 0.]
[ 0. 0. 0. 1. -6. 1. 0. 0. 0.]
[ 0. 0. 0. 0. 1. -6. 1. 0. 0.]
[ 0. 0. 0. 0. 0. 1. -6. 1. 0.]
[ 0. 0. 0. 0. 0. 0. 1. -6. 1.]
[ 0. 0. 0. 0. 0. 0. 0. 1. -6.]]
B :
[-160. -100. -100. -100. -100. -100. -100. -100. -140.]
X :
[31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57]
Resultados para U(x,y)
xi: [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
yj: [0. 0.01 0.02 0.03 0.04]
j, U[i,j]
4 [60. 40.67 30.61 26.75 25.51 25.18 25.24 25.76 27.41 31.71 40. ]
3 [60. 38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72 40. ]
2 [60. 35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39 40. ]
1 [60. 31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57 40. ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
Instrucciones en Python
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 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, 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,1): (U - U.subs(j,j-1))/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) 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) 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: # si es entero L_k = int(L_k) resultado['Lambda L_k'] = L_k discreta = edp.lhs # EDP discreta for 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) # 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]) # 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/2,dx) yj = np.arange(y0,n_iteraciones*dy+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] = 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'] 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) # ITERAR para plantear las ecuaciones en [i,j] resultado = {} # resultados en diccionario eq_itera = [] ; tamano = (n-2) A = np.zeros(shape=(tamano,tamano),dtype=float) B = np.zeros(tamano,dtype=float) for j_k in range(1,m,1): # no usar valores en bordes resultado[j_k] ={} for i_k in range(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 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 ]