EDP Parab\u00f3lica<\/p>\n\n\n\n
cont\u00ednua a discreta<\/a><\/p>\n\n\n\n algoritmo discreta<\/a><\/p>\n\n\n\n algoritmo explicito<\/a><\/p>\n<\/div>\n\n\n\n Referencia<\/strong><\/em>: Chapra 30.2 p888 pdf912, Burden 9Ed 12.2 p725, Rodr\u00edguez 10.2 p406<\/p>\n\n\n \\frac{\\partial ^2 u}{\\partial x ^2} = K\\frac{\\partial u}{\\partial t}<\/span>\n\n\n\n Valores de frontera: T<\/strong>a = 60, T<\/strong>b = 40, T<\/strong>0 = 25 EDP Parab\u00f3lica<\/p>\n\n\n\n cont\u00ednua a discreta<\/a><\/p>\n\n\n\n algoritmo discreta<\/a><\/p>\n\n\n\n algoritmo explicito<\/a><\/p>\n<\/div>\n\n\n\n El resultado del algoritmo para el ejercicio propuesto es:<\/p>\n\n\n\n El desarrollo anal\u00edtico comienza en escribir la ecuaci\u00f3n diferencial parcial parab\u00f3lica de la forma cont\u00ednua a la representaci\u00f3n discreta, siguiendo los criterios indicados en la parte te\u00f3rica usando derivadas expresadas en diferencias divididas.<\/p>\n\n\n\n La EDP se escribe en formato Sympy en el bloque de ingreso en dos partes: lado izquierdo (LHS) de la ecuaci\u00f3n y lado derecho (RHS), definiendo u como una funci\u00f3n de variables x,y.<\/p>\n\n\n La ecuaci\u00f3n f(x,y) se a\u00f1ade para ser usada con otros ejercicios que puede revisar en la secci\u00f3n de Evaluaciones.<\/p>\n\n\n\n Las expresiones para dif<\/strong>erencias divididas<\/strong> seleccionadas se ingresan como un diccionario:<\/p>\n\n\n Las condiciones iniciales en los extremos a y b, y la temperatura en la barra, pueden ser funciones matem\u00e1ticas, por lo que se propone usar el formato lambda de las variables x o y. En el ejercicio b\u00e1sico presentado en la parte te\u00f3rica, la barra tiene temperatura constante, por lo que escribe como una constante mas cero por la variable independiente, esta forma facilita la evaluaci\u00f3n de un vector xi por ejemplo, en caso de disponer una expresi\u00f3n diferente.<\/p>\n\n\n los dem\u00e1s par\u00e1metros del ejercicio se ingresan m\u00e1s adelante de forma semejante al ejercicio con Numpy.<\/p>\n\n\n\n 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\u00f3n 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\u00f3n usando \u03bb como referencia.<\/p>\n\n\n\n Al desarrollar el m\u00e9todo expl\u00edcito se requiere indicar el nodo u(i,j+1<\/strong><\/em>) a buscar<\/p>\n\n\n\n El valor de K se puede determinar en la expresi\u00f3n al reordenar para que el t\u00e9rmino de Dx de mayor grado sea la unidad. El valor de K ser\u00e1 el del t\u00e9rmino Dy con signo negativo, al encontrarse en el lado izquierdo (LHS) de la expresi\u00f3n, en la propuesta te\u00f3rica se encontraba como positivo del lado derecho.<\/p>\n\n\n\n El t\u00e9rmino de mayor grado se puede obtener al revisar la expresi\u00f3n por cada t\u00e9rmino de la suma, buscando el factor que contiene Dx o Dy en cada caso. Por ser un poco ordenado se a\u00f1ade la funci\u00f3n EDP Parab\u00f3lica<\/p>\n\n\n\n cont\u00ednua a discreta<\/a><\/p>\n\n\n\n algoritmo discreta<\/a><\/p>\n\n\n\n algoritmo explicito<\/a><\/p>\n<\/div>\n\n\n\n Para revisar los resultados de algunos pasos para obtener la expresi\u00f3n discreta de EDP, se usa un diccionario a ser usado en el bloque de salida. Para las expresiones de ecuaci\u00f3n se usa EDP Parab\u00f3lica<\/p>\n\n\n\n cont\u00ednua a discreta<\/a><\/p>\n\n\n\n algoritmo discreta<\/a><\/p>\n\n\n\n
\n\n\n\n1. Ejercicio<\/h2>\n\n\n\n
Longitud en x<\/strong> a = 0, b = 1, Constante K<\/strong>= 4
Tama\u00f1o de paso dx = 0.1, dt = dx\/10<\/p>\n\n\n\n
\n\n\n\n
\n\n\n\n2. EDP Parab\u00f3lica cont\u00ednua a discreta<\/h2>\n\n\n\n
EDP Parab\u00f3lica cont\u00ednua a discreta\n ordenDx : 2\n ordenDy : 1\n\n edp=0 :\n 2 \n \u2202 \u2202 \n\u2500\u2500\u2500(u(x, y)) - 4\u22c5\u2500\u2500(u(x, y)) = 0\n 2 \u2202y \n\u2202x \n K_ : 4\n\n Lambda_L :\n Dy \n\u2500\u2500\u2500\u2500\u2500\n 2\n4\u22c5Dx \n Lambda L_k : 0.250000000000000\n\n (discreta=0)*Dy**ordeny :\n 2\u22c5Dy\u22c5U(i, j) Dy\u22c5U(i - 1, j) Dy\u22c5U(i + 1, j)\n4\u22c5U(i, j) - 4\u22c5U(i, j + 1) - \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n 2 2 2 \n Dx Dx Dx \n\n discreta_L=0 :\n\u239b 2\u22c5Dy\u239e Dy\u22c5U(i - 1, j) Dy\u22c5U(i + 1, j) \n\u239c4 - \u2500\u2500\u2500\u2500\u239f\u22c5U(i, j) - 4\u22c5U(i, j + 1) + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 = 0\n\u239c 2 \u239f 2 2 \n\u239d Dx \u23a0 Dx Dx \n\n discreta :\n2\u22c5U(i, j) - 4\u22c5U(i, j + 1) + U(i - 1, j) + U(i + 1, j) = 0\n<\/code><\/pre>\n\n\n\n\n# INGRESO\nfxy = 0*x+0*y # f(x,y) = 0 , ecuacion complementaria\n# ecuacion edp : LHS=RHS\nLHS = sym.diff(u,x,2) \nRHS = 4*sym.diff(u,y,1) + fxy\nedp = sym.Eq(LHS-RHS,0)\n<\/pre><\/div>\n\n\n
\n# centrada, adelante \ndif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))\/(Dx**2),\n sym.diff(u,y,1): (U.subs(j,j+1)-U)\/Dy}\n<\/pre><\/div>\n\n\n\n# Valores de frontera\nfya = lambda y: 60 +0*y # izquierda\nfyb = lambda y: 40 +0*y # derecha\nfxc = lambda x: 25 +0*x # inferior, funci\u00f3n inicial\n<\/pre><\/div>\n\n\n
buscar = U.subs<\/strong><\/em>(j<\/strong><\/em>,j+1) # U[i,j+1<\/strong>] m\u00e9todo expl\u00edcito<\/span><\/code><\/pre>\n\n\n\nedp_coef_Dx()<\/code>.<\/p>\n\n\n\n
\n\n\n\n
\n\n\n\n3. Algoritmo cont\u00ednua a discreta en Python<\/h2>\n\n\n\n
sym.pprint()<\/code>.<\/p>\n\n\n\n# Ecuaciones Diferenciales Parciales Parab\u00f3licas\n# EDP Parab\u00f3licas cont\u00ednua a discreta con Sympy\nimport numpy as np\nimport sympy as sym\n\n# u(x,y) funciones continuas y variables simb\u00f3licas usadas\nt = sym.Symbol('t',real=True) # independiente\nx = sym.Symbol('x',real=True)\ny = sym.Symbol('y',real=True) \nu = sym.Function('u')(x,y) # funcion \nf = sym.Function('f')(x,y) # funcion complemento\nK = sym.Symbol('K',real=True)\n# U[i,j] funciones discretas y variables simb\u00f3licas usadas\ni = sym.Symbol('i',integer=True,positive=True) # indices\nj = sym.Symbol('j',integer=True,positive=True)\nDt = sym.Symbol('Dt',real=True,positive=True)\nDx = sym.Symbol('Dx',real=True,positive=True)\nDy = sym.Symbol('Dy',real=True,positive=True)\nL = sym.Symbol('L',real=True)\nU = sym.Function('U')(i,j)\n\n# INGRESO\nfxy = 0*x+0*y # f(x,y) = 0 , ecuacion complementaria\n# ecuacion edp : LHS=RHS\nLHS = sym.diff(u,x,2) \nRHS = 4*sym.diff(u,y,1) + fxy\nedp = sym.Eq(LHS-RHS,0)\n\n# centrada, adelante \ndif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))\/(Dx**2),\n sym.diff(u,y,1): (U.subs(j,j+1)-U)\/Dy}\n\nbuscar = U.subs(j,j+1) # U[i,j+1] m\u00e9todo expl\u00edcito\n\n# Valores de frontera\nfya = lambda y: 60 +0*y # izquierda\nfyb = lambda y: 40 +0*y # derecha\nfxc = lambda x: 25 +0*x # inferior, funci\u00f3n inicial\n\n# [a,b] dimensiones de la barra\na = 0 # longitud en x\nb = 1\ny0 = 0 # tiempo inicial, aumenta con dt en n iteraciones\n\n# muestreo en ejes, discreto\nx_tramos = 10\ndx = (b-a)\/x_tramos\ndy = dx\/10\nn_iteraciones = 4 # iteraciones en tiempo\n\nverdigitos = 2 # para mostrar en tabla de resultados\ncasicero = 1e-15 # para redondeo de t\u00e9rminos en ecuacion\n\n# PROCEDIMIENTO --------------------------------------\n\ndef edp_coef_Dx(edp,x,ordenx):\n ''' Extrae el coeficiente de la derivada Dx de ordenx,\n edp es la ecuaci\u00f3n como lhs=rhs\n '''\n coeff_x = 1.0 # valor inicial\n # separa cada t\u00e9rmino de suma\n term_suma = sym.Add.make_args(edp.lhs)\n for term_k in term_suma:\n if term_k.is_Mul: # mas de un factor\n factor_Mul = sym.Mul.make_args(term_k)\n coeff_temp = 1; coeff_usar = False\n # separa cada factor de t\u00e9rmino \n for factor_k in factor_Mul:\n if not(factor_k.is_Derivative):\n coeff_temp = coeff_temp*factor_k\n else: # factor con derivada de ordenx\n partes = factor_k.args\n if partes[1]==(x,ordenx):\n coeff_usar = True\n if coeff_usar==True:\n coeff_x = coeff_x*coeff_temp\n return(coeff_x)\n\ndef redondea_coef(ecuacion, precision=6,casicero = 1e-15):\n ''' redondea coeficientes de t\u00e9rminos suma de una ecuacion\n ecuaci\u00f3n como lhs=rhs\n '''\n tipo = type(ecuacion)\n tipo_eq = False\n if tipo == sym.core.relational.Equality:\n RHS = ecuacion.rhs # separa lado derecho\n ecuacion = ecuacion.lhs # analiza lado izquierdo\n tipo = type(ecuacion)\n tipo_eq = True\n\n if tipo == sym.core.add.Add: # t\u00e9rminos suma de ecuacion\n term_sum = sym.Add.make_args(ecuacion)\n ecuacion = sym.S.Zero # vacia\n for term_k in term_sum:\n # factor multiplicativo de termino suma\n term_mul = sym.Mul.make_args(term_k)\n producto = sym.S.One\n for factor in term_mul:\n if not(factor.has(sym.Symbol)): # es numerico\n factor = np.around(float(factor),precision)\n if (abs(factor)%1)<casicero: # si es entero\n factor = int(factor)\n producto = producto*factor\n ecuacion = ecuacion + producto\n \n if tipo == sym.core.mul.Mul: # termino \u00fanico, busca factores\n term_mul = sym.Mul.make_args(ecuacion)\n producto = sym.S.One\n for factor in term_mul:\n if not(factor.has(sym.Symbol)): # es numerico\n factor = np.around(float(factor),precision)\n if (abs(factor)%1)<casicero: # si es entero\n factor = int(factor)\n producto = producto*factor\n ecuacion = producto\n \n if tipo == float: # solo un numero\n if (abs(ecuacion)%1)<casicero: \n ecuacion = int(ecuacion)\n \n if tipo_eq==True: # era igualdad, integra lhs=rhs\n ecuacion = sym.Eq(ecuacion,RHS)\n\n return(ecuacion)\n\n\nresultado = {} # resultados en diccionario\n# orden de derivada por x, y\nedp_x = edp.subs(x,0)\nedp_y = edp.subs(y,0)\nordenDx = sym.ode_order(edp_x,u)\nordenDy = sym.ode_order(edp_y,u)\nresultado['ordenDx'] = ordenDx # guarda en resultados\nresultado['ordenDy'] = ordenDy\n# coeficiente derivada orden mayor a 1 (d2u\/dx2)\ncoeff_x = edp_coef_Dx(edp,x,ordenDx)\nLHS = edp.lhs # lado izquierdo de edp\nRHS = edp.rhs # lado derecho de edp\nif not(coeff_x==1):\n LHS = LHS\/coeff_x\n RHS = RHS\/coeff_x\n# toda la expresi\u00f3n a la izquierda\nedp = sym.Eq(LHS-RHS,0) \nK_ = -edp_coef_Dx(edp,y,ordenDy)\nif abs(K_)%1<casicero: # si es entero\n K_ = int(K_)\nresultado['edp=0'] = edp\nresultado['K_'] = K_\n\n# lambda en ecuacion edp\nlamb = (Dy**ordenDy)\/(Dx**ordenDx)\nif ordenDy==1 and ordenDx==2:\n lamb = lamb\/K_\nresultado['Lambda_L'] = lamb\n# valor de Lambda en ecuacion edp\nL_k = lamb.subs([(Dx,dx),(Dy,dy)])\nif abs(L_k)%1<casicero:\n L_k = int(L_k)\nresultado['Lambda L_k'] = L_k\n\ndiscreta = edp.lhs # EDP discreta\nfor derivada in dif_dividida: # reemplaza diferencia dividida\n discreta = discreta.subs(derivada,dif_dividida[derivada])\n\n# simplifica con Lambda_L\ndiscreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)\nresultado['(discreta=0)*Dy**ordeny'] = discreta_L\n\ndiscreta_L = discreta_L.subs(lamb,L)\ndiscreta_L = sym.collect(discreta_L,U) #agrupa u[,]\ndiscreta_L = sym.Eq(discreta_L,0)\nresultado['discreta_L=0'] = discreta_L\n \n# sustituye constantes en ecuaci\u00f3n a iterar\ndiscreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),(L,L_k)])\ndiscreta_L = redondea_coef(discreta_L)\nresultado['discreta'] = discreta_L\n\n# SALIDA\nnp.set_printoptions(precision=verdigitos)\nalgun_numero = [int,float,str,'Lambda L_k']\nprint('EDP Parab\u00f3lica cont\u00ednua a discreta')\nfor entrada in resultado:\n tipo = type(resultado[entrada])\n if tipo in algun_numero or entrada in algun_numero:\n print('',entrada,':',resultado[entrada])\n else:\n print('\\n',entrada,':')\n sym.pprint(resultado[entrada])\n<\/pre><\/div>\n\n\n
\n\n\n\n