1. Ejercicio - EDP Parabólicas
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:
Para la ecuación diferencial parcial parabólica:
\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 |
2. Método implícito - paso a paso
En éste caso se usan diferencias finitas centradas y hacia atrás; la línea de referencia es t1:
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 t1 para resolver el problema.
Se agrupan términos constantes:
\frac{\Delta t}{K}\left( \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} \right) = \left(\cancel{K} \frac{u_{i,j}-u_{i,j-1}}{\cancel{\Delta t}}\right) \frac{\cancel{\Delta t}}{\cancel{K}} \frac{\Delta t}{K(\Delta x)^2}\left( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \right) = u_{i,j}-u_{i,j-1} λ = \frac{\Delta t}{K (\Delta x)^2}Se agrupan términos semejantes, dejando el valor conocido a la derecha de la expresión.
\lambda \left( u_{i+1,j}-2u_{i,j}+u_{i-1,j} \right) = u_{i,j}-u_{i,j-1} \lambda u_{i+1,j}-2\lambda u_{i,j}+ \lambda u_{i-1,j} = u_{i,j}-u_{i,j-1} \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:
3. Desarrollo analítico - paso a paso
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(xi,tj).
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} P T_{A} + Qu_{1,1} + Ru_{2,1} = -T_{0} Qu_{1,1} + Ru_{2,1} = -T_{0} - P T_{A}i=2,j=1
Pu_{1,1} + Qu_{2,1} + Ru_{3,1} = -u_{2,0} Pu_{1,1} + Qu_{2,1} + Ru_{3,1} = -T_{0}i=3,j=1
Pu_{2,1} + Qu_{3,1} + Ru_{4,1} = -u_{3,0} Pu_{2,1} + Qu_{3,1} + RT_{B} = -T_{0} Pu_{2,1} + Qu_{3,1} = -T_{0} - RT_{B}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.
4. 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] ... 1.0
Tabla de resultados en malla EDP Parabólica
j, U[:,: 5 ], primeras iteraciones de 100
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, 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
fc = 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), dtype=float)
# u[i,j], valores iniciales
u[0,:]= Ta # Izquierda
u[ultimox,:] = Tb # derecha
# estado inicial de barra en x, Tc
fci = fc(xi) # f(x) en tiempo inicial
u[1:ultimox,0] = fci[1:ultimox]
# Calcula U para cada tiempo + dt
tj = np.arange(t0,(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
j_mostrar = 5
np.set_printoptions(precision=verdigitos)
print('EDP Parabólica - Método implícito ')
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])
5. Gráfica con Python
Para realizar la gráfica se aplica lo mismo que en el método explícito
# GRAFICA ------------
import matplotlib.pyplot as plt
titulo = 'EDP parabólica - Método Implícito'
etiq_x = 'x[i]'
etiq_y = 't[j]'
etiq_z = 'u[i,j]'
tramos = 10
# GRAFICA 2D ------
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')
# entorno de gráfica
plt.xlabel(etiq_x)
plt.ylabel(etiq_z)
plt.title(titulo)
plt.tight_layout()
# plt.show() #comentar para la siguiente gráfica
Sin embargo para la gráfica en 3D se ajustan los puntos de referencia agregados por las diferencias finitas.
# 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)
fig3D = plt.figure()
graf3D = fig3D.add_subplot(111, projection='3d')
graf3D.plot_wireframe(Xi,Yi,U,
color ='blue')
graf3D.plot(Xi[2,0],Yi[2,0],U[2,0],'o',color ='orange')
graf3D.plot(Xi[2,1],Yi[2,1],U[2,1],'o',color ='salmon')
graf3D.plot(Xi[2,2],Yi[2,2],U[2,2],'o',color ='salmon')
graf3D.plot(Xi[1,1],Yi[1,1],U[1,1],'o',color ='green',
label='U[i,j-1]')
# entorno de gráfica
graf3D.set_title(titulo)
graf3D.set_xlabel(etiq_x)
graf3D.set_ylabel(etiq_y)
graf3D.set_zlabel(etiq_z)
graf3D.legend()
graf3D.view_init(35, -45)
plt.tight_layout()
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