Referencia: Chapra 30.3 p893 pdf917, Burden 9Ed 12.2 p729, Rodriguez 10.2.4 p415
Se utiliza el mismo ejercicio presentado en método explícito.
\frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}En éste caso se usan diferencias finitas centradas y hacia atras; la línea de referencia es t1:
\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 t1 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(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}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.
Algoritmos en Python
para la solución con el método implícito.
# 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 import matplotlib.pyplot as plt # INGRESO # Valores de frontera Ta = 60 Tb = 40 T0 = 25 # longitud en x a = 0 b = 1 # Constante K K = 4 # Tamaño de paso dx = 0.1 dt = 0.01 # iteraciones en tiempo n = 100 # PROCEDIMIENTO # Valores de x xi = np.arange(a,b+dx,dx) m = len(xi) ultimox = m-1 # Resultados en tabla de u u = np.zeros(shape=(m,n), dtype=float) # valores iniciales de u[:,j] j = 0 ultimot = n-1 u[0,j]= Ta u[1:ultimox,j] = T0 u[ultimox,j] = Tb # factores P,Q,R lamb = dt/(K*dx**2) P = lamb Q = -1 -2*lamb R = lamb # Calcula U para cada tiempo + dt j = 1 while not(j>=n): u[0,j] = Ta u[m-1,j] = Tb # Matriz de ecuaciones tamano = m-2 A = np.zeros(shape=(tamano,tamano), dtype = float) B = np.zeros(tamano, dtype = float) for f in range(0,tamano,1): if (f>0): A[f,f-1]=P A[f,f] = Q if (f<(tamano-1)): A[f,f+1]=R B[f] = -u[f+1,j-1] B[0] = B[0]-P*u[0,j] B[tamano-1] = B[tamano-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,tamano,1): u[f+1,j] = C[f] # siguiente iteración j = j + 1 # SALIDA print('Tabla de resultados') np.set_printoptions(precision=2) print(u)
algunos valores:
Tabla de resultados [[ 60. 60. 60. ..., 60. 60. 60. ] [ 25. 31.01 35.25 ..., 57.06 57.09 57.11] [ 25. 26.03 27.49 ..., 54.22 54.26 54.31] ..., [ 25. 25.44 26.07 ..., 42.22 42.27 42.31] [ 25. 27.57 29.39 ..., 41.07 41.09 41.11] [ 40. 40. 40. ..., 40. 40. 40. ]]
Para realizar la gráfica se aplica lo mismo que en el método explícito
# Gráfica salto = int(n/10) if (salto == 0): salto = 1 for j in range(0,n,salto): vector = u[:,j] plt.plot(xi,vector) plt.plot(xi,vector, '.m') plt.xlabel('x[i]') plt.ylabel('t[j]') plt.title('Solución EDP parabólica') 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 pdf915, Burden 9Ed 12.2 p727, Rodriguez 10.2.2 pdf409 .
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:
http://blog.espol.edu.ec/analisisnumerico/recursos/resumen-python/tiempos-de-ejecucion/