EDP Parabólica [ contínua a discreta ] || [ sympy_explícito ] [sympy_ 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 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 de diferencias divididas seleccionadas se ingresan como un diccionario:
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) # 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().
El resultado del algoritmo para el ejercicio propuesto es:
EDP Parabólica contínua a discreta
ordenDx : 2
ordenDy : 1
edp=0 :
2
d d
---(u(x, y)) - 4*--(u(x, y)) = 0
2 dy
dx
K_ : 4.0
Lambda L :
0.25*Dy
-------
2
Dx
Lambda L_k : 0.24999999999999997
(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.0*U(i, j) - 4*U(i, j + 1) + 1.0*U(i - 1, j) + 1.0*U(i + 1, j) = 0
>>>
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 # funciones continuas y variables simbólicas usadas t = sym.Symbol('t',real=True) y = sym.Symbol('y',real=True) x = sym.Symbol('x',real=True) u = sym.Function('u')(x,y) K = sym.Symbol('K',real=True) f = sym.Function('f')(x,y) # funcion complemento # funciones discretas y variables simbólicas usadas i = sym.Symbol('i',integer=True,positive=True) 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) 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) # 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 # dimensiones de la barra a = 0 # longitud en x b = 1 y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones # discretiza, supone dx=dy x_tramos = 10 dx = (b-a)/x_tramos dy = dx/10 verdigitos = 2 # para mostrar en tabla de resultados n_iteraciones = 4 # iteraciones en tiempo # PROCEDIMIENTO def edp_coef_Dx(edp,x,ordenx): ''' Extrae el coeficiente de la derivada Dx de ordenx ''' coeff_x = 1.0 # 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=0 # 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==1: coeff_x = coeff_x*coeff_temp return(coeff_x) resultado={} # expresión todo a la izquierda LHS (Left Hand side) if not(edp.rhs==0): LHS = EDP.lhs RHS = EDP.rhs edp = sym.Eq(LHS,RHS) # para 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 resultado['ordenDy'] = ordenDy # convierte coeficiente d2u/dx2 a 1 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,0) K_ = -edp_coef_Dx(edp,y,ordenDy) resultado['edp=0'] = edp resultado['K_'] = float(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),(1.0,1)]) resultado['Lambda L_k'] = float(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 = 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),(1.0+0*x,1)]) discreta_L = discreta_L.subs(1.0+0*x,1) 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 ] [sympy_ 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 + 1) = 0.5*U(i, j) + 0.25*U(i - 1, j) + 0.25*U(i + 1, j) iterar i,j: j: 0 ; i: 1 ; ecuacion: 0 U(1, 1) = 0.25*U(0, 0) + 0.5*U(1, 0) + 0.25*U(2, 0) U(1, 1) = 33.75 j: 0 ; i: 2 ; ecuacion: 1 U(2, 1) = 0.25*U(1, 0) + 0.5*U(2, 0) + 0.25*U(3, 0) U(2, 1) = 25.0 j: 0 ; i: 3 ; ecuacion: 2 U(3, 1) = 0.25*U(2, 0) + 0.5*U(3, 0) + 0.25*U(4, 0) U(3, 1) = 25.0 j: 0 ; i: 4 ; ecuacion: 3 U(4, 1) = 0.25*U(3, 0) + 0.5*U(4, 0) + 0.25*U(5, 0) U(4, 1) = 25.0 j: 0 ; i: 5 ; ecuacion: 4 U(5, 1) = 0.25*U(4, 0) + 0.5*U(5, 0) + 0.25*U(6, 0) U(5, 1) = 25.0 j: 0 ; i: 6 ; ecuacion: 5 U(6, 1) = 0.25*U(5, 0) + 0.5*U(6, 0) + 0.25*U(7, 0) U(6, 1) = 25.0 j: 0 ; i: 7 ; ecuacion: 6 U(7, 1) = 0.25*U(6, 0) + 0.5*U(7, 0) + 0.25*U(8, 0) U(7, 1) = 25.0 j: 0 ; i: 8 ; ecuacion: 7 U(8, 1) = 0.25*U(7, 0) + 0.5*U(8, 0) + 0.25*U(9, 0) U(8, 1) = 25.0 j: 0 ; i: 9 ; ecuacion: 8 U(9, 1) = 0.25*U(8, 0) + 0.5*U(9, 0) + 0.25*U(10, 0) 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,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 ] [sympy_ 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 atras 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
d d
---(u(x, y)) - 4*--(u(x, y)) = 0
2 dy
dx
K_ : 4.0
Lambda_L :
0.25*Dy
-------
2
Dx
Lambda L_k : 0.24999999999999997
(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.0*L*U(i - 1, j) + 4.0*L*U(i + 1, j) + (-8.0*L - 4)*U(i, j) + 4*U(i, j - 1) =
0
discreta :
-6.0*U(i, j) + 4*U(i, j - 1) + 1.0*U(i - 1, j) + 1.0*U(i + 1, j) = 0
Método implícito EDP Parabólica
j: 1 ; i: 1
1.0*U(0, 1) + 4*U(1, 0) - 6.0*U(1, 1) + 1.0*U(2, 1) = 0
-6.0*U(1, 1) + 1.0*U(2, 1) = -160.0
j: 1 ; i: 2
1.0*U(1, 1) + 4*U(2, 0) - 6.0*U(2, 1) + 1.0*U(3, 1) = 0
1.0*U(1, 1) - 6.0*U(2, 1) + 1.0*U(3, 1) = -100.0
j: 1 ; i: 3
1.0*U(2, 1) + 4*U(3, 0) - 6.0*U(3, 1) + 1.0*U(4, 1) = 0
1.0*U(2, 1) - 6.0*U(3, 1) + 1.0*U(4, 1) = -100.0
j: 1 ; i: 4
1.0*U(3, 1) + 4*U(4, 0) - 6.0*U(4, 1) + 1.0*U(5, 1) = 0
1.0*U(3, 1) - 6.0*U(4, 1) + 1.0*U(5, 1) = -100.0
j: 1 ; i: 5
1.0*U(4, 1) + 4*U(5, 0) - 6.0*U(5, 1) + 1.0*U(6, 1) = 0
1.0*U(4, 1) - 6.0*U(5, 1) + 1.0*U(6, 1) = -100.0
j: 1 ; i: 6
1.0*U(5, 1) + 4*U(6, 0) - 6.0*U(6, 1) + 1.0*U(7, 1) = 0
1.0*U(5, 1) - 6.0*U(6, 1) + 1.0*U(7, 1) = -100.0
j: 1 ; i: 7
1.0*U(6, 1) + 4*U(7, 0) - 6.0*U(7, 1) + 1.0*U(8, 1) = 0
1.0*U(6, 1) - 6.0*U(7, 1) + 1.0*U(8, 1) = -100.0
j: 1 ; i: 8
1.0*U(7, 1) + 4*U(8, 0) - 6.0*U(8, 1) + 1.0*U(9, 1) = 0
1.0*U(7, 1) - 6.0*U(8, 1) + 1.0*U(9, 1) = -100.0
j: 1 ; i: 9
1.0*U(8, 1) + 4*U(9, 0) - 6.0*U(9, 1) + 1.0*U(10, 1) = 0
1.0*U(8, 1) - 6.0*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]
j: 2 ; i: 1
1.0*U(0, 2) + 4*U(1, 1) - 6.0*U(1, 2) + 1.0*U(2, 2) = 0
-6.0*U(1, 2) + 1.0*U(2, 2) = -184.020210038146
...
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 # funciones continuas y variables simbólicas usadas t = sym.Symbol('t',real=True) y = sym.Symbol('y',real=True) x = sym.Symbol('x',real=True) u = sym.Function('u')(x,y) K = sym.Symbol('K',real=True) f = sym.Function('f')(x,y) # funcion complemento # funciones discretas y variables simbólicas usadas i = sym.Symbol('i',integer=True,positive=True) 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) 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} # 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 # dimensiones de la barra a = 0 # longitud en x b = 1 y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones # discretiza x_tramos = 10 dx = (b-a)/x_tramos dy = dx/10 verdigitos = 2 # para mostrar en tabla de resultados n_iteraciones = 4 # iteraciones en tiempo # PROCEDIMIENTO def edp_coef_Dx(edp,x,ordenx): ''' Extrae el coeficiente de la derivada Dx de ordenx ''' coeff_x = 1.0 # 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=0 # 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==1: coeff_x = coeff_x*coeff_temp return(coeff_x) 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 = 0 # 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=1 if cambiar: term_k = term_k*L/resultado['Lambda_L'] discreta_L = discreta_L + term_k # simplifica unos con decimal a entero discreta_L = discreta_L.subs(1.0,1) return(discreta_L) resultado={} # expresión todo a la izquierda LHS (Left Hand side) if not(edp.rhs==0): LHS = EDP.lhs RHS = EDP.rhs edp = sym.Eq(LHS,RHS) # para 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 resultado['ordenDy'] = ordenDy # convierte coeficiente d2u/dx2 a 1 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,0) K_ = -edp_coef_Dx(edp,y,ordenDy) resultado['edp=0'] = edp resultado['K_'] = float(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),(1.0,1)]) resultado['Lambda L_k'] = float(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) 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),(1.0+0*x,1)]) 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,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 # 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 = 0 ; 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: 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 = {} 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 ] [sympy_ implícito ]