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 como:
| 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}=0 \lambda u_{i+1,j} +(-2\lambda-2) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1} +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 # d2u/dx2 + du/dt = f(x,y) import numpy as np # INGRESO fxy = lambda x,y: 0*x+0*y # f(x,y) = 0 , ecuacion de Poisson # Condiciones iniciales en los bordes, valores de frontera Ta = 60 # izquierda de la placa Tb = 25 # derecha de la placa #Tc = 50 # inferior fxc = lambda x: 50+0*x # f(x) inicial inferior # Td = 70 # superior fxd = lambda x: 70+0*x # f(x) 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 # maximo de iteraciones tolera = 0.0001 verdigitos = 2 # decimales a mostrar en tabla de resultados vertabla = True # ver iteraciones # coeficientes de U[x,t]. factores P,Q,R lamb = (dy**2)/(dx**2) P = lamb # izquierda P*U[i-1,j] Q = -2*lamb-2 # centro Q*U[i,j] R = lamb # derecha R*U[i+1,j] # PROCEDIMIENTO # iteraciones en xi,yj ancho,profundidad xi = np.arange(x0,xn+dx/2,dx) yj = np.arange(y0,yn+dy/2,dy) n = len(xi) m = len(yj) # Matriz u[xi,yj], tabla de resultados 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) # Tc inferior u[:,0] = fic fid = fxd(xi) # Td superior u[:,m-1] = fid # estado inicial dentro de la placa # promedio = (Ta+Tb+Tc+Td)/4 promedio = (Ta+Tb+np.max(fic)+np.max(fid))*(1/-Q) u[1:n-1,1:m-1] = promedio np.set_printoptions(precision=verdigitos) if vertabla == True: print('u inicial:') print(np.transpose(u)) # iterar puntos interiores U_3D = [np.copy(u)] U_error = ['--'] 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] = P*u[i-1,j]+R*u[i+1,j]+u[i,j-1]+u[i,j+1] u[i,j] = u[i,j] - fxy(xi[i],yj[j])*(dy**2) u[i,j] = (1/-Q)*u[i,j] diferencia = u - Uantes errado_u = np.max(np.abs(diferencia)) if (errado_u<tolera): converge=True U_3D.append(np.copy(u)) U_error.append(np.around(errado_u,verdigitos)) if (vertabla==True) and (itera<=3): np.set_printoptions(precision=verdigitos) print('itera:',itera-1,' ; ','u:') print(np.transpose(u)) print('errado_u: ', errado_u) if (vertabla==True) and (itera==4): print('... continua') # SALIDA print('Método Iterativo EDP Elíptica') print('iteraciones:',itera,' converge = ', converge) print('errado_u: ', errado_u) 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 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 ]