Ejercicio: 2Eva_IIT2017_T3 EDP parabólica con diferencias regresivas
\frac{dU}{dt} - \frac{1}{16} \frac{d^2U}{dx^2} = 0
Las diferencias finitas requidas en el enunciado son:
U'(x_i,t_j) = \frac{U(x_{i},t_j)-U(x_{i},t_{j-1})}{\Delta t} + O(\Delta t)
U''(x_i,t_j) = \frac{U(x_{i+1},t_j)-2U(x_{i},t_j)+U(x_{i-1},t_j)}{\Delta x^2} + O(\Delta x^2)
La indicación de regresiva es para la primera derivada, dependiente de tiempo t.
que al reemplazar en la ecuación sin el término de error, se convierte a.
\frac{U(x_{i},t_j)-U(x_{i},t_{j-1})}{\Delta t} - \frac{1}{16}\frac{U(x_{i+1},t_j)-2U(x_{i},t_j)+U(x_{i-1},t_j)}{\Delta x^2} =0
Se reordenan los términos de forma semejante al modelo planteado en el método básico:
\frac{\Delta t}{16\Delta x^2}[U(x_{i+1},t_j)-2U(x_{i},t_j)+U(x_{i-1},t_j)] = U(x_{i},t_j)-U(x_{i},t_{j-1})
Se simplifica haciendo que haciendo
\lambda = \frac{\Delta t}{16\Delta x^2}
Cambiando la nomenclatura con solo los índices para las variables x y t, ordenando en forma ascendente los índices previo a crear el algoritmo.
\lambda[U(i+1,j)-2U(i,j)+U(i-1,j)] = U(i,j)-U(i,j-1)
Se reordena la ecuación como modelo para el sistema de ecuaciones.
\lambda U(i+1,j)+(-2\lambda-1)U(i,j)+ \lambda U(i-1,j) = -U(i,j-1)
P U(i-1,j) + Q U(i,j) + R U(i+1,j) = -U(i,j-1)
Se calculan los valores constantes:
λ = dt/(16*dx2) = 0.05/[16*(1/3)2] = 0.028125
P = λ = 0.028125
Q = (-1-2λ) = (1-2*0.028125) = -1.05625
R = λ = 0.028125
Usando las condiciones del problema:
U(0,t) = U(1,t) = 0, entonces, Ta = 0, Tb = 0
Para los valores de la barra iniciales se debe usar un vector calculado como 2sin(π x) en cada valor de xi espaciados por hx = 1/3, x entre [0,1]
xi = [0,1/3, 2/3, 1]
U[xi,0] = [2sin (0*π), 2sin(π/3), 2sin(2π/3), 2sin(π)]
U[xi,0] = [0, 2sin(π/3), 2sin(2π/3), 0]
U[xi,0] = [0, 1.732050, 1.732050, 0]
Con lo que se puede plantear las ecuaciones:
j=1: i=1
0.028125 U(0,1) + (-1.05625) U(1,1) + 0.028125 U(2,1) = -U(1,0)
j=1: i=2
0.028125 U(1,1) + (-1.05625) U(2,1) + 0.028125 U(3,1) = -U(2,0)
y reemplazando los valores de la malla conocidos:
0.028125 (0) - 1.05625 U(1,1) + 0.028125 U(2,1) = -1.732050
0.028125 U(1,1) - 1.05625 U(2,1) + 0.028125 (0) = -1.732050
hay que resolver el sistema de ecuaciones:
-1.05625 U(1,1) + 0.028125 U(2,1) = -1.732050
0.028125 U(1,1) - 1.05625 U(2,1) = -1.732050
A = [[-1.05625 , 0.028125],
[ 0.028125, -1.05625 ]]
B = [-1.732050,-1.732050]
que resuelta con un método numérico:
[ 1.68, 1.68]
Por lo que la solución para una gráfica, con los índices de (fila,columna) como (t,x):
U = [[0, 1.732050, 1.732050, 0],
[0, 1.680000, 1,680000, 0]]
El error del procedimiento, tal como fué planteado es del orden de O(Δt) y O(Δx2), o error de truncamiento E = O(Δx2) + O(Δt). Δt debe ser menor que Δx en aproximadamente un orden de magnitud
Usando algoritmo en python.
Usando lo resuelto en clase y laboratorio, se comprueba la solución con el algoritmo, con hx y ht mas pequeños y más iteraciones:
# 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 = 0
Tb = 0
# longitud en x
a = 0
b = 1
# Constante K
K = 16
# Tamaño de paso
dx = 0.1
dt = 0.01
# temperatura en barra
tempbarra = lambda x: 2*np.sin(np.pi*x)
# iteraciones
n = 100
# PROCEDIMIENTO
# Valores de x
xi = np.arange(a,b+dx,dx)
m = len(xi)
# Resultados en tabla de u
u = np.zeros(shape=(m,n), dtype=float)
# valores iniciales de u
j=0
u[0,j] = Ta
for i in range(1,m-1,1):
u[i,j] = tempbarra(xi[i])
u[m-1,j] = Tb
# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = -1 -2*lamb
R = lamb
vector = np.array([P,Q,R])
tvector = len(vector)
# 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):
for c in range(0,tvector,1):
c1 = f+c-1
if(c1>=0 and c1<tamano):
A[f,c1] = vector[c]
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]
j=j+1 # siguiente iteración
# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)
# 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()
