{"id":1005,"date":"2017-10-03T09:05:33","date_gmt":"2017-10-03T14:05:33","guid":{"rendered":"http:\/\/blog.espol.edu.ec\/matg1013\/?p=1005"},"modified":"2026-03-14T19:51:41","modified_gmt":"2026-03-15T00:51:41","slug":"edp-parabolicas-metodo-explicito","status":"publish","type":"post","link":"https:\/\/blog.espol.edu.ec\/algoritmos101\/mn-u07\/edp-parabolicas-metodo-explicito\/","title":{"rendered":"7.1.1 EDP Parab\u00f3licas m\u00e9todo expl\u00edcito con Python"},"content":{"rendered":"\n
\n\n\n\n
\n

EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n

ejercicio<\/a><\/p>\n\n\n\n

M\u00e9todo expl\u00edcito<\/a>:<\/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

\n\n\n\n

1. Ejercicio - EDP Parab\u00f3lica<\/h2>\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\n

Siguiendo el tema anterior en EDP Parab\u00f3licas<\/a>, para resolver la parte num\u00e9rica considere como:<\/p>\n\n\n\n

Valores de frontera: T<\/strong>a = 60, T<\/strong>b = 40, T<\/strong>0 = 25
Longitud en x<\/strong>: a = 0, b = 1,  Constante K<\/strong>= 4
Tama\u00f1o de paso: dx = 0.1, dt = dx\/10<\/code><\/pre>\n\n\n\n
Para la ecuaci\u00f3n diferencial parcial parab\u00f3lica:<\/p>\n\n\n\n \\frac{\\partial ^2 u}{\\partial x ^2} = K\\frac{\\partial u}{\\partial t}<\/span>\n\n\n\n
\"EDP<\/figure>\n\n\n\n
\n\n\n\n
\n
EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n
ejercicio<\/a><\/p>\n\n\n\n
M\u00e9todo expl\u00edcito<\/a>:<\/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
\n\n\n\n
2. M\u00e9todo Expl\u00edcito - paso a paso<\/h2>\n\n\n\n
Se usan las diferencias divididas, donde se requieren dos valores para la derivada de orden uno y tres valores para la derivada de orden dos:<\/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+1} - u_{i,j} }{\\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. La l\u00ednea de referencia es el tiempo t0<\/sub><\/p>\n\n\n\n
\"\"<\/figure>\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+1}-u_{i,j}}{\\Delta t}<\/span>\n\n\n\n
De la gr\u00e1fica se destaca que en la f\u00f3rmula solo hay un valor desconocido, destacado por el punto en amarillo dentro del tri\u00e1ngulo. En la ecuaci\u00f3n se representa por U[i,j+1].<\/p>\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+1}-u_{i,j}}{\\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+1}-u_{i,j}<\/span>\n\n\n\n
Despejando la ecuaci\u00f3n, se agrupan t\u00e9rminos constantes:<\/p>\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+1}-u_{i,j}<\/span>\n\n\n \\lambda u_{i+1,j} - 2 \\lambda u_{i,j}+ \\lambda u_{i-1,j} = u_{i,j+1}-u_{i,j}<\/span>\n\n\n\n
quedando la ecuaci\u00f3n, con los t\u00e9rminos ordenados de izquierda a derecha como en la gr\u00e1fica:<\/p>\n\n\n u_{i,j+1} = \\lambda u_{i-1,j} +(1-2\\lambda)u_{i,j}+\\lambda u_{i+1,j}<\/span>\n\n\n\n
Al resolver se encuentra que que cada valor en un punto amarillo se calcula como una suma ponderada de valores conocidos, por lo que el desarrollo se conoce como el m\u00e9todo expl\u00edcito<\/strong>.<\/p>\n\n\n\n
La ponderaci\u00f3n est\u00e1 determinada por los t\u00e9rminos con los factores P, Q, y R.<\/p>\n\n\n\n
\n P = \\lambda <\/span>\n\n\n\n
;<\/p>\n\n\n\n
 Q = 1-2\\lambda <\/span><\/p>\n\n\n\n
;<\/p>\n\n\n\n
 R = \\lambda <\/span><\/p>\n<\/div>\n\n\n u_{i ,j+1} = Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j}<\/span>\n\n\n\n
Para el ejercicio, el valor num\u00e9rico de lambda es:<\/p>\n\n\n \u03bb = \\frac{\\Delta t}{K (\\Delta x)^2} = \\frac{0.1\/10}{4 (0.1)^2} = \\frac{1}{4} <\/span>\n\n\n\n
Con lo que P = 1\/4, Q = 1\/2 y R=1\/4.<\/p>\n\n\n u_{i ,j+1} = \\frac{u_{i-1,j}}{4} + \\frac{u_{i,j}}{2} + \\frac{u_{i+1,j}}{4}<\/span>\n\n\n\n
\n\n\n\n
\n
EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n
ejercicio<\/a><\/p>\n\n\n\n
M\u00e9todo expl\u00edcito<\/a>:<\/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
\n\n\n\n
3. Desarrollo anal\u00edtico - paso a paso<\/h2>\n\n\n\n
Se muestran tres iteraciones para j=0, que corresponden a las tres primeras casillas de la matriz U(i,j).<\/p>\n\n\n\n
j=0 ; i=1<\/strong><\/p>\n\n\n\n
          U(0, 0)   U(1, 0)   U(2, 0)\nU(1, 1) = \u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\n             4         2         4   \n           60     25     25\nU(1, 1) = \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\n           4       2      4  \nU(1, 1) = 33.75<\/code><\/pre>\n\n\n\n
j=0 ; i=2<\/strong><\/p>\n\n\n\n
j=0 ; i=2\n          U(1, 0)   U(2, 0)   U(3, 0)\nU(2, 1) = \u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\n             4         2         4   \n           25     25     25\nU(2, 1) = \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\n            4      2      4 \nU(2, 1) = 25.0<\/code><\/pre>\n\n\n\n
j=0 ; i=3<\/strong><\/p>\n\n\n\n
          U(2, 0)   U(3, 0)   U(4, 0)\nU(3, 1) = \u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\u2500\u2500\u2500\n             4         2         4   \n           25     25     25\nU(3, 1) = \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500 + \u2500\u2500\u2500\u2500\n            4      2      4    \nU(3, 1) = 25.0<\/code><\/pre>\n\n\n\n
F\u00f3rmulas que se desarrollan usando un algoritmo y considerando que al disminuir los valores de \u0394x y \u0394t la cantidad de operaciones aumenta.<\/p>\n\n\n\n
\"EDP<\/figure>\n\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.<\/p>\n\n\n\n
Revisar los criterios en<\/strong><\/em>: Chapra 30.2.1 p891 pdf915, Burden 9Ed 12.2 p727, Rodr\u00edguez 10.2.2 pdf409 .<\/p>\n\n\n \\lambda \\leq \\frac{1}{2} <\/span>\n\n\n\n
Cuando \u03bb \u2264 1\/2 se tiene como resultado una soluci\u00f3n donde los errores no crecen, sino que oscilan.<\/p>\n\n\n\n
Haciendo \u03bb \u2264 1\/4 asegura que la soluci\u00f3n no oscilar\u00e1.<\/p>\n\n\n\n
Tambi\u00e9n se sabe que con \u03bb= 1\/6 se tiende a minimizar los errores por truncamiento (v\u00e9ase Carnahan y cols., 1969).<\/p>\n\n\n\n
Para resolver la parte num\u00e9rica considere que:<\/p>\n\n\n\n
Valores de frontera: T<\/strong>a = 60, T<\/strong>b = 40, T<\/strong>0 = 25
Longitud en x<\/strong>: a = 0, b = 1,  Constante K<\/strong>= 4
Tama\u00f1o de paso: dx = 0.1, dt = dx\/10<\/code><\/pre>\n\n\n\n
\n\n\n\n
\n
EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n
ejercicio<\/a><\/p>\n\n\n\n
M\u00e9todo expl\u00edcito<\/a>:<\/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
\n\n\n\n
4. Algoritmo en Python para M\u00e9todo expl\u00edcito<\/h2>\n\n\n\n
Se muestra una tabla resumida de resultados a forma de ejemplo.<\/p>\n\n\n\n
M\u00e9todo expl\u00edcito EDP Parab\u00f3lica\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] ... 1.0\nTabla de resultados en malla EDP Parab\u00f3lica\nj, U[:,: 5 ], primeras iteraciones de  100\n5 [60.   44.21 32.93 27.29 25.41 25.05 25.18 25.98 28.4  33.23 40.  ]\n4 [60.   42.77 31.29 26.37 25.14 25.   25.06 25.59 27.7  32.62 40.  ]\n3 [60.   40.86 29.38 25.55 25.   25.   25.   25.23 26.88 31.8  40.  ]\n2 [60.   38.12 27.19 25.   25.   25.   25.   25.   25.94 30.62 40.  ]\n1 [60.   33.75 25.   25.   25.   25.   25.   25.   25.   28.75 40.  ]\n0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]\n>>> <\/code><\/pre>\n\n\n\n
Se presenta la propuesta del algoritmo EDP parab\u00f3licas para el m\u00e9todo expl\u00edcito.<\/p>\n\n\n
\n# EDP parab\u00f3licas d2u\/dx2  = K du\/dt\n# m\u00e9todo expl\u00edcito,usando diferencias divididas\nimport numpy as np\n\n# INGRESO\n# Valores de frontera\nTa = 60  # izquierda de barra\nTb = 40  # derecha de barra\n#Tc = 25 # estado inicial de barra en x\nfxc = lambda x: 25 + 0*x # f(x) en tiempo inicial\n\n# dimensiones de la barra\na = 0  # longitud en x\nb = 1\nt0 = 0 # tiempo inicial, aumenta con dt en n iteraciones\n\nK = 4     # Constante K\ndx = 0.1  # muestreo en x, tama\u00f1o de paso\ndt = dx\/10\n\nn = 100  # iteraciones en tiempo\nverdigitos = 2   # decimales a mostrar en tabla de resultados\n\n# coeficientes de U[x,t]. factores P,Q,R, \nlamb = dt\/(K*dx**2)\nP = lamb       # izquierda P*U[i-1,j]\nQ = 1 - 2*lamb # centro Q*U[i,j]\nR = lamb       # derecha R*U[i+1,j]\n\n# PROCEDIMIENTO\n# iteraciones en x, longitud\nxi = np.arange(a,b+dx\/2,dx)\nm = len(xi) \nultimox = m-1\nultimot = n-1\n\n# u[xi,tj], tabla de resultados\nu = np.zeros(shape=(m,n+1), dtype=float)\n\n# u[i,j], valores iniciales\nu[0,:] = Ta  # Izquierda\nu[ultimox,:] = Tb  # derecha\n# estado inicial de barra en x, Tc\nfic = fxc(xi) # f(x) en tiempo inicial\nu[1:ultimox,0] = fic[1:ultimox]\n\n# Calcula U para cada tiempo + dt\ntj = np.arange(t0,(n+1)*dt,dt)\nfor j  in range(0,n,1):\n    for i in range(1,ultimox,1):\n        # ecuacion discreta, entre [1,ultimox]\n        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]\n\n# SALIDA\nj_mostrar = 5\nnp.set_printoptions(precision=verdigitos)\nprint('M\u00e9todo expl\u00edcito EDP Parab\u00f3lica')\nprint('lambda: ',np.around(lamb,verdigitos))\nprint('x:',xi)\nprint('t:',tj[0:j_mostrar],'...',tj[-1])\nprint('Tabla de resultados en malla EDP Parab\u00f3lica')\nprint('j, U[:,:',j_mostrar,'], primeras iteraciones de ',n)\nfor j in range(j_mostrar,-1,-1):\n    print(j,u[:,j])\n<\/pre><\/div>\n\n\n
\n\n\n\n
\n
EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n
ejercicio<\/a><\/p>\n\n\n\n
M\u00e9todo expl\u00edcito<\/a>:<\/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
\n\n\n\n
5. Gr\u00e1fica con Python para M\u00e9todo expl\u00edcito<\/h2>\n\n\n\n
Si la cantidad de puntos aumenta al disminuir \u0394x y \u0394t, es mejor disminuir la cantidad de curvas a usar en el gr\u00e1fico para evitar superponerlas. Se usa el par\u00e1metro 'salto' para indicar cada cu\u00e1ntas curvas calculadas se incorporan en la gr\u00e1fica.<\/p>\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
Note que en la gr\u00e1fica se toma como coordenadas x vs t, mientras que para la soluci\u00f3n de la malla en la matriz las se usan filas y columnas. En la matriz el primer \u00edndice es fila y el segundo \u00edndice es columna.<\/p>\n\n\n\n
En el caso de una gr\u00e1fica que se realice usando los ejes x,t,U en tres dimensiones se requiere usar las instrucciones:<\/p>\n\n\n
\n# GRAFICA en 3D ------\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=(m,tramos),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)\nU = np.transpose(U)\n\nfig_3D = plt.figure()\ngraf_3D = fig_3D.add_subplot(111, projection='3d')\ngraf_3D.plot_wireframe(Xi,Yi,U, color ='blue')\ngraf_3D.plot(xi[0],tk[1],U[1,0],'o',color ='orange')\ngraf_3D.plot(xi[1],tk[1],U[1,1],'o',color ='green')\ngraf_3D.plot(xi[2],tk[1],U[1,2],'o',color ='green')\ngraf_3D.plot(xi[1],tk[2],U[2,1],'o',color ='salmon',\n             label='U[i,j+1]')\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
\n
EDP Parab\u00f3lica<\/a><\/p>\n\n\n\n
ejercicio<\/a><\/p>\n\n\n\n
M\u00e9todo expl\u00edcito<\/a>:<\/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
\n\n\n\n
Solving the heat equation | DE3. 3Blue1Brown. 16-junio-2019<\/p>\n\n\n\n
\n