Ejercicio: 2Eva_IIT2019_T3 EDP elíptica, placa en (1,1)
dada la ecuación del problema:
\frac{\delta ^2 u}{\delta x^2} + \frac{\delta ^2 u}{\delta y^2} = \frac{x}{y} + \frac{y}{x}
1 < x < 2
1 < y < 2
Se convierte a la versión discreta usando diferencias divididas centradas:
\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} = \frac{x_i}{y_j} + \frac{y_j}{x_i}
Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy2
\frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) +
+ \frac{\Delta y^2}{\Delta y^2}\Big(u[i,j-1]-2u[i,j]+u[i,j+1]\Big) =
=\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)
los tamaños de paso en ambos ejes son de igual valor, se simplifica la ecuación
\lambda= \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\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)
u[i-1,j]-4u[i,j]+u[i+1,j] +
+ u[i,j-1]+u[i,j+1] =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)
Iteraciones
que permite plantear las ecuaciones para cada punto en posición [i,j]
i=1, j=1
u[0,1]-4u[1,1]+u[2,1] +
+ u[1,0]+u[1,2] =(0.25)^2\Big( \frac{0.25}{0.25} + \frac{0.25}{0.25}\Big)
0.25ln(0.25)-4u[1,1]+u[2,1] +
+ 0.25ln(0.25)+u[1,2] = 0.125
-4u[1,1]+u[2,1] +u[1,2] = 0.125 - 2(0.25)ln(0.25)
-4u[1,1]+u[2,1] +u[1,2] = 0.8181
i=2, j=1
u[1,1]-4u[2,1]+u[3,1]+
+0.5 ln(0.5)+u[2,2]=0.15625
u[1,1]-4u[2,1]+u[3,1]+u[2,2]=0.15625-0.5 ln(0.5)
u[1,1]-4u[2,1]+u[3,1]+u[2,2]=0.8493
i=3, j=1
u[2,1]-4u[3,1]+u[4,1] +
+ u[3,0]+u[3,2] =(0.25)^2\Big( \frac{0.75}{0.25} + \frac{0.25}{0.75}\Big)
u[2,1]-4u[3,1]+u[4,1]+u[3,2] = 0.20833333-0.75 ln(0.75)
u[2,1]-4u[3,1]+u[4,1]+u[3,2] =0.4240
Tarea: continuar con el ejercicio hasta plantear todo el sistema de ecuaciones.
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(
[0.125 - 2(0.25)ln(0.25),
0.15625 - 0.5ln(0.5),
0.20833 - 0.75 ln(0.75),
...])
B = [0.8181,
0.8493,
0.4240,
...]
Algoritmo con Python
Con valores para la matriz solución:
iteraciones: 15
error entre iteraciones: 6.772297286980838e-05
solución para u:
[[0. 0.2789294 0.6081976 0.9793276 1.3862943]
[0.2789294 0.6978116 1.1792239 1.7127402 2.2907268]
[0.6081976 1.1792239 1.8252746 2.5338403 3.2958368]
[0.9793276 1.7127402 2.5338403 3.4280053 4.3846703]
[1.3862943 2.2907268 3.2958368 4.3846703 5.5451774]]
>>>
y algoritmo detallado:
# 2Eva_IIT2019_T2 EDO, problema de valor inicial
# método iterativo
import numpy as np
# INGRESO
# longitud en x
a = 1
b = 2
# longitud en y
c = 1
d = 2
# tamaño de paso
dx = 0.25
dy = 0.25
# funciones en los bordes de la placa
abajo = lambda x,y: x*np.log(x)
arriba = lambda x,y: x*np.log(4*(x**2))
izquierda = lambda x,y: y*np.log(y)
derecha = lambda x,y: 2*y*np.log(2*y)
# función de la ecuación
fxy = lambda x,y: x/y + y/x
# control de iteraciones
maxitera = 100
tolera = 0.0001
# PROCEDIMIENTO
# tamaño de la matriz
n = int((b-a)/dx)+1
m = int((d-c)/dy)+1
# vectores con valore de ejes
xi = np.linspace(a,b,n)
yj = np.linspace(c,d,m)
# matriz de puntos muestra
u = np.zeros(shape=(n,m),dtype=float)
# valores en los bordes
u[:,0] = abajo(xi,yj[0])
u[:,m-1] = arriba(xi,yj[m-1])
u[0,:] = izquierda(xi[0],yj)
u[n-1,:] = derecha(xi[n-1],yj)
# valores interiores la mitad en intervalo x,
# mitad en intervalo y, para menos iteraciones
mx = int(n/2)
my = int(m/2)
promedio = (u[mx,0]+u[mx,m-1]+u[0,my]+u[n-1,my])/4
u[1:n-1,1:m-1] = promedio
# método iterativo
itera = 0
converge = 0
while not(itera>=maxitera or converge==1):
itera = itera +1
# copia u para calcular errores entre iteraciones
nueva = np.copy(u)
for i in range(1,n-1):
for j in range(1,m-1):
# usar fórmula desarrollada para algoritmo
fij = (dy**2)*fxy(xi[i],yj[j])
u[i,j]=(u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-fij)/4
diferencia = nueva-u
erroru = np.linalg.norm(np.abs(diferencia))
if (erroru<tolera):
converge=1
# SALIDA
print('iteraciones: ',itera)
print('error entre iteraciones: ',erroru)
print('solución para u: ')
print(u)
# Gráfica
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
# matrices de ejes para la gráfica 3D
X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()
grafica = Axes3D(figura)
grafica.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()
