Unidad 7 EDP Ecuaciones Diferenciales Parciales

2023-10-112025-09-09

7.4 EDP con gráficos animados en Python

Solo para fines didácticos, y como complemento para los ejercicios presentados en la unidad para Ecuaciones Diferenciales Parciales, se presentan las instrucciones para las animaciones usadas en la presentación de los conceptos y ejercicios. Los algoritmos para animación NO son necesarios para realizar los ejercicios, que requieren una parte analítica con al menos tres iteraciones en papel y lápiz. Se lo adjunta como una herramienta didáctica de asistencia para las clases.

EDP Parabólica [ concepto ] [ ejercicio ] || explícito [ algoritmo ] gráfica [ 2D ] [ 3D ] animada [ 2D ] [ 3D ] || implícito gráfica [ 3D animada ] ||
EDP Elípticas [ concepto ] [ ejercicio ] iterativo [ algoritmo ] [ 3D ] [ 3D animada ]
..

1. EDP Parabólicas método explícito

resultados del algoritmo

Método explícito EDP Parabólica
lambda:  0.25
x: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0.   0.01 0.02 0.03 0.04] ... 1.0
Tabla de resultados en malla EDP Parabólica
j, U[:,: 5 ], primeras iteraciones de  100
5 [60.   44.21 32.93 27.29 25.41 25.05 25.18 25.98 28.4  33.23 40.  ]
4 [60.   42.77 31.29 26.37 25.14 25.   25.06 25.59 27.7  32.62 40.  ]
3 [60.   40.86 29.38 25.55 25.   25.   25.   25.23 26.88 31.8  40.  ]
2 [60.   38.12 27.19 25.   25.   25.   25.   25.   25.94 30.62 40.  ]
1 [60.   33.75 25.   25.   25.   25.   25.   25.   25.   28.75 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]

Algoritmo en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito,usando diferencias divididas
import numpy as np

# INGRESO
# Valores de frontera
Ta = 60  # izquierda de barra
Tb = 40  # derecha de barra
#Tc = 25 # estado inicial de barra en x
fxc = lambda x: 25 + 0*x # f(x) en tiempo inicial

# dimensiones de la barra
a = 0  # longitud en x
b = 1
t0 = 0 # tiempo inicial, aumenta con dt en n iteraciones

K = 4     # Constante K
dx = 0.1  # muestreo en x, tamaño de paso
dt = dx/10

n = 100  # iteraciones en tiempo
verdigitos = 2   # decimales a mostrar en tabla de resultados

# coeficientes de U[x,t]. factores P,Q,R, 
lamb = dt/(K*dx**2)
P = lamb       # izquierda P*U[i-1,j]
Q = 1 - 2*lamb # centro Q*U[i,j]
R = lamb       # derecha R*U[i+1,j]

# PROCEDIMIENTO
# iteraciones en x, longitud
xi = np.arange(a,b+dx/2,dx)
m = len(xi) 
ultimox = m-1
ultimot = n-1

# u[xi,tj], tabla de resultados
u = np.zeros(shape=(m,n+1), dtype=float)

# u[i,j], valores iniciales
u[0,:] = Ta  # Izquierda
u[ultimox,:] = Tb  # derecha
# estado inicial de barra en x, Tc
fic = fxc(xi) # f(x) en tiempo inicial
u[1:ultimox,0] = fic[1:ultimox]

# Calcula U para cada tiempo + dt
tj = np.arange(t0,(n+1)*dt,dt)
for j  in range(0,n,1):
    for i in range(1,ultimox,1):
        # ecuacion discreta, entre [1,ultimox]
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]

# SALIDA
j_mostrar = 5
np.set_printoptions(precision=verdigitos)
print('Método explícito EDP Parabólica')
print('lambda: ',np.around(lamb,verdigitos))
print('x:',xi)
print('t:',tj[0:j_mostrar],'...',tj[-1])
print('Tabla de resultados en malla EDP Parabólica')
print('j, U[:,:',j_mostrar,'], primeras iteraciones de ',n)
for j in range(j_mostrar,-1,-1):
    print(j,u[:,j])

1.1 Gráficas 2D, EDP Parabólica, método explícito

# GRAFICA 2D ------------
import matplotlib.pyplot as plt
tramos = 10
salto = int(n/tramos) # evita muchas líneas
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.',color='red')
plt.xlabel('x[i]')
plt.ylabel('u[j]')
plt.title('Solución EDP parabólica')
#plt.show()

1.3 Gráficas 3D, EDP Parabólica, método explícito

Esta sección es la base para la generación de las animaciones siguientes, por lo que debe ser incluida antes de las dos siguientes animaciones 2D o 3D. Al menos hasta la parte donde se realiza la transpuesta de U.

# GRAFICA 3D -----------
tj = np.arange(0,n*dt,dt)
tk = np.zeros(tramos,dtype=float)

# Extrae parte de la matriz U,acorde a los tramos
U = np.zeros(shape=(m,tramos),dtype=float)
for k in range(0,tramos,1):
    U[:,k] = u[:,k*salto]
    tk[k] = tj[k*salto]
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi,tk)
U = np.transpose(U)

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[0],tk[1],U[1,0],'o',color ='orange')
graf_3D.plot(xi[1],tk[1],U[1,1],'o',color ='green')
graf_3D.plot(xi[2],tk[1],U[1,2],'o',color ='green')
graf_3D.plot(xi[1],tk[2],U[2,1],'o',color ='salmon',
             label='U[i,j+1]')
graf_3D.set_title('EDP Parabólica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('t')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()

1.4 Gráficas 2D animada, EDP Parabólica, método explícito

La animación en 2D para la EDP parabólica método explícito se añade, luego de comentar el #plt.show() anterior

# GRAFICA 2D ANIMADA --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation
unmetodo = 'EDP Parabólica - Método explícito'
narchivo = 'EdpParabolica' # nombre archivo GIF

# Parametros de trama/foto
retardo = 500   # milisegundos entre tramas

# GRAFICA animada en fig_ani
fig_ani, graf_ani = plt.subplots()

# Lineas y puntos base
linea_f, = graf_ani.plot(xi,u[:,0])
graf_ani.axhline(0, color='k')  # Eje en x=0

maximoy = np.max(u) 
txt_y = maximoy*(1-0.1)
txt_f = graf_ani.text(xi[0],txt_y,'tiempo:',
                      horizontalalignment='left')
# Configura gráfica
graf_ani.set_xlim([xi[0],xi[ultimox]])
graf_ani.set_ylim([ np.min(u),maximoy])
graf_ani.set_title(unmetodo)
graf_ani.set_xlabel('x')
graf_ani.set_ylabel('u')
#graf_ani.legend()
graf_ani.grid()

# Nueva Trama
def unatrama(i,xi,u):
    # actualiza cada linea
    linea_f.set_ydata(u[:,i])
    linea_f.set_xdata(xi)
    linea_f.set_label('tiempo linea: '+str(i))
    if (i<=9): # color de la línea
        lineacolor = 'C'+str(i)
    else:
        numcolor = i%10
        lineacolor = 'C'+str(numcolor)
    linea_f.set_color(lineacolor)
    txt_f.set_text('tiempo['+str(i)+'] = ' + str(i*dt))
    return linea_f, txt_f,

# Limpia Trama anterior
def limpiatrama(): 
    linea_f.set_ydata(np.ma.array(xi, mask=True))
    linea_f.set_label('')
    txt_f.set_text('')
    return linea_f, txt_f,

# contador de tramas
i = np.arange(0,len(tk),1)
ani = animation.FuncAnimation(fig_ani, unatrama,
                              i, fargs=(xi,u),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif',writer='pillow')
#ani.save(narchivo+'.mp4')
plt.draw()
plt.show()

1.5 Gráficas 3D animada, EDP Parabólica, método explícito

Para el caso de animación 3D, se cambia la sección de animación del algoritmo anterior por:

# GRAFICA 3D ANIMADA --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation
unmetodo = 'EDP Parabólica - Método explícito'
narchivo = 'EdpParabolica2' # nombre archivo GIF

# Parametros de trama/foto
retardo = 500   # milisegundos entre tramas

# GRAFICA animada en fig_ani
fig_ani3D = plt.figure()
graf_ani3 = fig_ani3D.add_subplot(projection='3d')

# Lineas y puntos base
##graf_ani3.plot_wireframe(Xi,Yi,U,
##                         color ='blue',label='EDP Parabólica')
punto_a, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',color ='green')
punto_b, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',
                          color ='salmon',label='U[i,j+1]')
linea_c, = graf_ani3.plot(xi[0],tk[0],U[1,0],'-',color ='blue')
# Configura gráfica
graf_ani3.set_xlim(min(xi),max(xi))
graf_ani3.set_ylim(min(tk),max(tk))
graf_ani3.set_title(unmetodo)
graf_ani3.set_xlabel('x')
graf_ani3.set_ylabel('t')
graf_ani3.set_zlabel('U')
graf_ani3.legend()
graf_ani3.grid()
graf_ani3.view_init(35, -45)

mallaEDP = None
# Nueva Trama
def unatrama(i,xi,tk,U,k):
    f = i//k
    c = i%k
    # actualiza cada punto
    punto_a.set_data([[xi[c],xi[c+1],xi[c+2]],
                      [tk[f],tk[f],tk[f]]])
    punto_a.set_3d_properties([U[f,c],U[f,c+1],U[f,c+2]])
    punto_b.set_data([xi[c+1]],[tk[f+1]])
    punto_b.set_3d_properties([U[f+1,c+1]])
    linea_c.set_data(xi[0:c+2],[tk[f+1]]*(c+2))
    linea_c.set_3d_properties([U[f+1,0:c+2]])
    
    global mallaEDP
    if mallaEDP:
        mallaEDP.remove()
    mallaEDP = graf_ani3.plot_wireframe(Xi[0:f+1,:],
                                        Yi[0:f+1,:],
                                        U[0:f+1,:],
                                        color ='blue',
                                        label='EDP Parabólica')
    return (punto_a,punto_b,linea_c)


# contador de tramas
k = len(xi)-2
k2 = k*(len(tk)-1)
i = np.arange(0,k2,1)
ani = animation.FuncAnimation(fig_ani3D, unatrama,
                              i ,
                              fargs=(xi,tk,U,k),
                              interval=retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='pillow')
#ani.save(narchivo+'.mp4')
#plt.draw()
plt.show()

1.6 Gráfica 3D, EDP Parabólicas método implícito

Puesto que la solución de la tabla de resultados en U es la misma que el los dos métodos explícito e implícito, solo se adjunta la parte de la animación.

Para el caso del la gráfica 3D del método implícito se cambia la sección de gráfica animada por:

# GRAFICA 3D ANIMADA --------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation
unmetodo = 'EDP Parabólica - Método explícito'
narchivo = 'EdpParabolica2' # nombre archivo GIF

# Parametros de trama/foto
retardo = 500   # milisegundos entre tramas

# GRAFICA animada en fig_ani
fig_ani3D = plt.figure()
graf_ani3 = fig_ani3D.add_subplot(projection='3d')

# Lineas y puntos base
##graf_ani3.plot_wireframe(Xi,Yi,U,
##                         color ='blue',label='EDP Parabólica')
punto_a, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',color ='green')
punto_b, = graf_ani3.plot(xi[0],tk[0],U[1,0],'o',
                          color ='salmon',label='U[i,j+1]')
linea_c, = graf_ani3.plot(xi[0],tk[0],U[1,0],'-',color ='blue')
# Configura gráfica
graf_ani3.set_xlim(min(xi),max(xi))
graf_ani3.set_ylim(min(tk),max(tk))
graf_ani3.set_title(unmetodo)
graf_ani3.set_xlabel('x')
graf_ani3.set_ylabel('t')
graf_ani3.set_zlabel('U')
graf_ani3.legend()
graf_ani3.grid()
graf_ani3.view_init(35, -45)

mallaEDP = None
# Nueva Trama
def unatrama(i,xi,tk,U,k):
    f = i//k
    c = i%k
    # actualiza cada punto
    punto_a.set_data([[xi[c],xi[c+1],xi[c+2]],
                      [tk[f],tk[f],tk[f]]])
    punto_a.set_3d_properties([U[f,c],U[f,c+1],U[f,c+2]])
    punto_b.set_data([xi[c+1]],[tk[f+1]])
    punto_b.set_3d_properties([U[f+1,c+1]])
    linea_c.set_data(xi[0:c+2],[tk[f+1]]*(c+2))
    linea_c.set_3d_properties([U[f+1,0:c+2]])
    
    global mallaEDP
    if mallaEDP:
        mallaEDP.remove()
    mallaEDP = graf_ani3.plot_wireframe(Xi[0:f+1,:],
                                        Yi[0:f+1,:],
                                        U[0:f+1,:],
                                        color ='blue',
                                        label='EDP Parabólica')
    return (punto_a,punto_b,linea_c)


# contador de tramas
k = len(xi)-2
k2 = k*(len(tk)-1)
i = np.arange(0,k2,1)
ani = animation.FuncAnimation(fig_ani3D, unatrama,
                              i ,
                              fargs=(xi,tk,U,k),
                              interval=retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='pillow')
#ani.save(narchivo+'.mp4')
#plt.draw()
plt.show()

2. EDP Elípticas método iterativo

Al algoritmo desarrollado en EDP Elípticas método iterativo, se crean una lista para almacenar los valores de la matriz U por cada iteración, se usan para graficar cada iteración.

# 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/2,dx)
yj = np.arange(y0,yn+dy/2,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
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] = (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
    U_3D.append(np.copy(u))
    U_error.append(np.around(errado_u,verdigitos))
    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])

2.1 Gráfica 3D, EDP Elípticas método iterativo

# 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()

2.2 Gráfica 3D animada, EDP Elípticas método iterativo

# GRAFICA CON ANIMACION 3D--------
# import matplotlib.pyplot as plt
import matplotlib.animation as animation
unmetodo = 'EDP Eliptica - Método iterativo'
narchivo = 'EdpElipticaItera' # nombre archivo GIF

# Parametros de trama/foto
retardo = 500   # milisegundos entre tramas

# GRAFICA animada en fig_ani
fig_ani3D = plt.figure()
graf_ani3 = fig_ani3D.add_subplot(projection='3d')

# Lineas y puntos base
punto_a, = graf_ani3.plot(xi[0],yj[0],U[0,0],'.',color ='blue')

# Configura gráfica
graf_ani3.set_xlim(min(xi),max(xi))
graf_ani3.set_ylim(min(yj),max(yj))
graf_ani3.set_title(unmetodo)
graf_ani3.set_xlabel('x')
graf_ani3.set_ylabel('t')
graf_ani3.set_zlabel('U')
# graf_ani3.legend()
graf_ani3.grid()
graf_ani3.view_init(35, -45)

etiqueta_i ='itera: '+str(0)
etiqueta_err = 'errado['+str(0)+']: '+str(U_error[0]) 
txt_i = graf_ani3.text2D(0.05, 0.95, etiqueta_i,
                         transform=graf_ani3.transAxes)
txt_errado = graf_ani3.text2D(0.05, 0.90, etiqueta_err,
                              transform=graf_ani3.transAxes)

mallaEDP = None
# Nueva Trama
def unatrama(i,xi,yj,U_3D,U_error):

    etiqueta_i ='itera: '+str(i)
    txt_i.set_text(etiqueta_i)
    etiqueta_err = 'errado['+str(i)+']: '+str(U_error[i])
    txt_errado.set_text(etiqueta_err)
    global mallaEDP
    if mallaEDP:
        graf_ani3.collections.remove(mallaEDP)
    mallaEDP = graf_ani3.plot_wireframe(Xi,Yi,
                                        np.transpose(U_3D[i]),
                                        color ='blue')
    return (txt_i,txt_errado,)


# contador de tramas
k = 11 #len(U_3D)
i = np.arange(0,k,1)
ani = animation.FuncAnimation(fig_ani3D, unatrama,
                              i ,
                              fargs=(xi,yj,U_3D,U_error),
                              interval=retardo,
                              blit=False)
# Graba Archivo GIFAnimado y video
ani.save(narchivo+'_animado.gif', writer='pillow')
#ani.save(narchivo+'.mp4')
#plt.draw()
plt.show()

2023-10-082025-09-07

7.2.4EDP Elípticas - analítico implícito con Sympy-Python

EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]
..

4. EDP Elípticas - Método Implícito con Sympy-Python

Desarrollo analítico del método implicito para una ecuación diferencial parcial elíptica usando Sympy. El algoritmo reutiliza el algoritmo para la EDP Elíptica de contínua a discreta, la creación de la matriz de valores u_xy, y la función
edp_sustituyeValorU() para buscar los valores conocidos de la u(x,y).

El algoritmo usa la ecuación discreta para en cada iteración i,j reemplazar los valores de U conocidos. Los valores de U se escriben en una matriz u_xy, para diferenciar si el valor existe se usa la matriz de estados u_mask.

Con las ecuaciones de cada iteración se llena la matriz A de coeficientes y el vector B de las constantes. Al resolver el sistema de ecuaciones se obtienen todos los valores de la matriz U, completando el ejercicio. Se desarrola la solución al sistema de ecuaciones usando Sympy como alternativa a usar numpy con np.linalg.solve(A.B) que se encuentra como comentario entre las instrucciones.

 discreta :
-4⋅U(i, j) + U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j) = 0

Método Implícito - EDP Elíptica
j: 1 ; i: 1 ; ecuacion: 0
U(0, 1) + U(1, 0) - 4⋅U(1, 1) + U(1, 2) + U(2, 1) = 0
-4⋅U(1, 1) + U(1, 2) + U(2, 1) = -110.0
j: 1 ; i: 2 ; ecuacion: 1
U(1, 1) + U(2, 0) - 4⋅U(2, 1) + U(2, 2) + U(3, 1) = 0
U(1, 1) - 4⋅U(2, 1) + U(2, 2) + U(3, 1) = -50.0
j: 1 ; i: 3 ; ecuacion: 2
U(2, 1) + U(3, 0) - 4⋅U(3, 1) + U(3, 2) + U(4, 1) = 0
U(2, 1) - 4⋅U(3, 1) + U(3, 2) + U(4, 1) = -50.0
j: 1 ; i: 4 ; ecuacion: 3
U(3, 1) + U(4, 0) - 4⋅U(4, 1) + U(4, 2) + U(5, 1) = 0
U(3, 1) - 4⋅U(4, 1) + U(4, 2) + U(5, 1) = -50.0
j: 1 ; i: 5 ; ecuacion: 4
U(4, 1) + U(5, 0) - 4⋅U(5, 1) + U(5, 2) + U(6, 1) = 0
U(4, 1) - 4⋅U(5, 1) + U(5, 2) + U(6, 1) = -50.0
j: 1 ; i: 6 ; ecuacion: 5
U(5, 1) + U(6, 0) - 4⋅U(6, 1) + U(6, 2) + U(7, 1) = 0
U(5, 1) - 4⋅U(6, 1) + U(6, 2) + U(7, 1) = -50.0
j: 1 ; i: 7 ; ecuacion: 6
U(6, 1) + U(7, 0) - 4⋅U(7, 1) + U(7, 2) + U(8, 1) = 0
U(6, 1) - 4⋅U(7, 1) + U(7, 2) = -75.0
j: 2 ; i: 1 ; ecuacion: 7
U(0, 2) + U(1, 1) - 4⋅U(1, 2) + U(1, 3) + U(2, 2) = 0
U(1, 1) - 4⋅U(1, 2) + U(1, 3) + U(2, 2) = -60.0
j: 2 ; i: 2 ; ecuacion: 8
U(1, 2) + U(2, 1) - 4⋅U(2, 2) + U(2, 3) + U(3, 2) = 0
j: 2 ; i: 3 ; ecuacion: 9
U(2, 2) + U(3, 1) - 4⋅U(3, 2) + U(3, 3) + U(4, 2) = 0
j: 2 ; i: 4 ; ecuacion: 10
U(3, 2) + U(4, 1) - 4⋅U(4, 2) + U(4, 3) + U(5, 2) = 0
j: 2 ; i: 5 ; ecuacion: 11
U(4, 2) + U(5, 1) - 4⋅U(5, 2) + U(5, 3) + U(6, 2) = 0
j: 2 ; i: 6 ; ecuacion: 12
U(5, 2) + U(6, 1) - 4⋅U(6, 2) + U(6, 3) + U(7, 2) = 0
j: 2 ; i: 7 ; ecuacion: 13
U(6, 2) + U(7, 1) - 4⋅U(7, 2) + U(7, 3) + U(8, 2) = 0
U(6, 2) + U(7, 1) - 4⋅U(7, 2) + U(7, 3) = -25.0
j: 3 ; i: 1 ; ecuacion: 14
U(0, 3) + U(1, 2) - 4⋅U(1, 3) + U(1, 4) + U(2, 3) = 0
U(1, 2) - 4⋅U(1, 3) + U(1, 4) + U(2, 3) = -60.0
j: 3 ; i: 2 ; ecuacion: 15
U(1, 3) + U(2, 2) - 4⋅U(2, 3) + U(2, 4) + U(3, 3) = 0
j: 3 ; i: 3 ; ecuacion: 16
U(2, 3) + U(3, 2) - 4⋅U(3, 3) + U(3, 4) + U(4, 3) = 0
j: 3 ; i: 4 ; ecuacion: 17
U(3, 3) + U(4, 2) - 4⋅U(4, 3) + U(4, 4) + U(5, 3) = 0
j: 3 ; i: 5 ; ecuacion: 18
U(4, 3) + U(5, 2) - 4⋅U(5, 3) + U(5, 4) + U(6, 3) = 0
j: 3 ; i: 6 ; ecuacion: 19
U(5, 3) + U(6, 2) - 4⋅U(6, 3) + U(6, 4) + U(7, 3) = 0
j: 3 ; i: 7 ; ecuacion: 20
U(6, 3) + U(7, 2) - 4⋅U(7, 3) + U(7, 4) + U(8, 3) = 0
U(6, 3) + U(7, 2) - 4⋅U(7, 3) + U(7, 4) = -25.0
j: 4 ; i: 1 ; ecuacion: 21
U(0, 4) + U(1, 3) - 4⋅U(1, 4) + U(1, 5) + U(2, 4) = 0
U(1, 3) - 4⋅U(1, 4) + U(1, 5) + U(2, 4) = -60.0
j: 4 ; i: 2 ; ecuacion: 22
U(1, 4) + U(2, 3) - 4⋅U(2, 4) + U(2, 5) + U(3, 4) = 0
j: 4 ; i: 3 ; ecuacion: 23
U(2, 4) + U(3, 3) - 4⋅U(3, 4) + U(3, 5) + U(4, 4) = 0
j: 4 ; i: 4 ; ecuacion: 24
U(3, 4) + U(4, 3) - 4⋅U(4, 4) + U(4, 5) + U(5, 4) = 0
j: 4 ; i: 5 ; ecuacion: 25
U(4, 4) + U(5, 3) - 4⋅U(5, 4) + U(5, 5) + U(6, 4) = 0
j: 4 ; i: 6 ; ecuacion: 26
U(5, 4) + U(6, 3) - 4⋅U(6, 4) + U(6, 5) + U(7, 4) = 0
j: 4 ; i: 7 ; ecuacion: 27
U(6, 4) + U(7, 3) - 4⋅U(7, 4) + U(7, 5) + U(8, 4) = 0
U(6, 4) + U(7, 3) - 4⋅U(7, 4) + U(7, 5) = -25.0
j: 5 ; i: 1 ; ecuacion: 28
U(0, 5) + U(1, 4) - 4⋅U(1, 5) + U(1, 6) + U(2, 5) = 0
U(1, 4) - 4⋅U(1, 5) + U(2, 5) = -130.0
j: 5 ; i: 2 ; ecuacion: 29
U(1, 5) + U(2, 4) - 4⋅U(2, 5) + U(2, 6) + U(3, 5) = 0
U(1, 5) + U(2, 4) - 4⋅U(2, 5) + U(3, 5) = -70.0
j: 5 ; i: 3 ; ecuacion: 30
U(2, 5) + U(3, 4) - 4⋅U(3, 5) + U(3, 6) + U(4, 5) = 0
U(2, 5) + U(3, 4) - 4⋅U(3, 5) + U(4, 5) = -70.0
j: 5 ; i: 4 ; ecuacion: 31
U(3, 5) + U(4, 4) - 4⋅U(4, 5) + U(4, 6) + U(5, 5) = 0
U(3, 5) + U(4, 4) - 4⋅U(4, 5) + U(5, 5) = -70.0
j: 5 ; i: 5 ; ecuacion: 32
U(4, 5) + U(5, 4) - 4⋅U(5, 5) + U(5, 6) + U(6, 5) = 0
U(4, 5) + U(5, 4) - 4⋅U(5, 5) + U(6, 5) = -70.0
j: 5 ; i: 6 ; ecuacion: 33
U(5, 5) + U(6, 4) - 4⋅U(6, 5) + U(6, 6) + U(7, 5) = 0
U(5, 5) + U(6, 4) - 4⋅U(6, 5) + U(7, 5) = -70.0
j: 5 ; i: 7 ; ecuacion: 34
U(6, 5) + U(7, 4) - 4⋅U(7, 5) + U(7, 6) + U(8, 5) = 0
U(6, 5) + U(7, 4) - 4⋅U(7, 5) = -95.0

 A : 
 [[-4.  1.  0. ...  0.  0.  0.]
 [ 1. -4.  1. ...  0.  0.  0.]
 [ 0.  1. -4. ...  0.  0.  0.]
 ...
 [ 0.  0.  0. ... -4.  1.  0.]
 [ 0.  0.  0. ...  1. -4.  1.]
 [ 0.  0.  0. ...  0.  1. -4.]]

 B : 
 [-110.  -50.  -50.  -50.  -50.  -50.  -75.  -60.    0.    0.    0.    0.
    0.  -25.  -60.    0.    0.    0.    0.    0.  -25.  -60.    0.    0.
    0.    0.    0.  -25. -130.  -70.  -70.  -70.  -70.  -70.  -95.]
Resultados para U(x,y)
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 ]
 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.81 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.]

Instrucciones en Python

Las instrucciones completas con Sympy-Python son:

# Ecuaciones Diferenciales Parciales Elipticas
# EDP Elípticas contínua a discreta con Sympy
import numpy as np
import sympy as sym

# u(x,y) funciones continuas y variables simbólicas usadas
x = sym.Symbol('x',real=True)
y = sym.Symbol('y',real=True)
u = sym.Function('u')(x,y) # funcion
f = sym.Function('f')(x,y) # funcion complemento
# U[i,j] funciones discretas y variables simbólicas usadas
i  = sym.Symbol('i',integer=True,positive=True)
j  = sym.Symbol('j',integer=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion de Poisson
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) + sym.diff(u,y,2)
RHS = fxy
edp = sym.Eq(LHS-RHS,0)

# centrada, centrada, atras
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,2): (U.subs(j,j-1)-2*U+U.subs(j,j+1))/(Dy**2),
               sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}

# Condiciones iniciales en los bordes
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 25 +0*y  # derecha
fxc = lambda x: 50 +0*x  # inferior 
fxd = lambda x: 70 +0*x  # superior 

# dimensiones de la placa
x0 = 0    # longitud en x
xn = 2
y0 = 0    # longitud en y
yn = 1.5
# muestreo en ejes, discreto, supone dx=dy
dx = 0.25  # Tamaño de paso
dy = dx # supone dx=dy
iteramax = 100 # revisa convergencia
tolera = 0.0001

verdigitos = 2      # para mostrar en tabla de resultados
casicero = 1e-15    # para redondeo de términos en ecuacion

# PROCEDIMIENTO
def edp_discreta(edp,dif_dividida,x,y,u):
    ''' EDP contínua a discreta, usa diferencias divididas
        proporcionadas en parámetros, indica las variables x,y
        con función u de (x,y)
    '''
    resultado={}
    # expresión todo a la izquierda LHS (Left Hand side)
    LHS = edp.lhs
    RHS = edp.rhs
    if not(edp.rhs==0):
        LHS = LHS-RHS
        RHS = 0
        edp = sym.Eq(LHS,RHS)
    # orden de derivada por x, y
    edp_x = edp.subs(x,0)
    edp_y = edp.subs(y,0)
    ordenDx = sym.ode_order(edp_x,u)
    ordenDy = sym.ode_order(edp_y,u)
    resultado['ordenDx'] = ordenDx  # guarda en resultados
    resultado['ordenDy'] = ordenDy
    # coeficiente derivada orden mayor a 1 (d2u/dx2)
    coeff_x = edp_coef_Dx(edp,x,ordenDx)
    if not(coeff_x==1):
        LHS = LHS/coeff_x
        RHS = RHS/coeff_x
    edp = sym.Eq(LHS,RHS)
    K_ = edp_coef_Dx(edp,y,ordenDy)
    if abs(K_)%1<casicero: # si es entero
        K_ = int(K_)
    resultado['edp=0']  = edp
    resultado['K_']  = K_
    
    discreta = edp.lhs  # EDP discreta
    for derivada in dif_dividida: # reemplaza diferencia dividida
        discreta = discreta.subs(derivada,dif_dividida[derivada])
    resultado['discreta=0'] = discreta
    return (resultado)

def edp_coef_Dx(edp,x,ordenx):
    ''' Extrae el coeficiente de la derivada Dx de ordenx
    edp es la ecuación como lhs=rhs
    '''
    coeff_x = 1.0 # valor inicial
    # separa cada término de suma
    term_suma = sym.Add.make_args(edp.lhs)
    for term_k in term_suma:
        if term_k.is_Mul: # mas de un factor
            factor_Mul = sym.Mul.make_args(term_k)
            coeff_temp = 1; coeff_usar=False
            # separa cada factor de término 
            for factor_k in factor_Mul:
                if not(factor_k.is_Derivative):
                    coeff_temp = coeff_temp*factor_k
                else: # factor con derivada de ordenx
                    partes = factor_k.args
                    if partes[1]==(x,ordenx):
                        coeff_usar = 1
            if coeff_usar==True:
                coeff_x = coeff_x*coeff_temp
    return(coeff_x)

def redondea_coef(ecuacion, precision=6,casicero = 1e-15):
    ''' redondea coeficientes de términos suma de una ecuacion
    ecuación como lhs=rhs
    '''
    tipo = type(ecuacion)
    tipo_eq = False
    if tipo == sym.core.relational.Equality:
        RHS = ecuacion.rhs
        ecuacion = ecuacion.lhs
        tipo = type(ecuacion)
        tipo_eq = True

    if tipo == sym.core.add.Add: # términos suma de ecuacion
        term_sum = sym.Add.make_args(ecuacion)
        ecuacion = sym.S.Zero # vacia
        for term_k in term_sum:
            # factor multiplicativo de termino suma
            term_mul = sym.Mul.make_args(term_k)
            producto = sym.S.One
            for factor in term_mul:
                if not(factor.has(sym.Symbol)): # es numerico
                    factor = np.around(float(factor),precision)
                    if (abs(factor)%1)<casicero: # si es entero
                        factor = int(factor)
                producto = producto*factor
            ecuacion = ecuacion + producto
    
    if tipo == sym.core.mul.Mul: # termino único, busca factores
        term_mul = sym.Mul.make_args(ecuacion)
        producto = sym.S.One
        for factor in term_mul:
            if not(factor.has(sym.Symbol)): # es numerico
                factor = np.around(float(factor),precision)
                if (abs(factor)%1)<casicero: # si es entero
                    factor = int(factor)
            producto = producto*factor
        ecuacion = producto
    
    if tipo == float: # # solo un numero
        if (abs(ecuacion)%1)<casicero: 
            ecuacion = int(ecuacion)
    if tipo_eq==True: # era igualdad, integra lhs=rhs
        ecuacion = sym.Eq(ecuacion,RHS)
    
    return(ecuacion)

def edp_simplificaLamba(resultado,dx,dy):
    '''simplifica ecuacion con valores de lambda, dx y dy
    entregando la edp discreta simplificada
    '''
    discreta = resultado['discreta=0']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    K_ = resultado['K_']
    lamb = (Dy**ordenDy)/(Dx**ordenDx)
    if ordenDy==1 and ordenDx==2:
        lamb = lamb/K_
    resultado['Lambda_L'] = lamb
    # valor de Lambda en ecuacion edp
    L_k = lamb.subs([(Dx,dx),(Dy,dy)])
    if abs(L_k)%1<casicero: # si es entero
        L_k = int(L_k)
    resultado['Lambda L_k'] = L_k
    # simplifica con lambda L
    discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
    resultado['(discreta=0)*Dy**ordeny'] = discreta_L
    discreta_L = edp_sustituye_L(resultado)
    discreta_L = discreta_L.subs(lamb,L)
    discreta_L = sym.collect(discreta_L,U)
    discreta_L = sym.Eq(discreta_L,0)
    resultado['discreta_L = 0'] = discreta_L
    # sustituye constantes en ecuación a iterar
    discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),(L,L_k)])
    discreta_L = redondea_coef(discreta_L)
    resultado['discreta'] = discreta_L
    return(resultado)

def edp_sustituye_L(resultado):
    ''' sustituye lambda con Dy**ordeny/Dx**x/K_
    por L, al simplificar Lambda
    '''
    discreta = resultado['(discreta=0)*Dy**ordeny']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    discreta_L = 0
    # separa cada término de suma
    term_suma = sym.Add.make_args(discreta)
    for term_k in term_suma:
        # busca partes de L y cambia por valor L
        cambiar = 0 # por orden de derivada
        if term_k.has(Dx) and term_k.has(Dy):
            partes = term_k.args
            ordeny=1
            ordenx=1
            for unaparte in partes:    
                if unaparte.has(Dy):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordeny = partey[1]
                if unaparte.has(Dx):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordenx = partey[1]
            if (ordeny<=ordenDy and ordenx<=-ordenDx):
                cambiar=1
        if cambiar:
            term_k = term_k*L/resultado['Lambda_L']
        discreta_L = discreta_L + term_k
        # simplifica unos con decimal a entero 
        discreta_L = discreta_L.subs(1.0,1)
    return(discreta_L)

def edp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
    '''Sustituye en edp discreta los valores conocidos de U[i,j]
    tomados desde u_xy, marcados con u_mask
    u_mask indica si el valor se ha calculado con edp.
    '''
    LHS = discreta.lhs # lado izquierdo de ecuacion
    RHS = discreta.rhs # lado derecho
    # sustituye U[i,j] con valores conocidos
    A_diagonal = [] # lista de i,j para matriz de coeficientes A
    # Separa términos suma
    term_suma = sym.Add.make_args(LHS)
    for term_k in term_suma:
        # busca U[i,j] y cambia por valor uxt[i,j]
        cambiar = 0 ; cambiar_valor = 0 ; cambiar_factor = 0
        # separa cada factor de término
        factor_Mul = sym.Mul.make_args(term_k)
        for factor_k in factor_Mul:
            # busca U[i,j] en matriz uxt[i,j]
            if factor_k.is_Function:
                [i_k,j_k] = factor_k.args
                if not(np.isnan(u_xy[i_k,j_k])):
                    cambiar = u_mask[i_k,j_k]
                    cambiar_factor = factor_k
                    cambiar_valor = u_xy[i_k,j_k]
                else:
                    A_diagonal.append([i_k,j_k,term_k/factor_k])
        # encontró valor U[i,j],term_k va a la derecha de ecuación
        if cambiar:
            LHS = LHS - term_k
            term_ki = term_k.subs(cambiar_factor,cambiar_valor)
            RHS = RHS - term_ki
    discreta = sym.Eq(LHS,RHS)
    B_diagonal = RHS
    resultado = [discreta,A_diagonal,B_diagonal]
    return (resultado)

# PROCEDIMIENTO
# transforma edp continua a discreta
resultado = edp_discreta(edp,dif_dividida,x,y,u)
resultado = edp_simplificaLamba(resultado,dx,dy)
discreta = resultado['discreta']

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Elíptica contínua a discreta')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    else:
        print('\n',entrada,':')
        sym.pprint(resultado[entrada])

# ITERAR para cada i,j dentro de U ------------
# x[i] , y[j]  valor en posición en cada eje
xi = np.arange(x0,xn+dx/2,dx)
yj = np.arange(y0,yn+dy/2,dy)
n = len(xi)
m = len(yj)

# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u
# llena u con valores en fronteras
u_xy[0,:]   = fya(yj)  # izquierda Ta
u_xy[n-1,:] = fyb(yj)  # derecha   Tb
u_xy[:,0]   = fxc(xi)  # inferior  Tc
u_xy[:,m-1] = fxd(xi)  # superior  Td
u0 = np.copy(u_xy)     # matriz u inicial

# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:]  = True # izquierda
u_mask[-1,:] = True # derecha
u_mask[:,0]  = True # inferior
u_mask[:,-1] = True # superior

# Método implícito para EDP Elíptica
# ITERAR para plantear las ecuaciones en [i,j]
resultado = {}
eq_itera = [] ; tamano = (n-2)*(m-2)
A = np.zeros(shape=(tamano,tamano),dtype=float)
B = np.zeros(tamano,dtype=float)
for j_k in range(1,m-1,1): # no usar valores en bordes
    for i_k in range(1,n-1,1): 
        eq_conteo = (j_k - 1)*(n-2)+(i_k-1)
        discreta_ij = discreta.subs({i:i_k,j:j_k,
                                     x:xi[i_k],y:yj[j_k]})
        resultado[eq_conteo]= {'j':j_k, 'i':i_k,
                                 'discreta_ij': discreta_ij}
        # usa valores de frontera segun u_mask con True
        discreta_k = edp_sustituyeValorU(discreta_ij,
                                        xi,yj,u_xy,u_mask)
        discreta_ij = discreta_k[0]
        A_diagonal  = discreta_k[1] # lista de (i,j,coeficiente) 
        B_diagonal  = discreta_k[2]
        resultado[eq_conteo]['discreta_k'] = discreta_k[0]
        # añade ecuacion a resolver
        eq_itera.append(discreta_ij)
        # Aplica coeficientes de ecuacion en A y B:
        # A_diagonal tiene lista de (i,j,coeficiente) 
        for uncoeff in A_diagonal:
            columna = (uncoeff[1]-1)*(n-2)+(uncoeff[0]-1)
            fila = (j_k - 1)*(n-2)+(i_k-1)
            A[fila,columna] = uncoeff[2]
        B[eq_conteo] = float(B_diagonal) # discreta_ij.rhs
resultado['A'] = np.copy(A)
resultado['B'] = np.copy(B)

# resuelve el sistema de ecuaciones en eq_itera en Sympy
X_k = sym.solve(eq_itera)[0]
# actualiza uxt[i,j] , u_mask segun X_k en Sympy
for nodo_Uij in X_k: 
    [i_k,j_k] = nodo_Uij.args
    u_xy[i_k,j_k] = X_k[nodo_Uij]
    u_mask[i_k,j_k] = True
# resuelve el sistema A.X=B en Numpy
#X = np.linalg.solve(A,B)
# tarea: llevar valores X a u_xy

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('\nMétodo Implícito - EDP Elíptica')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    elif (tipo==dict):
        print('j:',resultado[entrada]['j'],'; '
              'i:',resultado[entrada]['i'],'; '
              'ecuacion:',entrada)
        sym.pprint(resultado[entrada]['discreta_ij'])
        if resultado[entrada]['discreta_k'].rhs!=0:
            sym.pprint(resultado[entrada]['discreta_k'])
    else:
        print('\n',entrada,': \n',resultado[entrada])
print('Resultados para U(x,y)')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k in range(m-1,-1,-1):
    print(j_k, (u_xy[:,j_k]))

EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]

2023-10-072025-09-07

7.2.3 EDP Elípticas - analítico iterativo con Sympy-Python

EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]
..

1. Ejercicio

Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0

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

(ecuación de Laplace, Ecuación de Poisson con f(x,y)=0)

EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]

2. EDP Elípticas contínua a discreta

El resultado de las operaciones en la expresión usando el algoritmo es:

EDP Elíptica contínua a discreta
 ordenDx : 2
 ordenDy : 2

 edp=0 :
  2              2             
 ∂              ∂              
───(u(x, y)) + ───(u(x, y)) = 0
  2              2             
∂x             ∂y              
 K_ : 1

 discreta=0 :
-2⋅U(i, j) + U(i, j - 1) + U(i, j + 1)   -2⋅U(i, j) + U(i - 1, j) + U(i + 1, j)
────────────────────────────────────── + ──────────────────────────────────────
                   2                                        2                 
                 Dy                                       Dx                  
 Lambda_L :
  2
Dy 
───
  2
Dx 
 Lambda L_k : 1

 (discreta=0)*Dy**ordeny :
                                             2             2                 2
                                         2⋅Dy ⋅U(i, j)   Dy ⋅U(i - 1, j)   Dy ⋅ U(i + 1, j)
-2⋅U(i, j) + U(i, j - 1) + U(i, j + 1) - ───────────── + ─────────────── + ───────────────
                                                2                2                 2   
                                              Dx               Dx                Dx  

 discreta_L = 0 :
L⋅U(i - 1, j) + L⋅U(i + 1, j) + (-2⋅L - 2)⋅U(i, j) + U(i, j - 1) + U(i, j + 1) = 0

 discreta :
-4⋅U(i, j) + U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j) = 0

discreta_iterativa
          U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j)
U(i, j) = ─────────────────────────────────────────────────────
                                    4

El desarrollo analítico inicia convirtiendo la ecuación diferencial parcial elíptica de la forma contínua a su representación discreta usando las expresiones de derivadas en diferencias divididas. El cambio de la expresión se realiza usando Sympy.

La EDP se escribe en formato Sympy en el bloque de ingreso. La ecuación EDP se compone del lado izquierdo (LHS) y lado derecho (RHS), indicando u como una función de las variables x,y.

# INGRESO
fxy = 0*x+0*y # f(x,y) = 0 , ecuacion de Poisson
# ecuacion: LHS=RHS
LHS = sym.diff(u,x,2) + sym.diff(u,y,2)
RHS = fxy
EDP = sym.Eq(LHS-RHS,0)

las diferencias divididas se ingresan como un diccionario:

# centrada, centrada, atras
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,2): (U.subs(j,j-1)-2*U+U.subs(j,j+1))/(Dy**2),
               sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}

las condiciones iniciales en los bordes, pueden ser funciones matemáticas, por lo que se usa el formato lambda. En el ejercicio básico presentado en la parte teórica, los bordes tienen temperaturas constantes, por lo que en ésta sección se escriben las condiciones como la constante mas cero por la variable independiente, facilitando la evaluación de un vector xi por ejemplo.

# Condiciones iniciales en los bordes
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 25 +0*y  # derecha
fxc = lambda x: 50 +0*x  # inferior, función inicial 
fxd = lambda x: 70 +0*x  # superior, función inicial

los demás parámetros del ejercicio se ingresan más adelante de forma semejante al ejercicio con Numpy.

Todas las expresiones se escriben al lado izquierdo (LHS) de la igualdad, la parte del lado derecho(RHS) se la prefiere mantener en cero. La expresión se organiza manteniendo el coeficiente en 1 para el termino de Dx de mayor orden. Se sustituyen las derivadas por diferencias divididas, para luego simplificar la expresión usando λ como referencia.

Instrucciones en Python

Los resultados al aplicar una operación a la expresión se guardan en un diccionario, para revisar cada paso intermedio al final

# Ecuaciones Diferenciales Parciales Elipticas
# EDP Elípticas contínua a discreta con Sympy
import numpy as np
import sympy as sym

# u(x,y) funciones continuas y variables simbólicas usadas
x = sym.Symbol('x',real=True)
y = sym.Symbol('y',real=True)
u = sym.Function('u')(x,y) # funcion
f = sym.Function('f')(x,y) # funcion complemento
# U[i,j] funciones discretas y variables simbólicas usadas
i  = sym.Symbol('i',integer=True,positive=True) # indices
j  = sym.Symbol('j',integer=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion de Poisson
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) + sym.diff(u,y,2)
RHS = fxy
edp = sym.Eq(LHS-RHS,0)

# centrada, centrada, atras
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,2): (U.subs(j,j-1)-2*U+U.subs(j,j+1))/(Dy**2),
               sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}

# Condiciones iniciales en los bordes
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 25 +0*y  # derecha
fxc = lambda x: 50 +0*x  # inferior 
fxd = lambda x: 70 +0*x  # superior 

# dimensiones de la placa
x0 = 0    # longitud en x
xn = 2
y0 = 0    # longitud en y
yn = 1.5
# muestreo en ejes, discreto, supone dx=dy
dx = 0.25  # Tamaño de paso
dy = dx # supone dx=dy
iteramax = 100 # revisa convergencia
tolera = 0.0001

verdigitos = 2      # para mostrar en tabla de resultados
casicero = 1e-15    # para redondeo de términos en ecuacion

# PROCEDIMIENTO
def edp_discreta(edp,dif_dividida,x,y,u):
    ''' EDP contínua a discreta, usa diferencias divididas
        proporcionadas en parámetros, indica las variables x,y
        con función u de (x,y)
    '''
    resultado={}
    # expresión todo a la izquierda LHS (Left Hand side)
    LHS = edp.lhs
    RHS = edp.rhs
    if not(edp.rhs==0):
        LHS = LHS-RHS
        RHS = 0
        edp = sym.Eq(LHS,RHS)
    # orden de derivada por x, y
    edp_x = edp.subs(x,0)
    edp_y = edp.subs(y,0)
    ordenDx = sym.ode_order(edp_x,u)
    ordenDy = sym.ode_order(edp_y,u)
    resultado['ordenDx'] = ordenDx  # guarda en resultados
    resultado['ordenDy'] = ordenDy
    # coeficiente derivada orden mayor a 1 (d2u/dx2)
    coeff_x = edp_coef_Dx(edp,x,ordenDx)
    if not(coeff_x==1):
        LHS = LHS/coeff_x
        RHS = RHS/coeff_x
    edp = sym.Eq(LHS,RHS)
    K_ = edp_coef_Dx(edp,y,ordenDy)
    if abs(K_)%1<casicero: # si es entero
        K_ = int(K_)
    resultado['edp=0']  = edp
    resultado['K_']  = K_
    
    discreta = edp.lhs  # EDP discreta
    for derivada in dif_dividida: # reemplaza diferencia dividida
        discreta = discreta.subs(derivada,dif_dividida[derivada])
    resultado['discreta=0'] = discreta
    return (resultado)

def edp_coef_Dx(edp,x,ordenx):
    ''' Extrae el coeficiente de la derivada Dx de ordenx
    edp es la ecuación como lhs=rhs
    '''
    coeff_x = 1.0 # valor inicial
    # separa cada término de suma
    term_suma = sym.Add.make_args(edp.lhs)
    for term_k in term_suma:
        if term_k.is_Mul: # mas de un factor
            factor_Mul = sym.Mul.make_args(term_k)
            coeff_temp = 1; coeff_usar=False
            # separa cada factor de término 
            for factor_k in factor_Mul:
                if not(factor_k.is_Derivative):
                    coeff_temp = coeff_temp*factor_k
                else: # factor con derivada de ordenx
                    partes = factor_k.args
                    if partes[1]==(x,ordenx):
                        coeff_usar = 1
            if coeff_usar==True:
                coeff_x = coeff_x*coeff_temp
    return(coeff_x)

def redondea_coef(ecuacion, precision=6,casicero = 1e-15):
    ''' redondea coeficientes de términos suma de una ecuacion
    ecuación como lhs=rhs
    '''
    tipo = type(ecuacion)
    tipo_eq = False
    if tipo == sym.core.relational.Equality:
        RHS = ecuacion.rhs
        ecuacion = ecuacion.lhs
        tipo = type(ecuacion)
        tipo_eq = True

    if tipo == sym.core.add.Add: # términos suma de ecuacion
        term_sum = sym.Add.make_args(ecuacion)
        ecuacion = sym.S.Zero # vacia
        for term_k in term_sum:
            # factor multiplicativo de termino suma
            term_mul = sym.Mul.make_args(term_k)
            producto = sym.S.One
            for factor in term_mul:
                if not(factor.has(sym.Symbol)): # es numerico
                    factor = np.around(float(factor),precision)
                    if (abs(factor)%1)<casicero: # si es entero
                        factor = int(factor)
                producto = producto*factor
            ecuacion = ecuacion + producto
    
    if tipo == sym.core.mul.Mul: # termino único, busca factores
        term_mul = sym.Mul.make_args(ecuacion)
        producto = sym.S.One
        for factor in term_mul:
            if not(factor.has(sym.Symbol)): # es numerico
                factor = np.around(float(factor),precision)
                if (abs(factor)%1)<casicero: # si es entero
                    factor = int(factor)
            producto = producto*factor
        ecuacion = producto
    
    if tipo == float: # # solo un numero
        if (abs(ecuacion)%1)<casicero: 
            ecuacion = int(ecuacion)
    if tipo_eq==True: # era igualdad, integra lhs=rhs
        ecuacion = sym.Eq(ecuacion,RHS)
    
    return(ecuacion)

def edp_simplificaLamba(resultado,dx,dy):
    '''simplifica ecuacion con valores de lambda, dx y dy
    entregando la edp discreta simplificada
    '''
    discreta = resultado['discreta=0']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    K_ = resultado['K_']
    lamb = (Dy**ordenDy)/(Dx**ordenDx)
    if ordenDy==1 and ordenDx==2:
        lamb = lamb/K_
    resultado['Lambda_L'] = lamb
    # valor de Lambda en ecuacion edp
    L_k = lamb.subs([(Dx,dx),(Dy,dy)])
    if abs(L_k)%1<casicero: # si es entero
        L_k = int(L_k)
    resultado['Lambda L_k'] = L_k
    # simplifica con lambda L
    discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
    resultado['(discreta=0)*Dy**ordeny'] = discreta_L
    discreta_L = edp_sustituye_L(resultado)
    discreta_L = discreta_L.subs(lamb,L)
    discreta_L = sym.collect(discreta_L,U)
    discreta_L = sym.Eq(discreta_L,0)
    resultado['discreta_L = 0'] = discreta_L
    # sustituye constantes en ecuación a iterar
    discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),(L,L_k)])
    discreta_L = redondea_coef(discreta_L)
    resultado['discreta'] = discreta_L
    return(resultado)

def edp_sustituye_L(resultado):
    ''' sustituye lambda con Dy**ordeny/Dx**x/K_
    por L, al simplificar Lambda
    '''
    discreta = resultado['(discreta=0)*Dy**ordeny']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    discreta_L = 0
    # separa cada término de suma
    term_suma = sym.Add.make_args(discreta)
    for term_k in term_suma:
        # busca partes de L y cambia por valor L
        cambiar = False # por orden de derivada
        if term_k.has(Dx) and term_k.has(Dy):
            partes = term_k.args
            ordeny=1
            ordenx=1
            for unaparte in partes:    
                if unaparte.has(Dy):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordeny = partey[1]
                if unaparte.has(Dx):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordenx = partey[1]
            if (ordeny<=ordenDy and ordenx<=-ordenDx):
                cambiar=True
        if cambiar==True:
            term_k = term_k*L/resultado['Lambda_L']
        discreta_L = discreta_L + term_k
        # simplifica unos con decimal a entero 
        discreta_L = redondea_coef(discreta_L)
    return(discreta_L)

# PROCEDIMIENTO -------------------------------
# transforma edp continua a discreta
resultado = edp_discreta(edp,dif_dividida,x,y,u)
resultado = edp_simplificaLamba(resultado,dx,dy)
discreta = resultado['discreta']

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Elíptica contínua a discreta')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    else:
        print('\n',entrada,':')
        sym.pprint(resultado[entrada])

EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]
..

3. EDP Elípticas - Método iterativo con Sympy-Python

El desarrollo del método iterativo para una ecuación diferencial parcial elíptica usando Sympy, sigue a continuación del anterior.

Se añaden las instrucciones para iterar la expresión discreta para cada i,j dentro de la matriz U creada para el ejercicio. Se aplican las condiciones iniciales a los bordes, para luego proceder a iterar.

Como en las actividades del curso se requiere realizar al menos 3 iteraciones para las expresiones del algoritmo con papel y lápiz, solo se presentarán los resultados para la primera fila en j=0

Los resultados del algoritmo luego del resultado del algoritmo anterior se presentan como:

discreta_iterativa
          U(i, j - 1) + U(i, j + 1) + U(i - 1, j) + U(i + 1, j)
U(i, j) = ─────────────────────────────────────────────────────
                                    4                          

 iterar i,j:
j:  1  ;  i:  1
          U(0, 1)   U(1, 0)   U(1, 2)   U(2, 1)
U(1, 1) = ─────── + ─────── + ─────── + ───────
             4         4         4         4   
U(1, 1) = 53.125
j:  1  ;  i:  2
          U(1, 1)   U(2, 0)   U(2, 2)   U(3, 1)
U(2, 1) = ─────── + ─────── + ─────── + ───────
             4         4         4         4   
U(2, 1) = 51.40625
j:  1  ;  i:  3
          U(2, 1)   U(3, 0)   U(3, 2)   U(4, 1)
U(3, 1) = ─────── + ─────── + ─────── + ───────
             4         4         4         4   
U(3, 1) = 50.9765625
j:  1  ;  i:  4
          U(3, 1)   U(4, 0)   U(4, 2)   U(5, 1)
U(4, 1) = ─────── + ─────── + ─────── + ───────
             4         4         4         4   
U(4, 1) = 50.869140625
j:  1  ;  i:  5
          U(4, 1)   U(5, 0)   U(5, 2)   U(6, 1)
U(5, 1) = ─────── + ─────── + ─────── + ───────
             4         4         4         4   
U(5, 1) = 50.84228515625
j:  1  ;  i:  6
          U(5, 1)   U(6, 0)   U(6, 2)   U(7, 1)
U(6, 1) = ─────── + ─────── + ─────── + ───────
             4         4         4         4   
U(6, 1) = 50.8355712890625
j:  1  ;  i:  7
          U(6, 1)   U(7, 0)   U(7, 2)   U(8, 1)
U(7, 1) = ─────── + ─────── + ─────── + ───────
             4         4         4         4   
U(7, 1) = 44.271392822265625

Método iterativo EDP Elíptica
continuar el desarrollo con: 
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 ]
 j, U[i,j]
2 [60.   51.25 51.25 51.25 51.25 51.25 51.25 51.25 25.  ]
1 [60.   53.12 51.41 50.98 50.87 50.84 50.84 44.27 25.  ]
0 [50. 50. 50. 50. 50. 50. 50. 50. 50.]

Se pueden sustituir los valores de las expresiones, evaluando cada u[i,j] y completando la matriz y para la gráfica. Esta parte queda como tarea.

Instrucciones en Python

Las instrucciones adicionales al algoritmo de EDP elíptica contínua a discreta empiezan con la creación de la matriz u_xy para los valores, los vectores xi,yj y una matriz u_mask que indica si se dispone de un valor calculado o conocido en ese punto (x,y) . Para el método iterativo, la matriz se rellena con el promedio de los valores máximos de las funciones dadas para los bordes de la placa.

Tarea: El algoritmo desarrolla el cálculo para j_k=1, solo la primera fila. Se deja como tarea desarrollar para toda la matriz.

# ITERAR para cada i,j dentro de U ------------
# x[i] , y[j]  valor en posición en cada eje
xi = np.arange(x0,xn+dx/2,dx)
yj = np.arange(y0,yn+dy/2,dy)
n = len(xi)
m = len(yj)

# Matriz U 
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u
# llena u con valores en fronteras
u_xy[0,:]   = fya(yj)  # izquierda Ta
u_xy[n-1,:] = fyb(yj)  # derecha   Tb
u_xy[:,0]   = fxc(xi)  # inferior  Tc
u_xy[:,m-1] = fxd(xi)  # superior  Td
u0 = np.copy(u_xy)     # matriz u inicial

# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:]  = True # izquierda
u_mask[-1,:] = True # derecha
u_mask[:,0]  = True # inferior
u_mask[:,-1] = True # superior

# valor inicial de iteración dentro de u
# promedio = (Ta+Tb+Tc+Td)/4
promedio = (np.max(fya(yj))+np.max(fyb(yj))+\
            np.max(fxc(xi))+np.max(fxd(xi)))/4
u_xy[1:n-1,1:m-1]   = promedio
u_mask[1:n-1,1:m-1] = True
u0 = np.copy(u_xy) # copia para revisión

def edp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
    '''Sustituye en edp discreta los valores conocidos de U[i,j]
    tomados desde u_xy, marcados con u_mask
    u_mask indica si el valor se ha calculado con edp.
    '''
    LHS = discreta.lhs # lado izquierdo de ecuacion
    RHS = discreta.rhs # lado derecho
    # sustituye U[i,j] con valores conocidos
    A_diagonal = [] # lista de i,j para matriz de coeficientes A
    # Separa términos suma
    term_suma = sym.Add.make_args(LHS)
    for term_k in term_suma:
        # busca U[i,j] y cambia por valor uxt[i,j]
        cambiar = False ; cambiar_valor = 0 ; cambiar_factor = 0
        # separa cada factor de término
        factor_Mul = sym.Mul.make_args(term_k)
        for factor_k in factor_Mul:
            # busca U[i,j] en matriz uxt[i,j]
            if factor_k.is_Function:
                [i_k,j_k] = factor_k.args
                if not(np.isnan(u_xy[i_k,j_k])):
                    cambiar = u_mask[i_k,j_k]
                    cambiar_factor = factor_k
                    cambiar_valor = u_xy[i_k,j_k]
                else:
                    A_diagonal.append([i_k,j_k,term_k/factor_k])
        # encontró valor U[i,j],term_k va a la derecha de ecuación
        if cambiar==True:
            LHS = LHS - term_k
            term_ki = term_k.subs(cambiar_factor,cambiar_valor)
            RHS = RHS - term_ki
    discreta = sym.Eq(LHS,RHS)
    B_diagonal = RHS
    resultado = [discreta,A_diagonal,B_diagonal]
    return (resultado)

# Método iterativo para EDP Elíptica
# separa término U[j,j] del centro,con valor no conocido
buscar = U.subs(j,j)
discreta = sym.solve(discreta,buscar)
discreta = sym.factor_terms(discreta[0])
discreta = sym.Eq(buscar,discreta)
resultado['discreta_iterativa'] = discreta

print('\ndiscreta_iterativa')
sym.pprint(discreta)

# Desarrollo de iteraciones en cada nodo i,j, para la fila j_k
# genera el valor en cada punto
# Nota: Solo la primera fila, Tarea desarrollar para la matriz
print('\n iterar i,j:')
j_k = 1
for i_k in range(1,n-1,1):
    print('j: ',j_k,' ; ','i: ',i_k)
    discreta_ij = discreta.subs({i:i_k,j:j_k,
                                 x:xi[i_k],y:yj[j_k]})
    sym.pprint(discreta_ij)
    RHS = discreta_ij.rhs
    discreta_ij = sym.Eq(RHS,0)
    # usa valores de frontera segun u_mask con True
    discreta_k = edp_sustituyeValorU(discreta_ij,
                                    xi,yj,u_xy,u_mask)
    u_xy[i_k,j_k] = sym.Float(-discreta_k[2])
    print(buscar.subs({i:i_k,j:j_k}),
          '=',u_xy[i_k,j_k])
    
# SALIDA
np.set_printoptions(precision=verdigitos)
print('\nMétodo iterativo EDP Elíptica')
print('continuar el desarrollo con: ')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k in range(2,-1,-1):
    print(j_k, (u_xy[:,j_k]))

EDP Elípticas [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]
..

Grafica en 3D para EDP Elípticas

A partir de la solución en matriz U_xy, Xi, Yi

# GRAFICA en 3D ------
import matplotlib.pyplot as plt
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi, yj)
U_xy = np.transpose(u_xy) # 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_xy,
                       color ='blue')
graf_3D.plot(Xi[1,0],Yi[1,0],U_xy[1,0],'o',color ='orange')
graf_3D.plot(Xi[1,1],Yi[1,1],U_xy[1,1],'o',color ='red',
             label='U[i,j]')
graf_3D.plot(Xi[1,2],Yi[1,2],U_xy[1,2],'o',color ='salmon')
graf_3D.plot(Xi[0,1],Yi[0,1],U_xy[0,1],'o',color ='green')
graf_3D.plot(Xi[2,1],Yi[2,1],U_xy[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 [ contínua a discreta ] || Sympy [ iterativo ] [ implícito ] [ gráfica_3D ]

2023-10-062025-09-07

7.1.4 EDP Parabólica - analítico implícito con Sympy-Python

EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]
..

4. EDP Parabólica - Método implícito con Sympy

Para el ejercicio presentado en el numeral 1, se desarrolla el método implícito.

Las diferencias finitas centradas y hacia atrás se definen en el diccionario dif_dividida.

Para desarrollar la expresión discreta, se incorpora la función edp_sustituye_L(), que busca el en cada término de la suma los componentes de λ y los sustituye por la variable L.

Para cada valor de j, se crean las ecuaciones desde la forma discreta de la EDP, moviendo a la derecha los valores conocidos para generar un sistema de ecuaciones A.X=B. Se resuelve el sistema de ecuaciones y se actualiza la matriz u[i,j] de resultados.

Por la extensión de los pasos a mostrar, para las iteraciones en cada fila, solo se presentan para la primera fila. En caso de requerir observar más filas, puede cambiar el parámetro en el bloque de salida, en el condicional:

    elif (tipo==dict):
        if entrada<2

Los resultados para el algoritmo son:

EDP Parabólica contínua a discreta
 ordenDx : 2
 ordenDy : 1

 edp=0 :
  2                             
 ∂               ∂              
───(u(x, y)) - 4⋅──(u(x, y)) = 0
  2              ∂y             
∂x                              
 K_ : 4

 Lambda_L :
  Dy 
─────
    2
4⋅Dx 
 Lambda L_k : 0.250000000000000

 (discreta=0)*Dy**ordeny :
                             2⋅Dy⋅U(i, j)   Dy⋅U(i - 1, j)   Dy⋅U(i + 1, j)
-4⋅U(i, j) + 4⋅U(i, j - 1) - ──────────── + ────────────── + ──────────────
                                   2               2                2      
                                 Dx              Dx               Dx       

 discreta_L=0 :
4⋅L⋅U(i - 1, j) + 4⋅L⋅U(i + 1, j) + (-8⋅L - 4)⋅U(i, j) + 4⋅U(i, j - 1) = 0

 discreta :
-6⋅U(i, j) + 4⋅U(i, j - 1) + U(i - 1, j) + U(i + 1, j) = 0

Método implícito EDP Parabólica
j: 1 ; i: 1
U(0, 1) + 4⋅U(1, 0) - 6⋅U(1, 1) + U(2, 1) = 0
-6⋅U(1, 1) + U(2, 1) = -160.0
j: 1 ; i: 2
U(1, 1) + 4⋅U(2, 0) - 6⋅U(2, 1) + U(3, 1) = 0
U(1, 1) - 6⋅U(2, 1) + U(3, 1) = -100.0
j: 1 ; i: 3
U(2, 1) + 4⋅U(3, 0) - 6⋅U(3, 1) + U(4, 1) = 0
U(2, 1) - 6⋅U(3, 1) + U(4, 1) = -100.0
j: 1 ; i: 4
U(3, 1) + 4⋅U(4, 0) - 6⋅U(4, 1) + U(5, 1) = 0
U(3, 1) - 6⋅U(4, 1) + U(5, 1) = -100.0
j: 1 ; i: 5
U(4, 1) + 4⋅U(5, 0) - 6⋅U(5, 1) + U(6, 1) = 0
U(4, 1) - 6⋅U(5, 1) + U(6, 1) = -100.0
j: 1 ; i: 6
U(5, 1) + 4⋅U(6, 0) - 6⋅U(6, 1) + U(7, 1) = 0
U(5, 1) - 6⋅U(6, 1) + U(7, 1) = -100.0
j: 1 ; i: 7
U(6, 1) + 4⋅U(7, 0) - 6⋅U(7, 1) + U(8, 1) = 0
U(6, 1) - 6⋅U(7, 1) + U(8, 1) = -100.0
j: 1 ; i: 8
U(7, 1) + 4⋅U(8, 0) - 6⋅U(8, 1) + U(9, 1) = 0
U(7, 1) - 6⋅U(8, 1) + U(9, 1) = -100.0
j: 1 ; i: 9
U(8, 1) + 4⋅U(9, 0) - 6⋅U(9, 1) + U(10, 1) = 0
U(8, 1) - 6⋅U(9, 1) = -140.0
A :
[[-6.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 1. -6.  1.  0.  0.  0.  0.  0.  0.]
 [ 0.  1. -6.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  1. -6.  1.  0.  0.  0.  0.]
 [ 0.  0.  0.  1. -6.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  1. -6.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  1. -6.  1.  0.]
 [ 0.  0.  0.  0.  0.  0.  1. -6.  1.]
 [ 0.  0.  0.  0.  0.  0.  0.  1. -6.]]
B :
[-160. -100. -100. -100. -100. -100. -100. -100. -140.]
X :
[31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57]
Resultados para U(x,y)
xi: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
yj: [0.   0.01 0.02 0.03 0.04]
 j, U[i,j]
4 [60.   40.67 30.61 26.75 25.51 25.18 25.24 25.76 27.41 31.71 40.  ]
3 [60.   38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72 40.  ]
2 [60.   35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39 40.  ]
1 [60.   31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]

Instrucciones en Python

Para el método implícito se añade la sección donde se reemplazan los valores en los bordes y en la barra en la matriz U(x,y) en el algoritmo identificada como u_xy . Para diferenciar si el valor existe se usa la matriz de estados u_mask con estados True/False.

También se incorpora la función edp_sustituyeValorU() que para cada término suma de la ecuación busca en la matriz u_xy si el valor U(i,j) existe y lo cambia. Luego pasa el valor al lado derecho de la ecuación par formar el sistema de ecuaciones.

Con los valores existentes antes de la iteración, se obtienen las ecuaciones por cada punto i,j de una fila para el sistema de ecuaciones en la forma A.x=B. El sistema se resuelve con las instrucciones de Numpy. El resultado se actualiza en la matriz u_xy y la matriz u_mask para indicar que el valor ya está disponible.

Compare los resultados con el método numérico y mida los tiempos de ejecución de cada algoritmo para escribir sus observaciones.

# Ecuaciones Diferenciales Parciales Parabólicas
# EDP Parabólicas contínua a discreta con Sympy
import numpy as np
import sympy as sym

# u(x,y) funciones continuas y variables simbólicas usadas
t = sym.Symbol('t',real=True) # independiente
x = sym.Symbol('x',real=True)
y = sym.Symbol('y',real=True) 
u = sym.Function('u')(x,y) # funcion
f = sym.Function('f')(x,y) # funcion complemento
K = sym.Symbol('K',real=True)
# U[i,j] funciones discretas y variables simbólicas usadas
i  = sym.Symbol('i',integer=True,positive=True) # indices
j  = sym.Symbol('j',integer=True,positive=True)
Dt = sym.Symbol('Dt',real=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion complementaria
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) 
RHS = 4*sym.diff(u,y,1) + fxy
edp = sym.Eq(LHS-RHS,0)

# centrada, atras
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,1): (U - U.subs(j,j-1))/Dy}

buscar = U.subs(j,j+1) # U[i,j+1] método explícito

# Valores de frontera
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 40 +0*y  # derecha
fxc = lambda x: 25 +0*x  # inferior, función inicial

# [a,b] dimensiones de la barra
a = 0  # longitud en x
b = 1
y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones

# muestreo en ejes, discreto
x_tramos = 10
dx  = (b-a)/x_tramos
dy  = dx/10
n_iteraciones = 4 # iteraciones en tiempo

verdigitos = 2      # para mostrar en tabla de resultados
casicero = 1e-15    # para redondeo de términos en ecuacion

# PROCEDIMIENTO

def edp_coef_Dx(edp,x,ordenx):
    ''' Extrae el coeficiente de la derivada Dx de ordenx,
    edp es la ecuación como lhs=rhs
    '''
    coeff_x = 1.0 # valor inicial
    # separa cada término de suma
    term_suma = sym.Add.make_args(edp.lhs)
    for term_k in term_suma:
        if term_k.is_Mul: # mas de un factor
            factor_Mul = sym.Mul.make_args(term_k)
            coeff_temp = 1; coeff_usar = False
            # separa cada factor de término 
            for factor_k in factor_Mul:
                if not(factor_k.is_Derivative):
                    coeff_temp = coeff_temp*factor_k
                else: # factor con derivada de ordenx
                    partes = factor_k.args
                    if partes[1]==(x,ordenx):
                        coeff_usar = True
            if coeff_usar==True:
                coeff_x = coeff_x*coeff_temp
    return(coeff_x)

def redondea_coef(ecuacion, precision=6,casicero = 1e-15):
    ''' redondea coeficientes de términos suma de una ecuacion
    ecuación como lhs=rhs
    '''
    tipo = type(ecuacion)
    tipo_eq = False
    if tipo == sym.core.relational.Equality:
        RHS = ecuacion.rhs  # separa lado derecho
        ecuacion = ecuacion.lhs # analiza lado izquierdo
        tipo = type(ecuacion)
        tipo_eq = True

    if tipo == sym.core.add.Add: # términos suma de ecuacion
        term_sum = sym.Add.make_args(ecuacion)
        ecuacion = sym.S.Zero  # vacia
        for term_k in term_sum:
            # factor multiplicativo de termino suma
            term_mul = sym.Mul.make_args(term_k)
            producto = sym.S.One
            for factor in term_mul:
                if not(factor.has(sym.Symbol)): # es numerico
                    factor = np.around(float(factor),precision)
                    if (abs(factor)%1)<casicero: # si es entero
                        factor = int(factor)
                producto = producto*factor
            ecuacion = ecuacion + producto
            
    if tipo == sym.core.mul.Mul: # termino único, busca factores
        term_mul = sym.Mul.make_args(ecuacion)
        producto = sym.S.One
        for factor in term_mul:
            if not(factor.has(sym.Symbol)): # es numerico
                factor = np.around(float(factor),precision)
                if (abs(factor)%1)<casicero: # si es entero
                    factor = int(factor)
            producto = producto*factor
        ecuacion = producto
        
    if tipo == float: # solo un numero
        if (abs(ecuacion)%1)<casicero: 
            ecuacion = int(ecuacion)
            
    if tipo_eq==True: # era igualdad, integra lhs=rhs
        ecuacion = sym.Eq(ecuacion,RHS)

    return(ecuacion)

def edp_sustituye_L(resultado):
    ''' sustituye lambda con Dy**ordeny/Dx**x/K_
    por L, al simplificar Lambda
    '''
    discreta = resultado['(discreta=0)*Dy**ordeny']
    ordenDy = resultado['ordenDy']
    ordenDx = resultado['ordenDx']
    discreta_L = 0
    # separa cada término de suma
    term_suma = sym.Add.make_args(discreta)
    for term_k in term_suma:
        # busca partes de L y cambia por valor L
        cambiar = False # por orden de derivada
        if term_k.has(Dx) and term_k.has(Dy):
            partes = term_k.args
            ordeny = 1
            ordenx = 1
            for unaparte in partes:    
                if unaparte.has(Dy):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordeny = partey[1]
                if unaparte.has(Dx):
                    if unaparte.is_Pow:
                        partey = unaparte.args
                        ordenx = partey[1]
            if (ordeny<=ordenDy and ordenx<=-ordenDx):
                cambiar = True
        if cambiar==True:
            term_k = term_k*L/resultado['Lambda_L']
        discreta_L = discreta_L + term_k
        # simplifica unos con decimal a entero 
        discreta_L = redondea_coef(discreta_L)
    return(discreta_L)


resultado = {} # resultados en diccionario
# orden de derivada por x, y
edp_x = edp.subs(x,0)
edp_y = edp.subs(y,0)
ordenDx = sym.ode_order(edp_x,u)
ordenDy = sym.ode_order(edp_y,u)
resultado['ordenDx'] = ordenDx  # guarda en resultados
resultado['ordenDy'] = ordenDy
# coeficiente derivada orden mayor a 1 (d2u/dx2)
coeff_x = edp_coef_Dx(edp,x,ordenDx)
LHS = edp.lhs  # lado izquierdo de edp
RHS = edp.rhs  # lado derecho de edp
if not(coeff_x==1):
    LHS = LHS/coeff_x
    RHS = RHS/coeff_x
# toda la expresión a la izquierda
edp = sym.Eq(LHS-RHS,0)
K_ = -edp_coef_Dx(edp,y,ordenDy)
if abs(K_)%1<casicero: # si es entero
    K_ = int(K_)
resultado['edp=0']  = edp
resultado['K_']  = K_

# lambda en ecuacion edp
lamb = (Dy**ordenDy)/(Dx**ordenDx)
if ordenDy==1 and ordenDx==2:
        lamb = lamb/K_
resultado['Lambda_L'] = lamb
# valor de Lambda en ecuacion edp
L_k = lamb.subs([(Dx,dx),(Dy,dy)])
if abs(L_k)%1<casicero: # si es entero
    L_k = int(L_k)
resultado['Lambda L_k'] = L_k

discreta = edp.lhs # EDP discreta
for derivada in dif_dividida:
    discreta = discreta.subs(derivada,dif_dividida[derivada])

# simplifica con Lambda L
discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
resultado['(discreta=0)*Dy**ordeny'] = discreta_L

discreta_L = edp_sustituye_L(resultado)
discreta_L = sym.collect(discreta_L,U) # agrupa u[,]
discreta_L = sym.Eq(discreta_L,0)
resultado['discreta_L=0'] = discreta_L
        
# sustituye constantes en ecuación a iterar
discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),(L,L_k)])
discreta_L = redondea_coef(discreta_L)
resultado['discreta'] = discreta_L

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Parabólica contínua a discreta')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    else:
        print('\n',entrada,':')
        sym.pprint(resultado[entrada])

# Método implícito EDP Parabólica ------------
# ITERAR para cada i,j dentro de U 
# x[i] , y[j]  valor en posición en cada eje
xi = np.arange(a,b+dx/2,dx)
yj = np.arange(y0,n_iteraciones*dy+dy/2,dy)
n = len(xi)
m = len(yj)
# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u
# llena u con valores en fronteras
u_xy[:,0]   = fxc(xi)  # inferior  Tc
u_xy[0,:]   = fya(yj)  # izquierda Ta
u_xy[n-1,:] = fyb(yj)  # derecha   Tb
u0 = np.copy(u_xy)     # matriz u inicial

# u_mask[i,j] con valores iniciales o calculados:True
u_mask = np.zeros(shape=(n,m),dtype=bool)
u_mask[0,:]  = True # izquierda Ta
u_mask[-1,:] = True # derecha   Tb
u_mask[:,0]  = True # inferior  Tc

# Método implícito para EDP Parabólica
# a usar en iteraciones
discreta = resultado['discreta']

def edp_sustituyeValorU(discreta,xi,yj,u_xy,u_mask):
    '''Sustituye en edp discreta los valores conocidos de U[i,j]
    tomados desde u_xy, marcados con u_mask
    u_mask indica si el valor se ha calculado con edp.
    '''
    LHS = discreta.lhs # lado izquierdo de ecuacion
    RHS = discreta.rhs # lado derecho
    # sustituye U[i,j] con valores conocidos
    A_diagonal = [] # lista de i,j para matriz de coeficientes A
    # Separa términos suma
    term_suma = sym.Add.make_args(LHS)
    for term_k in term_suma:
        # busca U[i,j] y cambia por valor uxt[i,j]
        cambiar = False ; cambiar_valor = 0 ; cambiar_factor = 0
        # separa cada factor de término
        factor_Mul = sym.Mul.make_args(term_k)
        for factor_k in factor_Mul:
            # busca U[i,j] en matriz uxt[i,j]
            if factor_k.is_Function:
                [i_k,j_k] = factor_k.args
                if not(np.isnan(u_xy[i_k,j_k])):
                    cambiar = u_mask[i_k,j_k]
                    cambiar_factor = factor_k
                    cambiar_valor = u_xy[i_k,j_k]
                else:
                    A_diagonal.append([i_k,j_k,term_k/factor_k])
        # encontró valor U[i,j],term_k va a la derecha de ecuación
        if cambiar == True:
            LHS = LHS - term_k
            term_ki = term_k.subs(cambiar_factor,cambiar_valor)
            RHS = RHS - term_ki
    discreta = sym.Eq(LHS,RHS)
    B_diagonal = RHS
    resultado = [discreta,A_diagonal,B_diagonal]
    return (resultado)

# ITERAR para plantear las ecuaciones en [i,j]
resultado = {} # resultados en diccionario
eq_itera = [] ; tamano = (n-2)
A = np.zeros(shape=(tamano,tamano),dtype=float)
B = np.zeros(tamano,dtype=float)
for j_k in range(1,m,1): # no usar valores en bordes
    resultado[j_k] ={}
    for i_k in range(1,n-1,1): 
        eq_conteo = i_k-1
        discreta_ij = discreta.subs({i:i_k,j:j_k,
                                     x:xi[i_k],y:yj[j_k]})
        resultado[j_k][eq_conteo]= {'j':j_k, 'i':i_k,
                                 'discreta_ij': discreta_ij}
        # usa valores de frontera segun u_mask con True
        discreta_k = edp_sustituyeValorU(discreta_ij,
                                        xi,yj,u_xy,u_mask)
        discreta_ij = discreta_k[0]
        A_diagonal  = discreta_k[1] # lista de (i,j,coeficiente) 
        B_diagonal  = discreta_k[2] # constante para B
        resultado[j_k][eq_conteo]['discreta_k'] = discreta_k[0]
        # añade ecuacion a resolver
        eq_itera.append(discreta_ij)
        # Aplica coeficientes de ecuacion en A y B:
        # A_diagonal tiene lista de (i,j,coeficiente) 
        for uncoeff in A_diagonal:
            columna = uncoeff[0]-1
            fila = eq_conteo
            A[fila,columna] = uncoeff[2]
        B[eq_conteo] = float(B_diagonal) # discreta_ij.rhs
    # resuelve el sistema A.X=B en Numpy
    X_k = np.linalg.solve(A,B)
    # actualiza uxt[i,j] , u_mask segun X_k
    u_xy[1:n-1,j_k] = X_k
    u_mask[1:n-1,j_k] = True
    # almacena resultados
    resultado[j_k]['A'] = np.copy(A)
    resultado[j_k]['B'] = np.copy(B)
    resultado[j_k]['X'] = np.copy(X_k)

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('\nMétodo implícito EDP Parabólica')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    elif (tipo==dict):
        if entrada<2:
            eq_conteo = resultado[entrada].keys()
            for una_eq in eq_conteo:
                if una_eq=='A' or una_eq=='B' or una_eq=='X':
                    print(una_eq,':')
                    print(resultado[entrada][una_eq])
                else:
                    print('j:',resultado[entrada][una_eq]['j'],'; '
                          'i:',resultado[entrada][una_eq]['i'])
                    sym.pprint(resultado[entrada][una_eq]['discreta_ij'])
                    sym.pprint(resultado[entrada][una_eq]['discreta_k'])
    else:
        print('\n',entrada,': \n',resultado[entrada])
print('Resultados para U(x,y)')
print('xi:',xi)
print('yj:',yj)
print(' j, U[i,j]')
for j_k in range(m-1,-1,-1):
    print(j_k, (u_xy[:,j_k]))

EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]

2023-10-052025-09-07

7.1.3 EDP Parabólica - analítico explícito con Sympy-Python

EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]
..

1. Ejercicio

Referencia: Chapra 30.2 p888 pdf912, Burden 9Ed 12.2 p725, Rodríguez 10.2 p406

\frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}

Valores de frontera: Ta = 60, Tb = 40, T0 = 25
Longitud en x a = 0, b = 1, Constante K= 4
Tamaño de paso dx = 0.1, dt = dx/10

EDP Parabólica [ contínua a discreta ] || [ sympy_explícito ] [sympy_ implícito ]
..

2. EDP Parabólica contínua a discreta

El resultado del algoritmo para el ejercicio propuesto es:

EDP Parabólica contínua a discreta
 ordenDx : 2
 ordenDy : 1

 edp=0 :
  2                             
 ∂               ∂              
───(u(x, y)) - 4⋅──(u(x, y)) = 0
  2              ∂y             
∂x                              
 K_ : 4

 Lambda_L :
  Dy 
─────
    2
4⋅Dx 
 Lambda L_k : 0.250000000000000

 (discreta=0)*Dy**ordeny :
                            2⋅Dy⋅U(i, j)   Dy⋅U(i - 1, j)   Dy⋅U(i + 1, j)
4⋅U(i, j) - 4⋅U(i, j + 1) - ──────────── + ────────────── + ──────────────
                                  2               2                2      
                                Dx              Dx               Dx       

 discreta_L=0 :
⎛    2⋅Dy⎞                           Dy⋅U(i - 1, j)   Dy⋅U(i + 1, j)    
⎜4 - ────⎟⋅U(i, j) - 4⋅U(i, j + 1) + ────────────── + ────────────── = 0
⎜      2 ⎟                                  2                2          
⎝    Dx  ⎠                                Dx               Dx           

 discreta :
2⋅U(i, j) - 4⋅U(i, j + 1) + U(i - 1, j) + U(i + 1, j) = 0

El desarrollo analítico comienza en escribir la ecuación diferencial parcial parabólica de la forma contínua a la representación discreta, siguiendo los criterios indicados en la parte teórica usando derivadas expresadas en diferencias divididas.

La EDP se escribe en formato Sympy en el bloque de ingreso en dos partes: lado izquierdo (LHS) de la ecuación y lado derecho (RHS), definiendo u como una función de variables x,y.

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion complementaria
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) 
RHS = 4*sym.diff(u,y,1) + fxy
edp = sym.Eq(LHS-RHS,0)

La ecuación f(x,y) se añade para ser usada con otros ejercicios que puede revisar en la sección de Evaluaciones.

Las expresiones para diferencias divididas seleccionadas se ingresan como un diccionario:

# centrada, adelante 
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,1): (U.subs(j,j+1)-U)/Dy}

Las condiciones iniciales en los extremos a y b, y la temperatura en la barra, pueden ser funciones matemáticas, por lo que se propone usar el formato lambda de las variables x o y. En el ejercicio básico presentado en la parte teórica, la barra tiene temperatura constante, por lo que escribe como una constante mas cero por la variable independiente, esta forma facilita la evaluación de un vector xi por ejemplo, en caso de disponer una expresión diferente.

# Valores de frontera
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 40 +0*y  # derecha
fxc = lambda x: 25 +0*x  # inferior, función inicial

los demás parámetros del ejercicio se ingresan más adelante de forma semejante al ejercicio con Numpy.

Al desarrollar el método explícito se requiere indicar el nodo u(i,j+1) a buscar

buscar = U.subs(j,j+1) # U[i,j+1] método explícito

El valor de K se puede determinar en la expresión al reordenar para que el término de Dx de mayor grado sea la unidad. El valor de K será el del término Dy con signo negativo, al encontrarse en el lado izquierdo (LHS) de la expresión, en la propuesta teórica se encontraba como positivo del lado derecho.

El término de mayor grado se puede obtener al revisar la expresión por cada término de la suma, buscando el factor que contiene Dx o Dy en cada caso. Por ser un poco ordenado se añade la función edp_coef_Dx().

Instrucciones en Python

Para revisar los resultados de algunos pasos para obtener la expresión discreta de EDP, se usa un diccionario a ser usado en el bloque de salida. Para las expresiones de ecuación se usa sym.pprint().

# Ecuaciones Diferenciales Parciales Parabólicas
# EDP Parabólicas contínua a discreta con Sympy
import numpy as np
import sympy as sym

# u(x,y) funciones continuas y variables simbólicas usadas
t = sym.Symbol('t',real=True) # independiente
x = sym.Symbol('x',real=True)
y = sym.Symbol('y',real=True) 
u = sym.Function('u')(x,y) # funcion 
f = sym.Function('f')(x,y) # funcion complemento
K = sym.Symbol('K',real=True)
# U[i,j] funciones discretas y variables simbólicas usadas
i  = sym.Symbol('i',integer=True,positive=True) # indices
j  = sym.Symbol('j',integer=True,positive=True)
Dt = sym.Symbol('Dt',real=True,positive=True)
Dx = sym.Symbol('Dx',real=True,positive=True)
Dy = sym.Symbol('Dy',real=True,positive=True)
L  = sym.Symbol('L',real=True)
U  = sym.Function('U')(i,j)

# INGRESO
fxy = 0*x+0*y  # f(x,y) = 0 , ecuacion complementaria
# ecuacion edp : LHS=RHS
LHS = sym.diff(u,x,2) 
RHS = 4*sym.diff(u,y,1) + fxy
edp = sym.Eq(LHS-RHS,0)

# centrada, adelante 
dif_dividida ={sym.diff(u,x,2): (U.subs(i,i-1)-2*U+U.subs(i,i+1))/(Dx**2),
               sym.diff(u,y,1): (U.subs(j,j+1)-U)/Dy}

buscar = U.subs(j,j+1) # U[i,j+1] método explícito

# Valores de frontera
fya = lambda y: 60 +0*y  # izquierda
fyb = lambda y: 40 +0*y  # derecha
fxc = lambda x: 25 +0*x  # inferior, función inicial

# [a,b] dimensiones de la barra
a = 0  # longitud en x
b = 1
y0 = 0 # tiempo inicial, aumenta con dt en n iteraciones

# muestreo en ejes, discreto
x_tramos = 10
dx  = (b-a)/x_tramos
dy  = dx/10
n_iteraciones = 4 # iteraciones en tiempo

verdigitos = 2      # para mostrar en tabla de resultados
casicero = 1e-15    # para redondeo de términos en ecuacion

# PROCEDIMIENTO --------------------------------------

def edp_coef_Dx(edp,x,ordenx):
    ''' Extrae el coeficiente de la derivada Dx de ordenx,
    edp es la ecuación como lhs=rhs
    '''
    coeff_x = 1.0 # valor inicial
    # separa cada término de suma
    term_suma = sym.Add.make_args(edp.lhs)
    for term_k in term_suma:
        if term_k.is_Mul: # mas de un factor
            factor_Mul = sym.Mul.make_args(term_k)
            coeff_temp = 1; coeff_usar = False
            # separa cada factor de término 
            for factor_k in factor_Mul:
                if not(factor_k.is_Derivative):
                    coeff_temp = coeff_temp*factor_k
                else: # factor con derivada de ordenx
                    partes = factor_k.args
                    if partes[1]==(x,ordenx):
                        coeff_usar = True
            if coeff_usar==True:
                coeff_x = coeff_x*coeff_temp
    return(coeff_x)

def redondea_coef(ecuacion, precision=6,casicero = 1e-15):
    ''' redondea coeficientes de términos suma de una ecuacion
    ecuación como lhs=rhs
    '''
    tipo = type(ecuacion)
    tipo_eq = False
    if tipo == sym.core.relational.Equality:
        RHS = ecuacion.rhs  # separa lado derecho
        ecuacion = ecuacion.lhs # analiza lado izquierdo
        tipo = type(ecuacion)
        tipo_eq = True

    if tipo == sym.core.add.Add: # términos suma de ecuacion
        term_sum = sym.Add.make_args(ecuacion)
        ecuacion = sym.S.Zero  # vacia
        for term_k in term_sum:
            # factor multiplicativo de termino suma
            term_mul = sym.Mul.make_args(term_k)
            producto = sym.S.One
            for factor in term_mul:
                if not(factor.has(sym.Symbol)): # es numerico
                    factor = np.around(float(factor),precision)
                    if (abs(factor)%1)<casicero: # si es entero
                        factor = int(factor)
                producto = producto*factor
            ecuacion = ecuacion + producto
            
    if tipo == sym.core.mul.Mul: # termino único, busca factores
        term_mul = sym.Mul.make_args(ecuacion)
        producto = sym.S.One
        for factor in term_mul:
            if not(factor.has(sym.Symbol)): # es numerico
                factor = np.around(float(factor),precision)
                if (abs(factor)%1)<casicero: # si es entero
                    factor = int(factor)
            producto = producto*factor
        ecuacion = producto
        
    if tipo == float: # solo un numero
        if (abs(ecuacion)%1)<casicero: 
            ecuacion = int(ecuacion)
            
    if tipo_eq==True: # era igualdad, integra lhs=rhs
        ecuacion = sym.Eq(ecuacion,RHS)

    return(ecuacion)


resultado = {} # resultados en diccionario
# orden de derivada por x, y
edp_x = edp.subs(x,0)
edp_y = edp.subs(y,0)
ordenDx = sym.ode_order(edp_x,u)
ordenDy = sym.ode_order(edp_y,u)
resultado['ordenDx'] = ordenDx  # guarda en resultados
resultado['ordenDy'] = ordenDy
# coeficiente derivada orden mayor a 1 (d2u/dx2)
coeff_x = edp_coef_Dx(edp,x,ordenDx)
LHS = edp.lhs  # lado izquierdo de edp
RHS = edp.rhs  # lado derecho de edp
if not(coeff_x==1):
    LHS = LHS/coeff_x
    RHS = RHS/coeff_x
# toda la expresión a la izquierda
edp = sym.Eq(LHS-RHS,0) 
K_ = -edp_coef_Dx(edp,y,ordenDy)
if abs(K_)%1<casicero: # si es entero
    K_ = int(K_)
resultado['edp=0']  = edp
resultado['K_']  = K_

# lambda en ecuacion edp
lamb = (Dy**ordenDy)/(Dx**ordenDx)
if ordenDy==1 and ordenDx==2:
        lamb = lamb/K_
resultado['Lambda_L'] = lamb
# valor de Lambda en ecuacion edp
L_k = lamb.subs([(Dx,dx),(Dy,dy)])
if abs(L_k)%1<casicero:
    L_k = int(L_k)
resultado['Lambda L_k'] = L_k

discreta = edp.lhs # EDP discreta
for derivada in dif_dividida: # reemplaza diferencia dividida
    discreta = discreta.subs(derivada,dif_dividida[derivada])

# simplifica con Lambda_L
discreta_L = sym.expand(discreta*(Dy**ordenDy),mul=True)
resultado['(discreta=0)*Dy**ordeny'] = discreta_L

discreta_L = discreta_L.subs(lamb,L)
discreta_L = sym.collect(discreta_L,U) #agrupa u[,]
discreta_L = sym.Eq(discreta_L,0)
resultado['discreta_L=0'] = discreta_L
        
# sustituye constantes en ecuación a iterar
discreta_L = discreta_L.subs([(Dx,dx),(Dy,dy),(L,L_k)])
discreta_L = redondea_coef(discreta_L)
resultado['discreta'] = discreta_L

# SALIDA
np.set_printoptions(precision=verdigitos)
algun_numero = [int,float,str,'Lambda L_k']
print('EDP Parabólica contínua a discreta')
for entrada in resultado:
    tipo = type(resultado[entrada])
    if tipo in algun_numero or entrada in algun_numero:
        print('',entrada,':',resultado[entrada])
    else:
        print('\n',entrada,':')
        sym.pprint(resultado[entrada])

EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]

3. EDP Parabólica - Método explícito con Sympy

El desarrollo del método explícito para una ecuación diferencial parcial parabólica usando Sympy, requiere añadir instrucciones para procesar la expresión obtenida en el algoritmo anterior con los valores de i,j.

Se añaden las instrucciones para iterar la expresión discreta para cada i,j dentro de la matriz U creada para el ejercicio. Se aplican las condiciones iniciales a los bordes de la matriz, para luego proceder a iterar.

Como en las actividades del curso se requiere realizar al menos 3 iteraciones para las expresiones del algoritmo con papel y lápiz, solo se presentarán los resultados de las expresiones para la primera fila en j=0.

El algoritmo para la EDP del ejemplo muestra el siguiente resultado:

Método explícito EDP Parabólica
discreta_itera:
              U(i, j)   U(i - 1, j)   U(i + 1, j)
U(i, j + 1) = ─────── + ─────────── + ───────────
                 2           4             4     

 iterar i,j:
j:  0  ;  i:  1  ;  ecuacion: 0
          U(0, 0)   U(1, 0)   U(2, 0)
U(1, 1) = ─────── + ─────── + ───────
             4         2         4   
U(1, 1) = 33.75
j:  0  ;  i:  2  ;  ecuacion: 1
          U(1, 0)   U(2, 0)   U(3, 0)
U(2, 1) = ─────── + ─────── + ───────
             4         2         4   
U(2, 1) = 25.0
j:  0  ;  i:  3  ;  ecuacion: 2
          U(2, 0)   U(3, 0)   U(4, 0)
U(3, 1) = ─────── + ─────── + ───────
             4         2         4   
U(3, 1) = 25.0
j:  0  ;  i:  4  ;  ecuacion: 3
          U(3, 0)   U(4, 0)   U(5, 0)
U(4, 1) = ─────── + ─────── + ───────
             4         2         4   
U(4, 1) = 25.0
j:  0  ;  i:  5  ;  ecuacion: 4
          U(4, 0)   U(5, 0)   U(6, 0)
U(5, 1) = ─────── + ─────── + ───────
             4         2         4   
U(5, 1) = 25.0
j:  0  ;  i:  6  ;  ecuacion: 5
          U(5, 0)   U(6, 0)   U(7, 0)
U(6, 1) = ─────── + ─────── + ───────
             4         2         4   
U(6, 1) = 25.0
j:  0  ;  i:  7  ;  ecuacion: 6
          U(6, 0)   U(7, 0)   U(8, 0)
U(7, 1) = ─────── + ─────── + ───────
             4         2         4   
U(7, 1) = 25.0
j:  0  ;  i:  8  ;  ecuacion: 7
          U(7, 0)   U(8, 0)   U(9, 0)
U(8, 1) = ─────── + ─────── + ───────
             4         2         4   
U(8, 1) = 25.0
j:  0  ;  i:  9  ;  ecuacion: 8
          U(8, 0)   U(9, 0)   U(10, 0)
U(9, 1) = ─────── + ─────── + ────────
             4         2         4    
U(9, 1) = 28.75
... continua calculando

x: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0.   0.01 0.02 0.03 0.04]
 j, U[i,j]
4 [60.   42.77 31.29 26.37 25.14 25.   25.06 25.59 27.7  32.62 40.  ]
3 [60.   40.86 29.38 25.55 25.   25.   25.   25.23 26.88 31.8  40.  ]
2 [60.   38.12 27.19 25.   25.   25.   25.   25.   25.94 30.62 40.  ]
1 [60.   33.75 25.   25.   25.   25.   25.   25.   25.   28.75 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]

Instrucciones en Python

# Método explícito EDP Parabólica ------------
# ITERAR para cada i,j dentro de U 
# x[i] , y[j]  valor en posición en cada eje
xi = np.arange(a,b+dx/2,dx)
yj = np.arange(y0,n_iteraciones*dy+dy,dy)
n = len(xi)
m = len(yj)
# Matriz U
u_xy = np.zeros(shape=(n,m),dtype = float)
u_xy = u_xy*np.nan # valor inicial dentro de u
# llena u con valores en fronteras
u_xy[:,0]   = fxc(xi)  # inferior  Tc
u_xy[0,:]   = fya(yj)  # izquierda Ta
u_xy[n-1,:] = fyb(yj)  # derecha   Tb
u0 = np.copy(u_xy)     # matriz u inicial

# busca el término no conocido
u_explicito = sym.solve(discreta_L,buscar)
u_explicito = sym.Eq(buscar,u_explicito[0])
resultado['discreta_itera'] = u_explicito

print('\nMétodo explícito EDP Parabólica')
print('discreta_itera:')
sym.pprint(u_explicito)

# iterando
print('\n iterar i,j:')
j_k = 0
for j_k in range(0,m-1,1):
    for i_k in range(1,n-1,1):
        eq_conteo = j_k*(n-2)+(i_k-1)
        discreta_ij = u_explicito.subs({i:i_k,j:j_k})
        valorij = 0 # calcula valor de nodo ij
        # separa cada término de suma
        term_suma = sym.Add.make_args(discreta_ij.rhs)
        for term_k in term_suma:
            term_ki = 1
            if term_k.is_Function:
                [i_f,j_f] = term_k.args
                term_ki = u_xy[i_f,j_f]
            elif term_k.is_Mul: # mas de un factor
                factor_Mul = sym.Mul.make_args(term_k)
                # separa cada factor de término 
                for factor_k in factor_Mul:
                    if factor_k.is_Function:
                        [i_f,j_f] = factor_k.args
                        term_ki = term_ki*u_xy[i_f,j_f]
                    else:
                        term_ki = term_ki*factor_k
            else:
                term_ki = term_k
            valorij = valorij + term_ki
        [i_f,j_f] = discreta_ij.lhs.args
        u_xy[i_f,j_f] = float(valorij) # actualiza matriz u
        # muestra las primeras m iteraciones
        if eq_conteo<(n-2):
            print('j: ',j_k,' ; ','i: ',i_k,
                  ' ; ','ecuacion:',eq_conteo)
            sym.pprint(discreta_ij)
            print(discreta_ij.lhs,'=',float(valorij))
print('... continua calculando')

# SALIDA
np.set_printoptions(precision=verdigitos)
print('\nx:',xi)
print('t:',yj)
print(' j, U[i,j]')
for j_k in range(m-1,-1,-1):
    print(j_k, (u_xy[:,j_k]))

EDP Parabólica [ contínua a discreta ] || sympy: [ explícito ] [ implícito ]

2017-10-082025-09-07

7.3 EDP hiperbólicas con Python

[ concepto ] [ analítico ] [ algoritmo ]
..

1. Concepto y ejercicio

Referencia: Chapra PT8.1 p860, Rodríguez 10.4 p435

Las Ecuaciones Diferenciales Parciales tipo hiperbólicas semejantes a la mostrada, para un ejemplo en particular, representa la ecuación de onda de una dimensión u[x,t], que representa el desplazamiento vertical de una cuerda.

\frac{\partial ^2 u}{\partial t^2}=c^2\frac{\partial ^2 u}{\partial x^2}

Los extremos de la cuerda de longitud unitaria están sujetos y referenciados a una posición 0 a la izquierda y 1 a la derecha.

u[x,t] , 0<x<1, t≥0
u(0,t) = 0 , t≥0
u(1,t) = 0 , t≥0

Al inicio, la cuerda está estirada por su punto central:

u(x,0) = \begin{cases} -0.5x &, 0\lt x\leq 0.5 \\ 0.5(x-1) &, 0.5\lt x \lt 1 \end{cases}

Se suelta la cuerda, con velocidad cero desde la posición inicial:

\frac{\partial u(x,0)}{\partial t} = 0

[ concepto ] [ analítico ] [ algoritmo ]
..

2. Desarrollo analítico

La solución se realiza de forma semejante al procedimiento para EDP parabólicas y elípticas. Se cambia a la forma discreta usando diferencias finitas divididas:

\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{(\Delta t)^2} =c^2 \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2}

El error es del orden $O(\Delta x)^2 + O(\Delta t)^2$ .
se reagrupa de la forma:

u_{i,j+1}-2u_{i,j}+u_{i,j-1} = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \big( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \big)

El término constante, por facilidad se simplifica con el valor de 1

\lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} =1

si c = 2 y Δx = 0.2, se deduce que Δt = 0.1

que al sustituir en la ecuación, se simplifica anulando el término u[i,j]:

u_{i,j+1}+u_{i,j-1} = u_{i+1,j}+u_{i-1,j}

en los que intervienen solo los puntos extremos. Despejando el punto superior del rombo:

u_{i,j+1} = u_{i+1,j}-u_{i,j-1}+u_{i-1,j}

Convergencia:

\lambda = \frac{c^2 (\Delta t)^2}{(\Delta x)^2} \leq 1

para simplificar aún más el problema, se usa la condición velocidad inicial de la cuerda igual a cero

\frac{\delta u_{i,0}}{\delta t}=\frac{u_{i,1}-u_{i,-1}}{2\Delta t} = 0

u_{i,-1}=u_{i,1}

que se usa para cuando el tiempo es cero, permite calcular los puntos para t[1]:

j=0

u_{i,1} = u_{i+1,0}-u_{i,-1}+u_{i-1,0}

u_{i,1} = u_{i+1,0}-u_{i,1}+u_{i-1,0}

2 u_{i,1} = u_{i+1,0}+u_{i-1,0}

u_{i,1} = \frac{u_{i+1,0}+u_{i-1,0}}{2}

Aplicando solo cuando j = 0

que al ponerlos en la malla se logra un sistema explícito para cada u[i,j], con lo que se puede resolver con un algoritmo:

[ concepto ] [ analítico ] [ algoritmo ]
..

3. Algoritmo con Python

# Ecuaciones Diferenciales Parciales
# Hiperbólica. Método explicito
import numpy as np

def cuerdainicio(xi):
    n = len(xi)
    y = np.zeros(n,dtype=float)
    for i in range(0,n,1):
        if (xi[i]<=0.5):
            y[i]=-0.5*xi[i]
        else:
            y[i]=0.5*(xi[i]-1)
    return(y)

# INGRESO
x0 = 0 # Longitud de cuerda
xn = 1
y0 = 0 # Puntos de amarre
yn = 0
t0 = 0 # tiempo inicial
c = 2  # constante EDP
# discretiza
tramosx = 16
tramost = 32
dx = (xn-x0)/tramosx 
dt = dx/c

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx/2,dx)
tj = np.arange(0,tramost*dt+dt/2,dt)
n = len(xi)
m = len(tj)

u = np.zeros(shape=(n,m),dtype=float)
u[:,0] = cuerdainicio(xi)
u[0,:] = y0
u[n-1,:] = yn
# Aplicando condición inicial
j = 0
for i in range(1,n-1,1):
    u[i,j+1] = (u[i+1,j]+u[i-1,j])/2
# Para los otros puntos
for j in range(1,m-1,1):
    for i in range(1,n-1,1):
        u[i,j+1] = u[i+1,j]-u[i,j-1]+u[i-1,j]

# SALIDA
np.set_printoptions(precision=2)
print('xi =')
print(xi)
print('tj =')
print(tj)
print('matriz u =')
print(u)

con lo que se obtienen los resultados numéricos, que para mejor interpretación se presentan en una gráfica estática y otra animada.

# GRAFICA
import matplotlib.pyplot as plt

for j in range(0,m,1):
    y = u[:,j]
    plt.plot(xi,y)

plt.title('EDP hiperbólica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

# **** GRÁFICO CON ANIMACION ***********
import matplotlib.animation as animation

# Inicializa parametros de trama/foto
retardo = 70   # milisegundos entre tramas
tramas  = m
maximoy = np.max(np.abs(u))
figura, ejes = plt.subplots()
plt.xlim([x0,xn])
plt.ylim([-maximoy,maximoy])

# lineas de función y polinomio en gráfico
linea_poli, = ejes.plot(xi,u[:,0], '-')
plt.axhline(0, color='k')  # Eje en x=0

plt.title('EDP hiperbólica')
# plt.legend()
# txt_x = (x0+xn)/2
txt_y = maximoy*(1-0.1)
texto = ejes.text(x0,txt_y,'tiempo:',
                  horizontalalignment='left')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()

# Nueva Trama
def unatrama(i,xi,u):
    # actualiza cada linea
    linea_poli.set_ydata(u[:,i])
    linea_poli.set_xdata(xi)
    linea_poli.set_label('tiempo linea: '+str(i))
    texto.set_text('tiempo['+str(i)+']')
    # color de la línea
    if (i<=9):
        lineacolor = 'C'+str(i)
    else:
        numcolor = i%10
        lineacolor = 'C'+str(numcolor)
    linea_poli.set_color(lineacolor)
    return linea_poli, texto

# Limpia Trama anterior
def limpiatrama():
    linea_poli.set_ydata(np.ma.array(xi, mask=True))
    linea_poli.set_label('')
    texto.set_text('')
    return linea_poli, texto

# Trama contador
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(figura,
                              unatrama,
                              i ,
                              fargs=(xi,u),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo video y GIFAnimado :
# ani.save('EDP_hiperbólica.mp4')
ani.save('EDP_hiperbolica.gif', writer='imagemagick')
plt.draw()
plt.show()

Una vez realizado el algoritmo, se pueden cambiar las condiciones iniciales de la cuerda y se observan los resultados.

Se recomienda realizar otros ejercicios de la sección de evaluaciones para EDP Hiperbólicas y observar los resultados con el algoritmo modificado.

[ concepto ] [ analítico ] [ algoritmo ]

2017-10-082025-09-07

7.2.2 EDP Elípticas método implícito con Python

EDP Elípticas [ concepto ] Método implícito: [ Analítico ] [ Algoritmo ]
..

1. EDP Elípticas: Método Implícito – Desarrollo Analítico

con el resultado desarrollado en EDP elípticas para:

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0

y 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} = 0

con lo que para el método implícito, se plantea un sistema de ecuaciones para determinar los valores en cada punto desconocido.

j=1, i =1

u_{2,1}-4u_{1,1}+u_{0,1} + u_{1,2} +u_{1,0} = 0

u_{2,1}-4u_{1,1}+Ta + u_{1,2} +Tc= 0

-4u_{1,1}+u_{2,1}+u_{1,2} = -(Tc+Ta)

j=1, i =2

u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +u_{2,0} = 0

u_{3,1}-4u_{2,1}+u_{1,1} + u_{2,2} +Tc = 0

u_{1,1}-4u_{2,1}+u_{3,1}+ u_{2,2}= -Tc

j=1, i=3

u_{4,1}-4u_{3,1}+u_{2,1} + u_{3,2} +u_{3,0} = 0

Tb-4u_{3,1}+u_{2,1} + u_{3,2} +Tc = 0

u_{2,1} -4u_{3,1} + u_{3,2} = -(Tc+Tb)

j=2, i=1

u_{2,2}-4u_{1,2}+u_{0,2} + u_{1,3} +u_{1,1} = 0

u_{2,2}-4u_{1,2}+Ta + u_{1,3} +u_{1,1} = 0

-4u_{1,2}+u_{2,2}+u_{1,1}+u_{1,3} = -Ta

j = 2, i = 2

u_{1,2}-4u_{2,2}+u_{3,2} + u_{2,3} +u_{2,1} = 0

j = 2, i = 3

u_{4,2}-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0

Tb-4u_{3,2}+u_{2,2} + u_{3,3} +u_{3,1} = 0

u_{2,2} -4u_{3,2}+ u_{3,3} +u_{3,1} = -Tb

j=3, i = 1

u_{2,3}-4u_{1,3}+u_{0,3} + u_{1,4} +u_{1,2} = 0

u_{2,3}-4u_{1,3}+Ta + Td +u_{1,2} = 0

-4u_{1,3}+u_{2,3}+u_{1,2} = -(Td+Ta)

j=3, i = 2

u_{3,3}-4u_{2,3}+u_{1,3} + u_{2,4} +u_{2,2} = 0

u_{3,3}-4u_{2,3}+u_{1,3} + Td +u_{2,2} = 0

+u_{1,3} -4u_{2,3}+u_{3,3} +u_{2,2} = -Td

j=3, i=3

u_{4,3}-4u_{3,3}+u_{2,3} + u_{3,4} +u_{3,2} = 0

Tb-4u_{3,3}+u_{2,3} + Td +u_{3,2} = 0

u_{2,3}-4u_{3,3}+u_{3,2} = -(Td+Tb)

con las ecuaciones se arma una matriz:

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([-(Tc+Ta),-Tc,-(Tc+Tb),
                  -Ta,   0,    -Tb,
              -(Td+Ta),-Td,-(Td+Tb)])

que al resolver el sistema de ecuaciones se obtiene:

>>> Xu
array([ 56.43,  55.71,  56.43,  60.  ,  60.  ,  60.  ,  63.57,  64.29,
        63.57])

ingresando los resultados a la matriz u:

xi=
[ 0.   0.5  1.   1.5  2. ]
yj=
[ 0.    0.38  0.75  1.12  1.5 ]
matriz u
[[ 60.    60.    60.    60.    60.  ]
 [ 50.    56.43  60.    63.57  70.  ]
 [ 50.    55.71  60.    64.29  70.  ]
 [ 50.    56.43  60.    63.57  70.  ]
 [ 60.    60.    60.    60.    60.  ]]
>>>

EDP Elípticas [ concepto ] Método implícito: [ Analítico ] [ Algoritmo ]

2. EDP Elípticas: Método Implícito – Desarrollo con algoritmo

Instrucciones en Python

# Ecuaciones Diferenciales Parciales
# Elipticas. Método implícito
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
tramosx = 4
tramosy = 4
dx = (xn-x0)/tramosx 
dy = (yn-y0)/tramosy 
maxitera = 100
tolera = 0.0001

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([-(Tc+Ta),-Tc,-(Tc+Tb),
              -Ta,0,-Tb,
              -(Td+Ta),-Td,-(Td+Tb)])


# PROCEDIMIENTO
# Resuelve sistema ecuaciones
Xu = np.linalg.solve(A,B)
[nx,mx] = np.shape(A)

xi = np.arange(x0,xn+dx/2,dx)
yj = np.arange(y0,yn+dy/2,dy)
n = len(xi)
m = len(yj)

u = np.zeros(shape=(n,m),dtype=float)
u[:,0]   = Tc
u[:,m-1] = Td
u[0,:]   = Ta
u[n-1,:] = Tb
u[1:1+3,1] = Xu[0:0+3]
u[1:1+3,2] = Xu[3:3+3]
u[1:1+3,3] = Xu[6:6+3]

# SALIDA
np.set_printoptions(precision=2)
print('xi=')
print(xi)
print('yj=')
print(yj)
print('matriz u')
print(u)

La gráfica de resultados se obtiene de forma semejante al ejercicio con método iterativo.

Se podría estandarizar un poco más el proceso para que sea realizado por el algoritmo y sea más sencillo generar la matriz con más puntos. Tarea.

2017-10-072025-09-06

7.2.1 EDP Elípticas método iterativo con Python

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}=0

y 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} = 0

que 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 ]

2017-10-062025-08-23

7.2 EDP Elípticas

EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]
..

1. EDP Elípticas

Referencia: Chapra 29.1 p866, Rodríguez 10.3 p425, Burden 12.1 p694

Las Ecuaciones Diferenciales Parciales tipo elípticas semejantes a la mostrada:

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{ \partial y^2} = 0

(ecuación de Laplace, Ecuación de Poisson con f(x,y)=0)

se interpreta como una placa metálica de dimensiones Lx y Ly, delgada con aislante que recubren las caras de la placa, y sometidas a condiciones en las fronteras:

Lx = dimensión x de placa metálica
Ly = dimensión y de placa metálica
u[0,y]  = Ta
u[Lx,y] = Tb
u[x,0]  = Tc
u[x,Ly] = Td

Para el planteamiento se usa una malla en la que cada nodo corresponden a los valores u[x_i,y_j]. Para simplificar la nomenclatura se usan los subíndices i para el eje de las x, j para el eje t, quedando u[i,j].

vista superior, plano xy:

La ecuación se cambia a la forma discreta, usando diferencias divididas que se sustituyen en la ecuación, por ejemplo:

\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}=0

Se agrupan los términos de los diferenciales:

\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 u_{i,j} + \lambda u_{i-1,j} + u_{i,j+1} - 2 u_{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}=0

Obteniendo así la solución numérica conceptual. La forma de resolver el problema determina el nombre del método a seguir.

Una forma de simplificar la expresión, es hacer λ=1, es decir $\lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = 1$ , se determina que los tamaño de paso Δx, Δy son iguales.

u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1} +u_{i,j-1} = 0

EDP Elípticas [ concepto ] Método [ iterativo ] [ implícito ]

2017-10-042025-09-07

7.1.2 EDP Parabólicas método implícito con Python

EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]

1. Ejercicio

Referencia: Chapra 30.3 p893 pdf917, Burden 9Ed 12.2 p729, Rodríguez 10.2.4 p415

Siguiendo el tema anterior en EDP Parabólicas, se tiene que:

\frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}

EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]

2. Desarrollo Analítico

En éste caso se usan diferencias finitas centradas y hacia atrás; la línea de referencia es t₁:

\frac{\partial^2 u}{\partial x^2} = \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2}

\frac{\partial u}{\partial t} = \frac{u_{i,j} - u_{i,j-1} }{\Delta t}

La selección de las diferencias divididas corresponden a los puntos que se quieren usar para el cálculo, se observan mejor en la gráfica de la malla.

Luego se sustituyen en la ecuación del problema, obteniendo:

\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = K\frac{u_{i,j}-u_{i,j-1}}{\Delta t}

De la gráfica se destaca que en la fórmula, dentro del triángulo solo hay DOS valores desconocidos, destacados por los punto en amarillo.
En la ecuación se representa por U[i,j] y U[i+1,j]. Por lo que será necesario crear un sistema de ecuaciones sobre toda la línea de tiempo t₁ para resolver el problema.

Despejando la ecuación, se agrupan términos constantes: $λ = \frac{\Delta t}{K (\Delta x)^2}$ .

\lambda u_{i-1,j} + (-1-2\lambda) u_{i,j} + \lambda u_{i+1,j} = -u_{i,j-1}

Los parámetro P, Q y R se determinan de forma semejante al método explícito:

P = \lambda

Q = -1-2\lambda

R = \lambda

Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j} = -u_{i,j-1}

Los valores en los extremos son conocidos, para los puntos intermedios se crea un sistema de ecuaciones para luego usar la forma Ax=B y resolver los valores para cada u(x_i,t_j).

Por ejemplo con cuatro tramos entre extremos se tiene que:
indice de tiempo es 1 e índice de x es 1.

i=1,j=1

Pu_{0,1} + Qu_{1,1} + Ru_{2,1} = -u_{1,0}

i=2,j=1

Pu_{1,1} + Qu_{2,1} + Ru_{3,1} = -u_{2,0}

i=3,j=1

Pu_{2,1} + Qu_{3,1} + Ru_{4,1} = -u_{3,0}

agrupando ecuaciones y sustituyendo valores conocidos:
$\begin{cases} Qu_{1,1} + Ru_{2,1} + 0 &= -T_{0}-PT_{A}\\Pu_{1,1} + Qu_{2,1} + Ru_{3,1} &= -T_{0}\\0+Pu_{2,1}+Qu_{3,1}&=-T_{0}-RT_{B}\end{cases}$

que genera la matriz a resolver:

\begin{bmatrix} Q && R && 0 && -T_{0}-PT_{A}\\P && Q && R && -T_{0}\\0 && P && Q && -T_{0}-RT_{B}\end{bmatrix}

Use alguno de los métodos de la unidad 3 para resolver el sistema y obtener los valores correspondientes.

Por la extensión de la solución es conveniente usar un algoritmo y convertir los pasos o partes pertinentes a funciones.

Tarea: Revisar y comparar con el método explícito.

EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]

Algoritmos en Python

Para la solución con el método implícito, se obtiene el mismo resultado en la gráfica y tabla. Aunque el algoritmo es diferente.

algunos valores:

iteración j: 1
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-40. -25. -25. -25. -25. -25. -25. -25. -35.]
resultado en j: 1 
 [31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57]
iteración j: 2
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-46.01 -26.03 -25.18 -25.03 -25.01 -25.01 -25.08 -25.44 -37.57]
resultado en j: 2 
 [35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39]
iteración j: 3
A:
 [[-1.5   0.25  0.    0.    0.    0.    0.    0.    0.  ]
 [ 0.25 -1.5   0.25  0.    0.    0.    0.    0.    0.  ]
 [ 0.    0.25 -1.5   0.25  0.    0.    0.    0.    0.  ]
 [ 0.    0.    0.25 -1.5   0.25  0.    0.    0.    0.  ]
 [ 0.    0.    0.    0.25 -1.5   0.25  0.    0.    0.  ]
 [ 0.    0.    0.    0.    0.25 -1.5   0.25  0.    0.  ]
 [ 0.    0.    0.    0.    0.    0.25 -1.5   0.25  0.  ]
 [ 0.    0.    0.    0.    0.    0.    0.25 -1.5   0.25]
 [ 0.    0.    0.    0.    0.    0.    0.    0.25 -1.5 ]]
B:
 [-50.25 -27.49 -25.55 -25.12 -25.03 -25.05 -25.24 -26.07 -39.39]
resultado en j: 3 
 [38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72]
EDP Parabólica - Método implícito 
lambda:  0.25
x: [0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
t: [0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 ] ... 1.0
Tabla de resultados en malla EDP Parabólica
j, U[:,: 10 ], primeras iteraciones
10 [60.   47.43 37.58 31.3  27.98 26.64 26.63 27.83 30.43 34.63 40.  ]
9 [60.   46.75 36.68 30.56 27.48 26.3  26.33 27.47 30.04 34.33 40.  ]
8 [60.   45.96 35.69 29.8  27.01 26.   26.06 27.12 29.6  33.99 40.  ]
7 [60.   45.01 34.6  29.02 26.57 25.73 25.81 26.76 29.13 33.58 40.  ]
6 [60.   43.87 33.39 28.24 26.16 25.51 25.58 26.41 28.6  33.09 40.  ]
5 [60.   42.45 32.06 27.48 25.8  25.32 25.39 26.07 28.03 32.48 40.  ]
4 [60.   40.67 30.61 26.75 25.51 25.18 25.24 25.76 27.41 31.71 40.  ]
3 [60.   38.34 29.06 26.09 25.28 25.09 25.13 25.47 26.74 30.72 40.  ]
2 [60.   35.25 27.49 25.55 25.12 25.03 25.05 25.24 26.07 29.39 40.  ]
1 [60.   31.01 26.03 25.18 25.03 25.01 25.01 25.08 25.44 27.57 40.  ]
0 [60. 25. 25. 25. 25. 25. 25. 25. 25. 25. 40.]
>>>

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método implícito
# Referencia: Chapra 30.3 p.895 pdf.917
#       Rodriguez 10.2.5 p.417
import numpy as np

# INGRESO
# Valores de frontera
Ta = 60
Tb = 40
#T0 = 25 # estado inicial de barra
fx = lambda x: 25 + 0*x # función inicial para T0
a = 0  # longitud en x
b = 1

K = 4     # Constante K
dx = 0.1  # Tamaño de paso
dt = dx/10

n = 100 # iteraciones en tiempo
verdigitos = 2   # decimales a mostrar en tabla de resultados

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx/2,dx)
fi = fx(xi)
m  = len(xi)
ultimox = m-1
ultimot = n-1
# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j = 0
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = -1 -2*lamb
R = lamb

# Calcula U para cada tiempo + dt
tj = np.arange(0,(n+1)*dt,dt)
j = 1
while not(j>=n):
    # Matriz de ecuaciones
    k = m-2
    A = np.zeros(shape=(k,k), dtype = float)
    B = np.zeros(k, dtype = float)
    for f in range(0,k,1):
        if (f>0):
            A[f,f-1]=P
        A[f,f] = Q
        if (f<(k-1)):
            A[f,f+1]=R
        B[f] = -u[f+1,j-1]
    B[0] = B[0]-P*u[0,j]
    B[k-1] = B[k-1]-R*u[m-1,j]
    # Resuelve sistema de ecuaciones
    C = np.linalg.solve(A, B)

    # copia resultados a u[i,j]
    for f in range(0,k,1):
        u[f+1,j] = C[f]

    # siguiente iteración
    j = j + 1

    # muestra 3 primeras iteraciones
    np.set_printoptions(precision=verdigitos)
    if j<(3+2): 
        print('iteración j:',j-1)
        print('A:\n',A)
        print('B:\n',B)
        print('resultado en j:',j-1,'\n',C)
        
# SALIDA
print('EDP Parabólica - Método implícito ')
print('lambda: ',np.around(lamb,verdigitos))
print('x:',xi)
print('t:',tj[0:len(xi)],'...',tj[-1])
print('Tabla de resultados en malla EDP Parabólica')
print('j, U[:,:',ultimox,'], primeras iteraciones')
for j in range(ultimox,-1,-1):
    print(j,u[:,j])

Para realizar la gráfica se aplica lo mismo que en el método explícito

# GRAFICA ------------
import matplotlib.pyplot as plt
tramos = 10
salto = int(n/tramos) # evita muchas líneas
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.',color='red')
plt.xlabel('x[i]')
plt.ylabel('u[j]')
plt.title('Solución EDP parabólica')
plt.show()

Sin embargo para la gráfica en 3D se ajustan los puntos de referencia agregados por las diferencias finitas.

# GRAFICA en 3D
# vector de tiempo, simplificando en tramos
tramos = 10
salto = int(n/tramos)
if (salto == 0):
    salto = 1
tj = np.arange(0,n*dt,dt)
tk = np.zeros(tramos,dtype=float)

# Extrae parte de la matriz U,acorde a los tramos
U = np.zeros(shape=(tramos,m),dtype=float)
for k in range(0,tramos,1):
    U[k,:] = u[:,k*salto]
    tk[k] = tj[k*salto]
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi,tk)

fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U,
                       color ='blue',label='EDP Parabólica')
graf_3D.plot(Xi[2,0],Yi[2,0],U[2,0],'o',color ='orange')
graf_3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')
graf_3D.plot(Xi[2,2],Yi[2,2],U[2,2],'o',color ='salmon')
graf_3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='green')
graf_3D.set_title('EDP Parabólica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('t')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()

Queda por revisar la convergencia y estabilidad de la solución a partir de los O(h) de cada aproximación usada. Revisar los criterios en Chapra 30.2.1 p891, Burden 9Ed 12.2 p727, Rodríguez 10.2.2 p409 .

Tarea o proyecto: Realizar la comparación de tiempos de ejecución entre los métodos explícitos e implícitos. La parte de animación funciona igual en ambos métodos. Los tiempos de ejecución se determinan usando las instrucciones descritas en el enlace: Tiempos de Ejecución en Python

EDP Parabólica [ concepto ] [ ejercicio ]
Método implícito: [ Analítico ] [ Algoritmo ]