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 ]