EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n ejercicio<\/a><\/p>\n\n\n\n M\u00e9todo impl\u00edcito<\/a><\/strong>:<\/p>\n\n\n\n Anal\u00edtico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n<\/div>\n\n\n\n Referencia<\/strong><\/em>: Chapra 30.3 p893 pdf917, Burden 9Ed 12.2 p729, Rodr\u00edguez 10.2.4 p415<\/p>\n\n\n\n Siguiendo el tema anterior en EDP Parab\u00f3licas<\/a>, se tiene que:<\/p>\n\n\n\n \\frac{\\partial ^2 u}{\\partial x ^2} = K\\frac{\\partial u}{\\partial t}<\/span>\n\n\n\n ..<\/p>\n\n\n\n EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n ejercicio<\/a><\/p>\n\n\n\n M\u00e9todo impl\u00edcito<\/a><\/strong>:<\/p>\n\n\n\n Anal\u00edtico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n<\/div>\n\n\n\n En \u00e9ste caso se usan diferencias finitas centradas y hacia atr\u00e1s; la l\u00ednea de referencia es t1<\/sub>:<\/p>\n\n\n \\frac{\\partial^2 u}{\\partial x^2} = \\frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\\Delta x)^2} <\/span>\n\n\n \\frac{\\partial u}{\\partial t} = \\frac{u_{i,j} - u_{i,j-1} }{\\Delta t} <\/span>\n\n\n\n La selecci\u00f3n de las diferencias divididas corresponden a los puntos que se quieren usar para el c\u00e1lculo, se observan mejor en la gr\u00e1fica de la malla.<\/p>\n\n\n\n Luego se sustituyen en la ecuaci\u00f3n del problema, obteniendo:<\/p>\n\n\n \\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\\Delta x)^2} = K\\frac{u_{i,j}-u_{i,j-1}}{\\Delta t}<\/span>\n\n\n\n De la gr\u00e1fica se destaca que en la f\u00f3rmula, dentro del tri\u00e1ngulo solo hay DOS valores desconocidos<\/strong>, destacados por los punto en amarillo.<\/p>\n\n\n\n En la ecuaci\u00f3n se representa por U[i,j] y U[i+1,j]. Por lo que ser\u00e1 necesario crear un sistema de ecuaciones sobre toda la l\u00ednea de tiempo t1<\/sub> para resolver el problema.<\/p>\n\n\n\n Despejando la ecuaci\u00f3n, se agrupan t\u00e9rminos constantes: <\/p>\n\n\n\n \u03bb = \\frac{\\Delta t}{K (\\Delta x)^2} <\/span>\n\n\n \\lambda \\left( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \\right) = u_{i,j}-u_{i,j-1}<\/span>\n\n\n \\lambda u_{i+1,j}-2\\lambda u_{i,j}+ \\lambda u_{i-1,j} = u_{i,j}-u_{i,j-1}<\/span>\n\n\n \\lambda u_{i-1,j} + (-1-2\\lambda) u_{i,j} + \\lambda u_{i+1,j} = -u_{i,j-1} <\/span>\n\n\n\n Los par\u00e1metro P, Q y R se determinan de forma semejante al m\u00e9todo expl\u00edcito:<\/p>\n\n\n P = \\lambda <\/span>\n\n\n Q = -1-2\\lambda <\/span>\n\n\n R = \\lambda <\/span>\n\n\n Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j} = -u_{i,j-1} <\/span>\n\n\n\n EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n ejercicio<\/a><\/p>\n\n\n\n M\u00e9todo impl\u00edcito<\/a><\/strong>:<\/p>\n\n\n\n Anal\u00edtico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n<\/div>\n\n\n\n Los valores en los extremos son conocidos, para los puntos intermedios se crea un sistema de ecuaciones para luego usar la forma Ax=B y resolver los valores para cada u(xi<\/sub>,tj<\/sub>).<\/p>\n\n\n\n Por ejemplo con cuatro tramos entre extremos se tiene que: i=1,j=1<\/p>\n\n\n Pu_{0,1} + Qu_{1,1} + Ru_{2,1} = -u_{1,0} <\/span>\n\n\n P T_{A} + Qu_{1,1} + Ru_{2,1} = -T_{0} <\/span>\n\n\n Qu_{1,1} + Ru_{2,1} = -T_{0} - P T_{A} <\/span>\n\n\n\n i=2,j=1<\/p>\n\n\n Pu_{1,1} + Qu_{2,1} + Ru_{3,1} = -u_{2,0} <\/span>\n\n\n Pu_{1,1} + Qu_{2,1} + Ru_{3,1} = -T_{0} <\/span>\n\n\n\n i=3,j=1<\/p>\n\n\n Pu_{2,1} + Qu_{3,1} + Ru_{4,1} = -u_{3,0} <\/span>\n\n\n Pu_{2,1} + Qu_{3,1} + RT_{B} = -T_{0} <\/span>\n\n\n Pu_{2,1} + Qu_{3,1} = -T_{0} - RT_{B} <\/span>\n\n\n\n agrupando ecuaciones y sustituyendo valores conocidos:<\/p>\n\n\n\n\\begin{cases} Qu_{1,1} + Ru_{2,1} + 0 &= -T_{0}-PT_{A}\\\\Pu_{1,1} + Qu_{2,1} + Ru_{3,1} &= -T_{0}\\\\0+Pu_{2,1}+Qu_{3,1}&=-T_{0}-RT_{B}\\end{cases}<\/span>\n\n\n\n que genera la matriz a resolver:<\/p>\n\n\n\\begin{bmatrix} Q && R && 0 && -T_{0}-PT_{A}\\\\P && Q && R && -T_{0}\\\\0 && P && Q && -T_{0}-RT_{B}\\end{bmatrix}<\/span>\n\n\n\n Use alguno de los m\u00e9todos de la unidad 3 para resolver el sistema y obtener los valores correspondientes.<\/p>\n\n\n\n Por la extensi\u00f3n de la soluci\u00f3n es conveniente usar un algoritmo y convertir los pasos o partes pertinentes a funciones.<\/p>\n\n\n\n Tarea<\/strong>: Revisar y comparar con el m\u00e9todo expl\u00edcito.<\/p>\n\n\n\n EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n ejercicio<\/a><\/p>\n\n\n\n M\u00e9todo impl\u00edcito<\/a><\/strong>:<\/p>\n\n\n\n Anal\u00edtico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n<\/div>\n\n\n\n Para la soluci\u00f3n con el m\u00e9todo impl\u00edcito, se obtiene el mismo resultado en la gr\u00e1fica y tabla. Aunque el algoritmo es diferente.<\/p>\n\n\n\n Instrucciones en Python<\/p>\n\n\n EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n ejercicio<\/a><\/p>\n\n\n\n M\u00e9todo impl\u00edcito<\/a><\/strong>:<\/p>\n\n\n\n Anal\u00edtico<\/a><\/p>\n\n\n\n Algoritmo<\/a><\/p>\n\n\n\n gr\u00e1fica<\/a><\/p>\n<\/div>\n\n\n\n Para realizar la gr\u00e1fica se aplica lo mismo que en el m\u00e9todo expl\u00edcito<\/p>\n\n\n Sin embargo para la gr\u00e1fica en 3D se ajustan los puntos de referencia agregados por las diferencias finitas.<\/p>\n\n\n Queda por revisar la convergencia y estabilidad<\/strong><\/em> de la soluci\u00f3n a partir de los O(h) de cada aproximaci\u00f3n usada. Revisar los criterios en Chapra 30.2.1 p891, Burden 9Ed 12.2 p727, Rodr\u00edguez 10.2.2 p409 .<\/p>\n\n\n\n Tarea o proyecto<\/strong>: Realizar la comparaci\u00f3n de tiempos de ejecuci\u00f3n entre los m\u00e9todos expl\u00edcitos e impl\u00edcitos. La parte de animaci\u00f3n funciona igual en ambos m\u00e9todos. Los tiempos de ejecuci\u00f3n se determinan usando las instrucciones descritas en el enlace: Tiempos de Ejecuci\u00f3n en Python<\/a><\/p>\n\n\n\n
\n\n\n\n1. Ejercicio - EDP Parab\u00f3licas<\/h2>\n\n\n\n
![\"EDP](\"http:\/\/blog.espol.edu.ec\/algoritmos101\/files\/2017\/10\/EdpParabImpli_animado.gif\")
<\/figure>\n\n\n\n
\n\n\n\n
\n\n\n\n2. M\u00e9todo impl\u00edcito - paso a paso<\/h2>\n\n\n\n
![\"Barra](\"http:\/\/blog.espol.edu.ec\/algoritmos101\/files\/2017\/10\/BarraMetalica04.png\")
<\/figure>\n\n\n\n![\"EDP](\"http:\/\/blog.espol.edu.ec\/algoritmos101\/files\/2017\/10\/EDP_MImplicito02.png\")
<\/figure>\n\n\n \\frac{\\Delta t}{K}\\left( \\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\\Delta x)^2} \\right) = \\cancel{K} \\frac{u_{i,j}-u_{i,j-1}}{\\cancel{\\Delta t}}\\frac{\\cancel{\\Delta t}}{\\cancel{K}}<\/span>\n\n\n \\frac{\\Delta t}{K(\\Delta x)^2}\\left( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \\right) = u_{i,j}-u_{i,j-1}<\/span>\n\n\n\n
\n\n\n\n
\n\n\n\n3. Desarrollo anal\u00edtico - paso a paso<\/h2>\n\n\n\n
indice de tiempo es 1 e \u00edndice de x es 1.<\/p>\n\n\n\n
\n\n\n\n
\n\n\n\n4. Algoritmos en Python<\/h2>\n\n\n\n
![\"EDP](\"http:\/\/blog.espol.edu.ec\/algoritmos101\/files\/2017\/10\/EDP_Parabolica01.png\")
<\/figure>\n\n\n\n
algunos valores:<\/p>\n\n\n\niteraci\u00f3n j: 1\nA:\n [[-1.5 0.25 0. 0. 0. 0. 0. 0. 0. ]\n [ 0.25 -1.5 0.25 0. 0. 0. 0. 0. 0. ]\n [ 0. 0.25 -1.5 0.25 0. 0. 0. 0. 0. ]\n [ 0. 0. 0.25 -1.5 0.25 0. 0. 0. 0. ]\n [ 0. 0. 0. 0.25 -1.5 0.25 0. 0. 0. ]\n [ 0. 0. 0. 0. 0.25 -1.5 0.25 0. 0. ]\n [ 0. 0. 0. 0. 0. 0.25 -1.5 0.25 0. ]\n [ 0. 0. 0. 0. 0. 0. 0.25 -1.5 0.25]\n [ 0. 0. 0. 0. 0. 0. 0. 0.25 -1.5 ]]\nB:\n [-40. -25. -25. -25. -25. -25. -25. -25. -35.]\nresultado en j: 1 \n [31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57]\niteraci\u00f3n j: 2\nA:\n [[-1.5 0.25 0. 0. 0. 0. 0. 0. 0. ]\n [ 0.25 -1.5 0.25 0. 0. 0. 0. 0. 0. ]\n [ 0. 0.25 -1.5 0.25 0. 0. 0. 0. 0. ]\n [ 0. 0. 0.25 -1.5 0.25 0. 0. 0. 0. ]\n [ 0. 0. 0. 0.25 -1.5 0.25 0. 0. 0. ]\n [ 0. 0. 0. 0. 0.25 -1.5 0.25 0. 0. ]\n [ 0. 0. 0. 0. 0. 0.25 -1.5 0.25 0. ]\n [ 0. 0. 0. 0. 0. 0. 0.25 -1.5 0.25]\n [ 0. 0. 0. 0. 0. 0. 0. 0.25 -1.5 ]]\nB:\n [-46.01 -26.03 -25.18 -25.03 -25.01 -25.01 -25.08 -25.44 -37.57]\nresultado en j: 2 \n [35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39]\niteraci\u00f3n j: 3\nA:\n [[-1.5 0.25 0. 0. 0. 0. 0. 0. 0. ]\n [ 0.25 -1.5 0.25 0. 0. 0. 0. 0. 0. ]\n [ 0. 0.25 -1.5 0.25 0. 0. 0. 0. 0. ]\n [ 0. 0. 0.25 -1.5 0.25 0. 0. 0. 0. ]\n [ 0. 0. 0. 0.25 -1.5 0.25 0. 0. 0. ]\n [ 0. 0. 0. 0. 0.25 -1.5 0.25 0. 0. ]\n [ 0. 0. 0. 0. 0. 0.25 -1.5 0.25 0. ]\n [ 0. 0. 0. 0. 0. 0. 0.25 -1.5 0.25]\n [ 0. 0. 0. 0. 0. 0. 0. 0.25 -1.5 ]]\nB:\n [-50.25 -27.49 -25.55 -25.12 -25.03 -25.05 -25.24 -26.07 -39.39]\nresultado en j: 3 \n [38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72]\nEDP Parab\u00f3lica - M\u00e9todo impl\u00edcito \nlambda: 0.25\nx: [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]\nt: [0. 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ] ... 1.0\nTabla de resultados en malla EDP Parab\u00f3lica\nj, U[:,: 10 ], primeras iteraciones\n10 [60. 47.43 37.58 31.3 27.98 26.64 26.63 27.83 30.43 34.63 40. ]\n9 [60. 46.75 36.68 30.56 27.48 26.3 26.33 27.47 30.04 34.33 40. ]\n8 [60. 45.96 35.69 29.8 27.01 26. 26.06 27.12 29.6 33.99 40. ]\n7 [60. 45.01 34.6 29.02 26.57 25.73 25.81 26.76 29.13 33.58 40. ]\n6 [60. 43.87 33.39 28.24 26.16 25.51 25.58 26.41 28.6 33.09 40. ]\n5 [60. 42.45 32.06 27.48 25.8 25.32 25.39 26.07 28.03 32.48 40. ]\n4 [60. 40.67 30.61 26.75 25.51 25.18 25.24 25.76 27.41 31.71 40. ]\n3 [60. 38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72 40. ]\n2 [60. 35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39 40. ]\n1 [60. 31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57 40. ]\n0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]\n>>> <\/code><\/pre>\n\n\n\n\n# EDP parab\u00f3licas d2u\/dx2 = K du\/dt\n# m\u00e9todo impl\u00edcito\n# Referencia: Chapra 30.3 p.895 pdf.917\n# Rodriguez 10.2.5 p.417\nimport numpy as np\n\n# INGRESO\n# Valores de frontera\nTa = 60\nTb = 40\n#T0 = 25 # estado inicial de barra\nfx = lambda x: 25 + 0*x # funci\u00f3n inicial para T0\na = 0 # longitud en x\nb = 1\n\nK = 4 # Constante K\ndx = 0.1 # Tama\u00f1o de paso\ndt = dx\/10\n\nn = 100 # iteraciones en tiempo\nverdigitos = 2 # decimales a mostrar en tabla de resultados\n\n# PROCEDIMIENTO\n# iteraciones en longitud\nxi = np.arange(a,b+dx\/2,dx)\nfi = fx(xi)\nm = len(xi)\nultimox = m-1\nultimot = n-1\n# Resultados en tabla u[x,t]\nu = np.zeros(shape=(m,n), dtype=float)\n\n# valores iniciales de u[:,j]\nj = 0\nu[0,:]= Ta\nu[1:ultimox,j] = fi[1:ultimox]\nu[ultimox,:] = Tb\n\n# factores P,Q,R\nlamb = dt\/(K*dx**2)\nP = lamb\nQ = -1 -2*lamb\nR = lamb\n\n# Calcula U para cada tiempo + dt\ntj = np.arange(0,(n+1)*dt,dt)\nj = 1\nwhile not(j>=n):\n # Matriz de ecuaciones\n k = m-2\n A = np.zeros(shape=(k,k), dtype = float)\n B = np.zeros(k, dtype = float)\n for f in range(0,k,1):\n if (f>0):\n A[f,f-1]=P\n A[f,f] = Q\n if (f<(k-1)):\n A[f,f+1]=R\n B[f] = -u[f+1,j-1]\n B[0] = B[0]-P*u[0,j]\n B[k-1] = B[k-1]-R*u[m-1,j]\n # Resuelve sistema de ecuaciones\n C = np.linalg.solve(A, B)\n\n # copia resultados a u[i,j]\n for f in range(0,k,1):\n u[f+1,j] = C[f]\n\n # siguiente iteraci\u00f3n\n j = j + 1\n\n # muestra 3 primeras iteraciones\n np.set_printoptions(precision=verdigitos)\n if j<(3+2): \n print('iteraci\u00f3n j:',j-1)\n print('A:\\n',A)\n print('B:\\n',B)\n print('resultado en j:',j-1,'\\n',C)\n \n# SALIDA\nprint('EDP Parab\u00f3lica - M\u00e9todo impl\u00edcito ')\nprint('lambda: ',np.around(lamb,verdigitos))\nprint('x:',xi)\nprint('t:',tj[0:len(xi)],'...',tj[-1])\nprint('Tabla de resultados en malla EDP Parab\u00f3lica')\nprint('j, U[:,:',ultimox,'], primeras iteraciones')\nfor j in range(ultimox,-1,-1):\n print(j,u[:,j])\n<\/pre><\/div>\n\n\n
\n\n\n\n
\n\n\n\n5. Gr\u00e1fica con Python<\/h2>\n\n\n\n
\n# GRAFICA ------------\nimport matplotlib.pyplot as plt\ntramos = 10\nsalto = int(n\/tramos) # evita muchas l\u00edneas\nif (salto == 0):\n salto = 1\nfor j in range(0,n,salto):\n vector = u[:,j]\n plt.plot(xi,vector)\n plt.plot(xi,vector, '.',color='red')\nplt.xlabel('x[i]')\nplt.ylabel('u[j]')\nplt.title('Soluci\u00f3n EDP parab\u00f3lica')\nplt.show()\n<\/pre><\/div>\n\n\n\n# GRAFICA en 3D\n# vector de tiempo, simplificando en tramos\ntramos = 10\nsalto = int(n\/tramos)\nif (salto == 0):\n salto = 1\ntj = np.arange(0,n*dt,dt)\ntk = np.zeros(tramos,dtype=float)\n\n# Extrae parte de la matriz U,acorde a los tramos\nU = np.zeros(shape=(tramos,m),dtype=float)\nfor k in range(0,tramos,1):\n U[k,:] = u[:,k*salto]\n tk[k] = tj[k*salto]\n# Malla para cada eje X,Y\nXi, Yi = np.meshgrid(xi,tk)\n\nfig_3D = plt.figure()\ngraf_3D = fig_3D.add_subplot(111, projection='3d')\ngraf_3D.plot_wireframe(Xi,Yi,U,\n color ='blue',label='EDP Parab\u00f3lica')\ngraf_3D.plot(Xi[2,0],Yi[2,0],U[2,0],'o',color ='orange')\ngraf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')\ngraf_3D.plot(Xi[2,2],Yi[2,2],U[2,2],'o',color ='salmon')\ngraf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='green')\ngraf_3D.set_title('EDP Parab\u00f3lica')\ngraf_3D.set_xlabel('x')\ngraf_3D.set_ylabel('t')\ngraf_3D.set_zlabel('U')\ngraf_3D.legend()\ngraf_3D.view_init(35, -45)\nplt.show()\n<\/pre><\/div>\n\n\n
\n\n\n\n