s2Eva_IT2018_T3 EDP Eliptica

Generar las ecuaciones a resolver usando diferencias finitas divididas centradas:

\frac{\delta ^2 u}{\delta x^2} + \frac{\delta ^2 u}{\delta y^2} = 2(x^2+y^2)

por facilidad se sustituye también en la forma discreta de la ecuación:

f[i,j] = f(x_i,y_j) = 2(x_{i} ^2+y_{j} ^2)
\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} = f[i,j]
\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] = \Delta y^2 f[i,j]

dado que hx = hy = 1/3

\frac{\Delta y^2}{\Delta x^2} = 1
(u[i-1,j]-2u[i,j]+u[i+1,j]) + + u[i,j-1]-2u[i,j]+u[i,j+1] = = \Delta y^2 f[i,j]
u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] = \Delta y^2 f[i,j]
4u[i,j] = u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1]-\Delta y^2 f[i,j]
u[i,j] = \frac{1}{4} \Big(u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1]-\Delta y^2 f[i,j] \Big)

La última ecuación puede ser usada de forma iterativa, para lo cual hay que definir los valores iniciales de la matriz u.

Al conocer el rango de operación para los ejes x, y, hx, hy se realizan los cálculos para:

1. Evaluar los valores para cada eje x[i], y[j]

x[i]:
[ 0.    0.33  0.67  1.  ]
y[j]:
[ 0.    0.33  0.67  1.  ]

2. Evaluar en cada punto generando una matriz f(i,j):

f[i,j]:
[[ 0.    0.22  0.89  2.  ]
 [ 0.22  0.44  1.11  2.22]
 [ 0.89  1.11  1.78  2.89]
 [ 2.    2.22  2.89  4.  ]]

3. Se evaluan las funciones indicadas para la frontera y se tiene la matriz inicial para u:

matriz inicial u[i,j]:
[[ 1.    1.33  1.67  2.  ]
 [ 1.33  0.    0.    2.44]
 [ 1.67  0.    0.    3.11]
 [ 2.    2.44  3.11  4.  ]]

con lo que se puede trabajar cada punto i,j de forma iterativa, teniendo como resultado para la matriz u:

resultado para u, iterando: 
converge =  1
[[ 1.    1.33  1.67  2.  ]
 [ 1.33  1.68  2.05  2.44]
 [ 1.67  2.05  2.53  3.11]
 [ 2.    2.44  3.11  4.  ]]

La gráfica usando una mayor resolución para tener una idea de la solución:


Los resultados se obtienen usando las siguientes instrucciones:

# 2da Evaluación I Término 2018
# Tema 3. EDP Eliptica
import numpy as np

# INGRESO
# ejes x,y
x0 = 0 ; xn = 1 ; hx = (1/3)# (1/3)/10
y0 = 0 ; yn = 1 ; hy = (1/3) # (1/3)/10
# Fronteras
fux0 = lambda x: x+1
fu0y = lambda y: y+1
fux1 = lambda x: x**2 + x + 2
fu1y = lambda y: y**2 + y + 2

fxy = lambda x,y: 2*(x**2+y**2)

# PROCEDIMIENTO
xi = np.arange(x0,xn+hx,hx)
yj = np.arange(y0,yn+hy,hy)
n = len(xi)
m = len(yj)
# funcion f[xi,yi]
fij = np.zeros(shape=(n,m), dtype = float)
for i in range(0,n,1):
    for j in range(0,m,1):
        fij[i,j]=fxy(xi[i],yj[j])
# matriz inicial u[i,j]
u = np.zeros(shape=(n,m), dtype = float)
u[:,0] = fux0(xi)
u[0,:] = fu0y(yj)
u[:,m-1] = fux1(xi)
u[n-1,:] = fu1y(yj)

uinicial = u.copy()

# Calcular de forma iterativa
maxitera = 100
tolera = 0.0001
# valor inicial de iteración
promedio = (np.max(u)+np.min(u))/2
u[1:n-1,1:m-1] = promedio
# iterar puntos interiores
itera = 0
converge = 0
erroru = 2*tolera # para calcular al menos una matriz
while not(erroru=maxitera):
    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]-(hy**2)*fij[i,j])/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
if (erroru<tolera):
    converge=1

# SALIDA
np.set_printoptions(precision=2)
print('x[i]:')
print(xi)
print('y[j]:')
print(yj)
print('f[i,j]:')
print(fij)
print('matriz inicial u[i,j]:')
print(uinicial)
print('resultado para u, iterando: ')
print('converge = ', converge)
print('iteraciones = ', itera)
print(u)

para obtener la gráfica se debe añadir:

# Gráfica
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

X, Y = np.meshgrid(xi, yj)
figura = plt.figure()
ax = Axes3D(figura)
U = np.transpose(u) # ajuste de índices fila es x
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()