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} = 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 numerico 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.
# Ecuaciones Diferenciales Parciales # Elipticas. Método iterativo 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 dx = 0.25 dy = 0.25 maxitera = 100 tolera = 0.0001 # 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) # valores en fronteras u[0,:] = Ta u[n-1,:] = Tb u[:,0] = Tc u[:,m-1] = Td # valor inicial de iteración promedio = (Ta+Tb+Tc+Td)/4 u[1:n-1,1:m-1] = promedio # iterar puntos interiores itera = 0 converge = 0 while not(itera>=maxitera or converge==1): itera = itera +1 nueva = 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 = nueva-u erroru = np.linalg.norm(np.abs(diferencia)) if (erroru<tolera): converge = 1 # SALIDA np.set_printoptions(precision=2) print('converge = ', converge) print('xi=') print(xi) print('yj=') print(yj) print('matriz u') print(u)
Un ejemplo de resultados:
converge = 1 xi= [ 0. 0.25 0.5 0.75 1. 1.25 1.5 1.75 2. ] yi= [ 0. 0.25 0.5 0.75 1. 1.25 1.5 ] matriz u [[ 50. 60. 60. 60. 60. 60. 70. ] [ 50. 55.6 58.23 60. 61.77 64.4 70. ] [ 50. 54.15 57.34 60. 62.66 65.85 70. ] [ 50. 53.67 56.97 60. 63.03 66.33 70. ] [ 50. 53.55 56.87 60. 63.13 66.45 70. ] [ 50. 53.67 56.97 60. 63.03 66.33 70. ] [ 50. 54.15 57.34 60. 62.66 65.85 70. ] [ 50. 55.6 58.23 60. 61.77 64.4 70. ] [ 50. 60. 60. 60. 60. 60. 70. ]] >>>
Cuyos valores se interpretan mejor en una gráfica, en este caso 3D:
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 graficada se obtiene como la transpuesta de u
# Gráfica import matplotlib.pyplot as plt from matplotlib import cm from mpl_toolkits.mplot3d import Axes3D X, Y = np.meshgrid(xi, yj) U = np.transpose(u) # ajuste de índices fila es x figura = plt.figure() ax = Axes3D(figura) ax.plot_surface(X, Y, U, rstride=1, cstride=1, cmap=cm.Reds) plt.title('EDP elíptica') plt.xlabel('x') plt.ylabel('y') plt.show()