Referencia: Chapra 30.2 p888 pdf912, Burden 9Ed 12.2 p725, Rodriguez 10.2 p406
Siguiendo el tema anterior en 9.1, se tiene que:
\frac{\partial ^2 u}{\partial x ^2} = K\frac{\partial u}{\partial t}Se usan las diferencias divididas, donde se requieren dos valores para la derivada de orden uno y tres valores para la derivada de orden dos:
\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+1} - u_{i,j} }{\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. La línea de referencia es el tiempo t0
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+1}-u_{i,j}}{\Delta t}De la gráfica se destaca que en la fórmula solo hay un valor desconocido, destacado por el punto en amarillo dentro del triángulo. En la ecuación se representa por U[i,j+1].
Despejando la ecuación, se agrupan términos constantes:
λ = \frac{\Delta t}{K (\Delta x)^2}quedando la ecuación, con los términos ordenados de izquierda a derecha como en la gráfica:
u_{i,j+1} = \lambda u_{i-1,j} +(1-2\lambda)u_{i,j}+\lambda u_{i+1,j}Al resolver se encuentra que que cada valor en un punto amarillo se calcula como una suma ponderada de valores conocidos, por lo que el desarrollo se conoce como el método explícito.
La ponderación está determinada por los términos P, Q, y R.
\lambda = \frac{\Delta t}{K(\Delta x)^2} P = \lambda Q = 1-2\lambda R = \lambda u_{i ,j+1} = Pu_{i-1,j} + Qu_{i,j} + Ru_{i+1,j}Fórmulas que se desarrollan usando un algoritmo y considerando que al disminuir los valores de Δx y Δt la cantidad de operaciones aumenta.
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 .
\lambda \leq \frac{1}{2}Cuando λ ≤ 1/2 se tiene como resultado una solución donde los errores no crecen, sino que oscilan.
Haciendo λ ≤ 1/4 asegura que la solución no oscilará.
También se sabe que con λ= 1/6 se tiende a minimizar los errores por truncamiento (véase Carnahan y cols., 1969).
Para resolver la parte numérica suponga:
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
Algoritmo en Python
Se presenta la propuesta en algoritmo para el método explícito.
# EDP parabólicas d2u/dx2 = K du/dt # método explícito,usando diferencias divididas # Referencia: Chapra 30.2 p.888 pdf.912 # Rodriguez 10.2 p.406 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 = dx/10 # iteraciones en tiempo n = 200 # PROCEDIMIENTO # iteraciones en longitud xi = np.arange(a,b+dx,dx) m = len(xi) ultimox = m-1 # Resultados en tabla u[x,t] u = np.zeros(shape=(m,n), dtype=float) # valores iniciales de u[:,j] j=0 ultimot = n-1 u[0,:]= Ta u[1:ultimox,j] = T0 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 j = 0 while not(j>=ultimot): # igual con lazo for for i in range(1,ultimox,1): u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j] j=j+1 # SALIDA print('Tabla de resultados') np.set_printoptions(precision=2) print(u)
Se muestra una tabla resumida de resultados a forma de ejemplo.
Tabla de resultados [[ 60. 60. 60. ..., 60. 60. 60. ] [ 25. 33.75 38.12 ..., 57.93 57.93 57.93] [ 25. 25. 27.19 ..., 55.86 55.86 55.87] ..., [ 25. 25. 25.94 ..., 43.86 43.86 43.87] [ 25. 28.75 30.62 ..., 41.93 41.93 41.93] [ 40. 40. 40. ..., 40. 40. 40. ]] >>>
Si la cantidad de puntos aumenta al disminuir Δx y Δt, es mejor disminuir la cantidad de curvas a usar en el gráfico para evitar superponerlas. Se usa el parámetro ‘salto’ para indicar cada cuántas curvas calculadas se incorporan en la gráfica.
# 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, '.r') plt.xlabel('x[i]') plt.ylabel('t[j]') plt.title('Solución EDP parabólica') plt.show()
Note que en la gráfica se toma como coordenadas x vs t, mientras que para la solución de la malla en la matriz las se usan filas y columnas. En la matriz el primer índice es fila y el segúndo índice es columna.
Tarea o Proyecto: Realizar la animación de los cambios de temperatura en el tiempo.
