U4videos Interpolación polinómica
U3videos Sistemas de ecuaciones
U2videos Raíces de ecuaciones en una variable
U1videos Polinomio de Taylor con Python
7.3 EDP hiperbólicas con Python
Referencia: Chapra PT8.1 p860, Rodríguez 10.4 p435
Las Ecuaciones Diferenciales Parciales tipo hiperbólicas semejantes a la mostrada, para un ejemplo en particular, representa la ecuación de onda de una dimensión u[x,t], que representa el desplazamiento vertical de una cuerda.
\frac{\partial ^2 u}{\partial t^2}=c^2\frac{\partial ^2 u}{\partial x^2}Los extremos de la cuerda de longitud unitaria están sujetos y referenciados a una posición 0 a la izquierda y 1 a la derecha.
u[x,t] , 0<x<1, t≥0
u(0,t) = 0 , t≥0
u(1,t) = 0 , t≥0
Al inicio, la cuerda está estirada por su punto central:
u(x,0) = \begin{cases} -0.5x &, 0\lt x\leq 0.5 \\ 0.5(x-1) &, 0.5\lt x \lt 1 \end{cases}Se suelta la cuerda, con velocidad cero desde la posición inicial:
\frac{\partial u(x,0)}{\partial t} = 0La solución se realiza de forma semejante al procedimiento para EDP parabólicas y elípticas. Se cambia a la forma discreta usando diferencias finitas divididas:
\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta t)^2} =c^2 \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}El error es del orden O(\Delta x)^2 + O(\Delta t)^2.
se reagrupa de la forma:
El término constante, por facilidad se simplifica con el valor de 1
\lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} =1si c = 2 y Δx = 0.2, se deduce que Δt = 0.1
que al sustituir en la ecuación, se simplifica anulando el término u[i,j]:
u_{i,j+1}+u_{i,j-1} = u_{i+1,j}+u_{i-1,j}en los que intervienen solo los puntos extremos. Despejando el punto superior del rombo:
u_{i,j+1} = u_{i+1,j}-u_{i,j-1}+u_{i-1,j}Convergencia:
\lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \leq 1para simplificar aún más el problema, se usa la condición velocidad inicial de la cuerda igual a cero
\frac{\delta u_{i,0}}{\delta t}=\frac{u_{i,1}-u_{i,-1}}{2\Delta t} = 0 u_{i,-1}=u_{i,1}que se usa para cuando el tiempo es cero, permite calcular los puntos para t[1]:
j=0
u_{i,1} = u_{i+1,0}-u_{i,-1}+u_{i-1,0} u_{i,1} = u_{i+1,0}-u_{i,1}+u_{i-1,0} 2 u_{i,1} = u_{i+1,0}+u_{i-1,0} u_{i,1} = \frac{u_{i+1,0}+u_{i-1,0}}{2}Aplicando solo cuando j = 0
que al ponerlos en la malla se logra un sistema explícito para cada u[i,j], con lo que se puede resolver con un algoritmo:
Algoritmo en Python:
# Ecuaciones Diferenciales Parciales # Hiperbólica. Método explicito import numpy as np def cuerdainicio(xi): n = len(xi) y = np.zeros(n,dtype=float) for i in range(0,n,1): if (xi[i]<=0.5): y[i]=-0.5*xi[i] else: y[i]=0.5*(xi[i]-1) return(y) # INGRESO x0 = 0 xn = 1 # Longitud de cuerda y0 = 0 yn = 0 # Puntos de amarre t0 = 0 c = 2 # discretiza tramosx = 16 tramost = 32 dx = (xn-x0)/tramosx dt = dx/c # PROCEDIMIENTO xi = np.arange(x0,xn+dx,dx) tj = np.arange(0,tramost*dt+dt,dt) n = len(xi) m = len(tj) u = np.zeros(shape=(n,m),dtype=float) u[:,0] = cuerdainicio(xi) u[0,:] = y0 u[n-1,:] = yn # Aplicando condición inicial j = 0 for i in range(1,n-1,1): u[i,j+1] = (u[i+1,j]+u[i-1,j])/2 # Para los otros puntos for j in range(1,m-1,1): for i in range(1,n-1,1): u[i,j+1] = u[i+1,j]-u[i,j-1]+u[i-1,j] # SALIDA np.set_printoptions(precision=2) print('xi =') print(xi) print('tj =') print(tj) print('matriz u =') print(u)
con lo que se obtienen los resultados numéricos, que para mejor interpretación se presentan en una gráfica estática y otra animada.
# GRAFICA import matplotlib.pyplot as plt for j in range(0,m,1): y = u[:,j] plt.plot(xi,y) plt.title('EDP hiperbólica') plt.xlabel('x') plt.ylabel('y') plt.show() # **** GRÁFICO CON ANIMACION *********** import matplotlib.animation as animation # Inicializa parametros de trama/foto retardo = 70 # milisegundos entre tramas tramas = m maximoy = np.max(np.abs(u)) figura, ejes = plt.subplots() plt.xlim([x0,xn]) plt.ylim([-maximoy,maximoy]) # lineas de función y polinomio en gráfico linea_poli, = ejes.plot(xi,u[:,0], '-') plt.axhline(0, color='k') # Eje en x=0 plt.title('EDP hiperbólica') # plt.legend() # txt_x = (x0+xn)/2 txt_y = maximoy*(1-0.1) texto = ejes.text(x0,txt_y,'tiempo:', horizontalalignment='left') plt.xlabel('x') plt.ylabel('y') plt.grid() # Nueva Trama def unatrama(i,xi,u): # actualiza cada linea linea_poli.set_ydata(u[:,i]) linea_poli.set_xdata(xi) linea_poli.set_label('tiempo linea: '+str(i)) texto.set_text('tiempo['+str(i)+']') # color de la línea if (i<=9): lineacolor = 'C'+str(i) else: numcolor = i%10 lineacolor = 'C'+str(numcolor) linea_poli.set_color(lineacolor) return linea_poli, texto # Limpia Trama anterior def limpiatrama(): linea_poli.set_ydata(np.ma.array(xi, mask=True)) linea_poli.set_label('') texto.set_text('') return linea_poli, texto # Trama contador i = np.arange(0,tramas,1) ani = animation.FuncAnimation(figura, unatrama, i , fargs=(xi,u), init_func=limpiatrama, interval=retardo, blit=True) # Graba Archivo video y GIFAnimado : # ani.save('EDP_hiperbólica.mp4') ani.save('EDP_hiperbolica.gif', writer='imagemagick') plt.draw() plt.show()
Una vez realizado el algoritmo, se pueden cambiar las condiciones iniciales de la cuerda y se observan los resultados.
Se recomienda realizar otros ejercicios de la sección de evaluaciones para EDP Hiperbólicas y observar los resultados con el algoritmo modificado.
7.2.2 EDP Elípticas método implícito con Python
EDP Elípticas [ concepto ] Método implícito: [ Analítico ] [ Algoritmo ]
1. EDP Elípticas: Método Implícito – Desarrollo Analítico
con el resultado desarrollado en EDP elípticas para:
\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1
se puede plantear que:
u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0con lo que para el método implícito, se plantea un sistema de ecuaciones para determinar los valores en cada punto desconocido.
j=1, i =1
u_{2,1}-4u_{1,1}+u_{0,1} + u_{1,2} +u_{1,0} = 0 u_{2,1}-4u_{1,1}+Ta + u_{1,2} +Tc= 0 -4u_{1,1}+u_{2,1}+u_{1,2} = -(Tc+Ta)j=1, i =2
u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +u_{2,0} = 0 u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +Tc = 0 u_{1,1}-4u_{2,1}+u_{3,1}+ u_{2,2}= -Tcj=1, i=3
u_{4,1}-4u_{3,1}+u_{2,1} + u_{3,2} +u_{3,0} = 0 Tb-4u_{3,1}+u_{2,1} + u_{3,2} +Tc = 0 u_{2,1} -4u_{3,1} + u_{3,2} = -(Tc+Tb)j=2, i=1
u_{2,2}-4u_{1,2}+u_{0,2} + u_{1,3} +u_{1,1} = 0 u_{2,2}-4u_{1,2}+Ta + u_{1,3} +u_{1,1} = 0 -4u_{1,2}+u_{2,2}+u_{1,1}+u_{1,3} = -Taj = 2, i = 2
u_{1,2}-4u_{2,2}+u_{3,2} + u_{2,3} +u_{2,1} = 0j = 2, i = 3
u_{4,2}-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0 Tb-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0 u_{2,2} -4u_{3,2}+ u_{3,3} +u_{3,1} = -Tbj=3, i = 1
u_{2,3}-4u_{1,3}+u_{0,3} + u_{1,4} +u_{1,2} = 0 u_{2,3}-4u_{1,3}+Ta + Td +u_{1,2} = 0 -4u_{1,3}+u_{2,3}+u_{1,2} = -(Td+Ta)j=3, i = 2
u_{3,3}-4u_{2,3}+u_{1,3} + u_{2,4} +u_{2,2} = 0 u_{3,3}-4u_{2,3}+u_{1,3} + Td +u_{2,2} = 0 +u_{1,3} -4u_{2,3}+u_{3,3} +u_{2,2} = -Tdj=3, i=3
u_{4,3}-4u_{3,3}+u_{2,3} + u_{3,4} +u_{3,2} = 0 Tb-4u_{3,3}+u_{2,3} + Td +u_{3,2} = 0 u_{2,3}-4u_{3,3}+u_{3,2} = -(Td+Tb)con las ecuaciones se arma una matriz:
A = np.array([ [-4, 1, 0, 1, 0, 0, 0, 0, 0], [ 1,-4, 1, 0, 1, 0, 0, 0, 0], [ 0, 1,-4, 0, 0, 1, 0, 0, 0], [ 1, 0, 0,-4, 1, 0, 1, 0, 0], [ 0, 1, 0, 1,-4, 1, 0, 1, 0], [ 0, 0, 1, 0, 1,-4, 0, 0, 1], [ 0, 0, 0, 1, 0, 0,-4, 1, 0], [ 0, 0, 0, 0, 1, 0, 1,-4, 1], [ 0, 0, 0, 0, 0, 1, 0, 1,-4], ]) B = np.array([-(Tc+Ta),-Tc,-(Tc+Tb), -Ta, 0, -Tb, -(Td+Ta),-Td,-(Td+Tb)])
que al resolver el sistema de ecuaciones se obtiene:
>>> Xu array([ 56.43, 55.71, 56.43, 60. , 60. , 60. , 63.57, 64.29, 63.57])
ingresando los resultados a la matriz u:
xi= [ 0. 0.5 1. 1.5 2. ] yj= [ 0. 0.38 0.75 1.12 1.5 ] matriz u [[ 60. 60. 60. 60. 60. ] [ 50. 56.43 60. 63.57 70. ] [ 50. 55.71 60. 64.29 70. ] [ 50. 56.43 60. 63.57 70. ] [ 60. 60. 60. 60. 60. ]] >>>
EDP Elípticas [ concepto ] Método implícito: [ Analítico ] [ Algoritmo ]
1. EDP Elípticas: Método Implícito – Desarrollo Analítico
Instrucciones en Python
# Ecuaciones Diferenciales Parciales # Elipticas. Método implícito import numpy as np # INGRESO # Condiciones iniciales en los bordes Ta = 60 Tb = 60 Tc = 50 Td = 70 # dimensiones de la placa x0 = 0 xn = 2 y0 = 0 yn = 1.5 # discretiza, supone dx=dy tramosx = 4 tramosy = 4 dx = (xn-x0)/tramosx dy = (yn-y0)/tramosy maxitera = 100 tolera = 0.0001 A = np.array([ [-4, 1, 0, 1, 0, 0, 0, 0, 0], [ 1,-4, 1, 0, 1, 0, 0, 0, 0], [ 0, 1,-4, 0, 0, 1, 0, 0, 0], [ 1, 0, 0,-4, 1, 0, 1, 0, 0], [ 0, 1, 0, 1,-4, 1, 0, 1, 0], [ 0, 0, 1, 0, 1,-4, 0, 0, 1], [ 0, 0, 0, 1, 0, 0,-4, 1, 0], [ 0, 0, 0, 0, 1, 0, 1,-4, 1], [ 0, 0, 0, 0, 0, 1, 0, 1,-4], ]) B = np.array([-(Tc+Ta),-Tc,-(Tc+Tb), -Ta,0,-Tb, -(Td+Ta),-Td,-(Td+Tb)]) # PROCEDIMIENTO # Resuelve sistema ecuaciones Xu = np.linalg.solve(A,B) [nx,mx] = np.shape(A) xi = np.arange(x0,xn+dx,dx) yj = np.arange(y0,yn+dy,dy) n = len(xi) m = len(yj) u = np.zeros(shape=(n,m),dtype=float) u[:,0] = Tc u[:,m-1] = Td u[0,:] = Ta u[n-1,:] = Tb u[1:1+3,1] = Xu[0:0+3] u[1:1+3,2] = Xu[3:3+3] u[1:1+3,3] = Xu[6:6+3] # SALIDA np.set_printoptions(precision=2) print('xi=') print(xi) print('yj=') print(yj) print('matriz u') print(u)
La gráfica de resultados se obtiene de forma semejante al ejercicio con método iterativo.
Se podría estandarizar un poco más el proceso para que sea realizado por el algoritmo y sea más sencillo generar la matriz con más puntos. Tarea.
7.2.1 EDP Elípticas método iterativo con Python
EDP Elípticas [ concepto ] [ ejercicio ]
Método iterativo: [ Analítico ] [ Algoritmo ]
1. Ejercicio
Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694
Siguiendo el tema anterior en EDP Elípticas, para resolver la parte numérica asuma
Valores de frontera: Ta = 60, Tb = 25, Tc = 50, Td = 70
Longitud en x0 = 0, xn = 2, y0 = 0, yn = 1.5
Tamaño de paso dx = 0.25, dy = dx
iteramax=100, tolera = 0.0001
Para la ecuación diferencial parcial Elíptica
\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0
EDP Elípticas [ concepto ] [ ejercicio ]
Método iterativo: [ Analítico ] [ Algoritmo ]
..
2. EDP Elípticas: Método iterativo – Desarrollo Analítico
A partir del resultado desarrollado en EDP elípticas:
\frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0y con el supuesto que: \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1
se puede plantear que:
u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0que reordenando para un punto central desconocido se convierte a:
u_{i,j} = \frac{1}{4} \big[ u_{i+1,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} \big]con lo que se interpreta que cada punto central es el resultado del promedio de los puntos alrededor del rombo formado en la malla.
El cálculo numérico se puede realizar de forma iterativa haciendo varias pasadas en la matriz, promediando cada punto. Para revisar las iteraciones se controla la convergencia junto con un máximo de iteraciones.
EDP Elípticas [ concepto ] [ ejercicio ]
Método iterativo: [ Analítico ] [ Algoritmo ]
3. Algoritmo en Python
Para el ejercicio dado, se presenta el resultado para tres iteraciones y el resultado final con un gráfico en 3D:
u inicial: [[50. 50. 50. 50. 50. 50. 50. 50. 50. ] [60. 51.25 51.25 51.25 51.25 51.25 51.25 51.25 25. ] [60. 51.25 51.25 51.25 51.25 51.25 51.25 51.25 25. ] [60. 51.25 51.25 51.25 51.25 51.25 51.25 51.25 25. ] [60. 51.25 51.25 51.25 51.25 51.25 51.25 51.25 25. ] [60. 51.25 51.25 51.25 51.25 51.25 51.25 51.25 25. ] [70. 70. 70. 70. 70. 70. 70. 70. 70. ]] itera: 0 ; u: [[50. 50. 50. 50. 50. 50. 50. 50. 50. ] [60. 53.12 51.41 50.98 50.87 50.84 50.84 44.27 25. ] [60. 53.91 51.95 51.36 51.18 51.13 51.12 42.91 25. ] [60. 54.1 52.14 51.5 51.3 51.23 51.21 42.59 25. ] [60. 54.15 52.2 51.55 51.34 51.27 51.24 42.52 25. ] [60. 58.85 58.07 57.72 57.58 57.52 57.5 48.76 25. ] [70. 70. 70. 70. 70. 70. 70. 70. 70. ]] errado_u: 8.728094100952148 itera: 1 ; u: [[50. 50. 50. 50. 50. 50. 50. 50. 50. ] [60. 53.83 51.69 50.98 50.75 50.68 49.02 41.73 25. ] [60. 54.97 52.54 51.55 51.18 51.05 48.55 39.47 25. ] [60. 55.31 52.89 51.82 51.39 51.23 48.4 38.85 25. ] [60. 56.59 54.78 53.91 53.54 53.38 50.45 40.76 25. ] [60. 61.17 60.91 60.6 60.42 60.33 57.38 48.29 25. ] [70. 70. 70. 70. 70. 70. 70. 70. 70. ]] errado_u: 3.7443900108337402 itera: 2 ; u: [[50. 50. 50. 50. 50. 50. 50. 50. 50. ] [60. 54.17 51.92 51.06 50.73 50.2 47.62 40.52 25. ] [60. 55.5 52.97 51.76 51.23 50.3 46.45 37.7 25. ] [60. 56.25 53.95 52.75 52.19 51.07 46.71 37.54 25. ] [60. 58.05 56.71 55.9 55.47 54.33 49.8 40.16 25. ] [60. 62.24 62.39 62.18 61.99 60.93 57.25 48.1 25. ] [70. 70. 70. 70. 70. 70. 70. 70. 70. ]] errado_u: 2.0990657806396484 ... continua Método Iterativo EDP Elíptica iteraciones: 41 converge = True xi: [0. 0.25 0.5 0.75 1. 1.25 1.5 1.75 2. ] yj: [0. 0.25 0.5 0.75 1. 1.25 1.5 ] Tabla de resultados en malla EDP Elíptica j, U[i,j] 6 [70. 70. 70. 70. 70. 70. 70. 70. 70.] 5 [60. 64.02 64.97 64.71 63.62 61.44 57.16 47.96 25. ] 4 [60. 61.1 61.14 60.25 58.35 54.98 49.23 39.67 25. ] 3 [60. 59.23 58.25 56.8 54.53 50.89 45.13 36.48 25. ] 2 [60. 57.56 55.82 54.19 52.09 48.92 43.91 36.14 25. ] 1 [60. 55.21 53.27 52.05 50.73 48.78 45.46 39.15 25. ] 0 [50. 50. 50. 50. 50. 50. 50. 50. 50.] >>>
cuyos valores se interpretan mejor en una gráfica, en este caso 3D:
Instrucciones en Python
# Ecuaciones Diferenciales Parciales # Elipticas. Método iterativo import numpy as np # INGRESO # Condiciones iniciales en los bordes Ta = 60 # izquierda Tb = 25 # derecha #Tc = 50 # inferior fxc = lambda x: 50+0*x # función inicial inferior # Td = 70 # superior fxd = lambda x: 70+0*x # función inicial superior # dimensiones de la placa x0 = 0 # longitud en x xn = 2 y0 = 0 # longitud en y yn = 1.5 # discretiza, supone dx=dy dx = 0.25 # Tamaño de paso dy = dx iteramax = 100 tolera = 0.0001 verdigitos = 2 # decimales a mostrar en tabla de resultados # PROCEDIMIENTO xi = np.arange(x0,xn+dx,dx) yj = np.arange(y0,yn+dy,dy) n = len(xi) m = len(yj) # Matriz u u = np.zeros(shape=(n,m),dtype = float) # llena u con valores en fronteras u[0,:] = Ta # izquierda u[n-1,:] = Tb # derecha fic = fxc(xi) u[:,0] = fic # Tc inferior fid = fxd(xi) u[:,m-1] = fid # Td superior # valor inicial de iteración dentro de u # promedio = (Ta+Tb+Tc+Td)/4 promedio = (Ta+Tb+np.max(fic)+np.max(fid))/4 u[1:n-1,1:m-1] = promedio np.set_printoptions(precision=verdigitos) print('u inicial:') print(np.transpose(u)) # iterar puntos interiores itera = 0 converge = False while (itera<=iteramax and converge==False): itera = itera +1 Uantes = np.copy(u) for i in range(1,n-1): for j in range(1,m-1): u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1])/4 diferencia = u - Uantes errado_u = np.max(np.abs(diferencia)) if (errado_u<tolera): converge=True if itera<=3: np.set_printoptions(precision=verdigitos) print('itera:',itera-1,' ; ','u:') print(np.transpose(u)) print('errado_u: ', errado_u) if itera==4: print('... continua') # SALIDA print('Método Iterativo EDP Elíptica') print('iteraciones:',itera,' converge = ', converge) print('xi:', xi) print('yj:', yj) print('Tabla de resultados en malla EDP Elíptica') print('j, U[i,j]') for j in range(m-1,-1,-1): print(j,u[:,j])
La gráfica de resultados requiere ajuste de ejes, pues el índice de filas es el eje x, y las columnas es el eje y. La matriz y sus datos en la gráfica se obtiene como la transpuesta de u
# GRAFICA en 3D ------ import matplotlib.pyplot as plt # Malla para cada eje X,Y Xi, Yi = np.meshgrid(xi, yj) U = np.transpose(u) # ajuste de índices fila es x fig_3D = plt.figure() graf_3D = fig_3D.add_subplot(111, projection='3d') graf_3D.plot_wireframe(Xi,Yi,U, color ='blue') graf_3D.plot(Xi[1,0],Yi[1,0],U[1,0],'o',color ='orange') graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='red', label='U[i,j]') graf_3D.plot(Xi[1,2],Yi[1,2],U[1,2],'o',color ='salmon') graf_3D.plot(Xi[0,1],Yi[0,1],U[0,1],'o',color ='green') graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon') graf_3D.set_title('EDP elíptica') graf_3D.set_xlabel('x') graf_3D.set_ylabel('y') graf_3D.set_zlabel('U') graf_3D.legend() graf_3D.view_init(35, -45) plt.show()
EDP Elípticas [ concepto ] [ ejercicio ]
Método iterativo: [ Analítico ] [ Algoritmo ]
7.2 EDP Elípticas
EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]
1. EDP Elípticas
Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694
Las Ecuaciones Diferenciales Parciales tipo elípticas semejantes a la mostrada:
\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0(ecuación de Laplace, Ecuación de Poisson con f(x,y)=0)
se interpreta como una placa metálica de dimensiones Lx y Ly, delgada con aislante que recubren las caras de la placa, y sometidas a condiciones en las fronteras:
Lx = dimensión x de placa metálica Ly = dimensión y de placa metálica u[0,y] = Ta u[Lx,y] = Tb u[x,0] = Tc u[x,Ly] = Td
Para el planteamiento se usa una malla en la que cada nodo corresponden a los valores u[xi,yj]. Para simplificar la nomenclatura se usan los subíndices i para el eje de las x y j para el eje t, quedando u[i,j].
Se discretiza la ecuación usando diferencias divididas que se sustituyen en la ecuación, por ejemplo:
\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta y)^2}=0Se agrupan los términos de los diferenciales:
\frac{(\Delta y)^2}{(\Delta x)^2} \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0 \lambda \Big( u_{i+1,j}-2u_{i,j} +u_{i-1,j} \Big)+ u_{i,j+1}-2u_{i,j}+u_{i,j-1}=0con lo que se simplifican los valores como uno solo \lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1 . Por facilidad de lo que se realiza se supone que lambda tiene valor de 1 o los delta son iguales.
u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0Obteniendo así la solución numérica conceptual. La forma de resolver el problema determina el nombre del método a seguir.
EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]
7.1.2 EDP Parabólicas método implícito con Python
EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]
1. Ejercicio
Referencia: Chapra 30.3 p893 pdf917, Burden 9Ed 12.2 p729, Rodríguez 10.2.4 p415
Siguiendo el tema anterior en EDP Parabólicas, se tiene que:
\frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]
2. Desarrollo Analítico
En éste caso se usan diferencias finitas centradas y hacia atrás; la línea de referencia es t1:
\frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2} \frac{\partial u}{\partial t} = \frac{u_{i,j} - u_{i,j-1} }{\Delta t}La selección de las diferencias divididas corresponden a los puntos que se quieren usar para el cálculo, se observan mejor en la gráfica de la malla.
Luego se sustituyen en la ecuación del problema, obteniendo:
\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}De la gráfica se destaca que en la fórmula, dentro del triángulo solo hay DOS valores desconocidos, destacados por los punto en amarillo.
En la ecuación se representa por U[i,j] y U[i+1,j]. Por lo que será necesario crear un sistema de ecuaciones sobre toda la línea de tiempo t1 para resolver el problema.
Despejando la ecuación, se agrupan términos constantes: λ = \frac{\Delta t}{K (\Delta x)^2} .
\lambda u_{i-1,j} + (-1-2\lambda) u_{i,j} + \lambda u_{i+1,j} = -u_{i,j-1}Los parámetro P, Q y R se determinan de forma semejante al método explícito:
P = \lambda Q = -1-2\lambda R = \lambda Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j} = -u_{i,j-1}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,tj).
Por ejemplo con cuatro tramos entre extremos se tiene que:
indice de tiempo es 1 e índice de x es 1.
i=1,j=1
Pu_{0,1} + Qu_{1,1} + Ru_{2,1} = -u_{1,0}i=2,j=1
Pu_{1,1} + Qu_{2,1} + Ru_{3,1} = -u_{2,0}i=3,j=1
Pu_{2,1} + Qu_{3,1} + Ru_{4,1} = -u_{3,0}agrupando ecuaciones y sustituyendo valores conocidos:
\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}
que genera la matriz a resolver:
\begin{bmatrix} Q && R && 0 && -T_{0}-PT_{A}\\P && Q && R && -T_{0}\\0 && P && Q && -T_{0}-RT_{B}\end{bmatrix}Use alguno de los métodos de la unidad 3 para resolver el sistema y obtener los valores correspondientes.
Por la extensión de la solución es conveniente usar un algoritmo y convertir los pasos o partes pertinentes a funciones.
Tarea: Revisar y comparar con el método explícito.
EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]
Algoritmos en Python
Para la solución con el método implícito, se obtiene el mismo resultado en la gráfica y tabla. Aunque el algoritmo es diferente.
iteración j: 1 A: [[-1.5 0.25 0. 0. 0. 0. 0. 0. 0. ] [ 0.25 -1.5 0.25 0. 0. 0. 0. 0. 0. ] [ 0. 0.25 -1.5 0.25 0. 0. 0. 0. 0. ] [ 0. 0. 0.25 -1.5 0.25 0. 0. 0. 0. ] [ 0. 0. 0. 0.25 -1.5 0.25 0. 0. 0. ] [ 0. 0. 0. 0. 0.25 -1.5 0.25 0. 0. ] [ 0. 0. 0. 0. 0. 0.25 -1.5 0.25 0. ] [ 0. 0. 0. 0. 0. 0. 0.25 -1.5 0.25] [ 0. 0. 0. 0. 0. 0. 0. 0.25 -1.5 ]] B: [-40. -25. -25. -25. -25. -25. -25. -25. -35.] resultado en j: 1 [31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57] iteración j: 2 A: [[-1.5 0.25 0. 0. 0. 0. 0. 0. 0. ] [ 0.25 -1.5 0.25 0. 0. 0. 0. 0. 0. ] [ 0. 0.25 -1.5 0.25 0. 0. 0. 0. 0. ] [ 0. 0. 0.25 -1.5 0.25 0. 0. 0. 0. ] [ 0. 0. 0. 0.25 -1.5 0.25 0. 0. 0. ] [ 0. 0. 0. 0. 0.25 -1.5 0.25 0. 0. ] [ 0. 0. 0. 0. 0. 0.25 -1.5 0.25 0. ] [ 0. 0. 0. 0. 0. 0. 0.25 -1.5 0.25] [ 0. 0. 0. 0. 0. 0. 0. 0.25 -1.5 ]] B: [-46.01 -26.03 -25.18 -25.03 -25.01 -25.01 -25.08 -25.44 -37.57] resultado en j: 2 [35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39] iteración j: 3 A: [[-1.5 0.25 0. 0. 0. 0. 0. 0. 0. ] [ 0.25 -1.5 0.25 0. 0. 0. 0. 0. 0. ] [ 0. 0.25 -1.5 0.25 0. 0. 0. 0. 0. ] [ 0. 0. 0.25 -1.5 0.25 0. 0. 0. 0. ] [ 0. 0. 0. 0.25 -1.5 0.25 0. 0. 0. ] [ 0. 0. 0. 0. 0.25 -1.5 0.25 0. 0. ] [ 0. 0. 0. 0. 0. 0.25 -1.5 0.25 0. ] [ 0. 0. 0. 0. 0. 0. 0.25 -1.5 0.25] [ 0. 0. 0. 0. 0. 0. 0. 0.25 -1.5 ]] B: [-50.25 -27.49 -25.55 -25.12 -25.03 -25.05 -25.24 -26.07 -39.39] resultado en j: 3 [38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72] EDP Parabólica - Método implícito lambda: 0.25 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 0.05 0.06 0.07 0.08 0.09 0.1 ] ... 1.0 Tabla de resultados en malla EDP Parabólica j, U[:,: 10 ], primeras iteraciones 10 [60. 47.43 37.58 31.3 27.98 26.64 26.63 27.83 30.43 34.63 40. ] 9 [60. 46.75 36.68 30.56 27.48 26.3 26.33 27.47 30.04 34.33 40. ] 8 [60. 45.96 35.69 29.8 27.01 26. 26.06 27.12 29.6 33.99 40. ] 7 [60. 45.01 34.6 29.02 26.57 25.73 25.81 26.76 29.13 33.58 40. ] 6 [60. 43.87 33.39 28.24 26.16 25.51 25.58 26.41 28.6 33.09 40. ] 5 [60. 42.45 32.06 27.48 25.8 25.32 25.39 26.07 28.03 32.48 40. ] 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
# EDP parabólicas d2u/dx2 = K du/dt # método implícito # Referencia: Chapra 30.3 p.895 pdf.917 # Rodriguez 10.2.5 p.417 import numpy as np # INGRESO # Valores de frontera Ta = 60 Tb = 40 #T0 = 25 # estado inicial de barra fx = lambda x: 25 + 0*x # función inicial para T0 a = 0 # longitud en x b = 1 K = 4 # Constante K dx = 0.1 # Tamaño de paso dt = dx/10 n = 100 # iteraciones en tiempo verdigitos = 2 # decimales a mostrar en tabla de resultados # PROCEDIMIENTO # iteraciones en longitud xi = np.arange(a,b+dx,dx) fi = fx(xi) m = len(xi) ultimox = m-1 ultimot = n-1 # Resultados en tabla u[x,t] u = np.zeros(shape=(m,n), dtype=float) # valores iniciales de u[:,j] j = 0 u[0,:]= Ta u[1:ultimox,j] = fi[1:ultimox] u[ultimox,:] = Tb # factores P,Q,R lamb = dt/(K*dx**2) P = lamb Q = -1 -2*lamb R = lamb # Calcula U para cada tiempo + dt tj = np.arange(0,(n+1)*dt,dt) j = 1 while not(j>=n): # Matriz de ecuaciones k = m-2 A = np.zeros(shape=(k,k), dtype = float) B = np.zeros(k, dtype = float) for f in range(0,k,1): if (f>0): A[f,f-1]=P A[f,f] = Q if (f<(k-1)): A[f,f+1]=R B[f] = -u[f+1,j-1] B[0] = B[0]-P*u[0,j] B[k-1] = B[k-1]-R*u[m-1,j] # Resuelve sistema de ecuaciones C = np.linalg.solve(A, B) # copia resultados a u[i,j] for f in range(0,k,1): u[f+1,j] = C[f] # siguiente iteración j = j + 1 # muestra 3 primeras iteraciones np.set_printoptions(precision=verdigitos) if j<(3+2): print('iteración j:',j-1) print('A:\n',A) print('B:\n',B) print('resultado en j:',j-1,'\n',C) # SALIDA print('EDP Parabólica - Método implícito ') print('lambda: ',np.around(lamb,verdigitos)) print('x:',xi) print('t:',tj[0:len(xi)],'...',tj[-1]) print('Tabla de resultados en malla EDP Parabólica') print('j, U[:,:',ultimox,'], primeras iteraciones') for j in range(ultimox,-1,-1): print(j,u[:,j])
Para realizar la gráfica se aplica lo mismo que en el método explícito
# GRAFICA ------------ import matplotlib.pyplot as plt tramos = 10 salto = int(n/tramos) # evita muchas líneas if (salto == 0): salto = 1 for j in range(0,n,salto): vector = u[:,j] plt.plot(xi,vector) plt.plot(xi,vector, '.',color='red') plt.xlabel('x[i]') plt.ylabel('u[j]') plt.title('Solución EDP parabólica') plt.show()
Sin embargo para la gráfica en 3D se ajustan los puntos de referencia agregados por las diferencias finitas.
# GRAFICA en 3D # vector de tiempo, simplificando en tramos tramos = 10 salto = int(n/tramos) if (salto == 0): salto = 1 tj = np.arange(0,n*dt,dt) tk = np.zeros(tramos,dtype=float) # Extrae parte de la matriz U,acorde a los tramos U = np.zeros(shape=(tramos,m),dtype=float) for k in range(0,tramos,1): U[k,:] = u[:,k*salto] tk[k] = tj[k*salto] # Malla para cada eje X,Y Xi, Yi = np.meshgrid(xi,tk) fig_3D = plt.figure() graf_3D = fig_3D.add_subplot(111, projection='3d') graf_3D.plot_wireframe(Xi,Yi,U, color ='blue',label='EDP Parabólica') graf_3D.plot(Xi[2,0],Yi[2,0],U[2,0],'o',color ='orange') graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon') graf_3D.plot(Xi[2,2],Yi[2,2],U[2,2],'o',color ='salmon') graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='green') graf_3D.set_title('EDP Parabólica') graf_3D.set_xlabel('x') graf_3D.set_ylabel('t') graf_3D.set_zlabel('U') graf_3D.legend() graf_3D.view_init(35, -45) plt.show()
Queda por revisar la convergencia y estabilidad de la solución a partir de los O(h) de cada aproximación usada. Revisar los criterios en Chapra 30.2.1 p891, Burden 9Ed 12.2 p727, Rodríguez 10.2.2 p409 .
Tarea o proyecto: Realizar la comparación de tiempos de ejecución entre los métodos explícitos e implícitos. La parte de animación funciona igual en ambos métodos. Los tiempos de ejecución se determinan usando las instrucciones descritas en el enlace: Tiempos de Ejecución en Python
EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]