Ejercicio: 3Eva2018TI_T3 EDP Parabólica, temperatura en varilla
Se generan las ecuaciones usando diferencias finitas divididas centradas y hacia adelante
\frac{d^2u}{dx^2} + \frac{Kr}{\rho C} = K\frac{du}{dt}\frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2}+ + \frac{Kr}{\rho C} = K \frac{u[i,j+1]-u[i,j]}{\Delta t}
\frac{\Delta t}{K\Delta x^2} \Big[u[i-1,j]-2u[i,j]+u[i+1,j] \Big] + + \frac{Kr}{\rho C} \frac{\Delta t}{K} = u[i,j+1]-u[i,j]
Se sustituye :
\lambda = \frac{\Delta t}{K\Delta x^2} \gamma = \frac{Kr}{\rho C} \frac{\Delta t}{K} = \frac{r\Delta t}{\rho C}simplificando a:
\lambda \Big[u[i-1,j]-2u[i,j]+u[i+1,j] \Big] + +\gamma = u[i,j+1]-u[i,j]
despejando para u[i,j+1]:
\lambda u[i-1,j]-2\lambda u[i,j]+\lambda u[i+1,j] + +\gamma = u[i,j+1]-u[i,j]u[i,j+1] =\lambda u[i-1,j]+(1-2\lambda) u[i,j]+ +\lambda u[i+1,j] +\gamma
con lo que se tiene una forma explicita de encontrar los valores de la ecuación.
La gráfica se realizó para 20 valores de t con tamaño de paso dt
u[:,t] para t = 0.07500000000000001
[ 0. 0.81 1.23 1.35 1.23 0.81 0. ]algunos valores de u[i,j]
Tabla de resultados
[[ 0. 0. 0. 0. 0. 0. ...]
[ 0.5 0.66 0.74 0.81 0.87 0.92 ...]
[ 0.87 0.99 1.12 1.23 1.33 1.41 ...]
[ 1. 1.11 1.23 1.35 1.47 1.57 ...]
[ 0.87 0.99 1.12 1.23 1.33 1.41 ...]
[ 0.5 0.66 0.74 0.81 0.87 0.92 ...]
[ 0. 0. 0. 0. 0. 0. ...]]Algoritmo en Python
# 3ra Evaluación I Término 2018
# Tema 3. EDP Parabólica, Temperatura en varilla
# método explícito, usando diferencias finitas
import numpy as np
import matplotlib.pyplot as plt
# INGRESO
# Constantes
L = 1.5
K = 1.04
ro = 10.6
C = 0.056
r = 5.0
# Tamaño de paso
dx = 0.25
dt = 0.025
# longitud en x
a = 0
b = L
# iteraciones en tiempo
n = 20
# Valores de frontera
Ta = 0
Tb = 0
ux0 = lambda x: np.sin(np.pi*x/L)
# 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,j]= Ta
u[1:ultimox,j] = ux0(xi[1:ultimox])
u[ultimox,j] = Tb
# factores P,Q,R
lamb = dt/(K*(dx**2))
gama = r*dt/(ro*C)
P = lamb
Q = 1 - 2*lamb
R = lamb
# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot):
u[0,j+1] = Ta
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] + gama
u[m-1,j+1] = Tb
j=j+1
# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)
print('u[:,t] para t =', 3*dt)
print(u[:,3])
# Gráfica
salto = 1 # 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('U[x,tj]')
plt.title('Solución EDP parabólica')
plt.show()
algunos valores de u[i,j]
Tabla de resultados
[[ 0. 0. 0. 0. 0. 0. ...]
[ 0.5 0.66 0.74 0.81 0.87 0.92 ...]
[ 0.87 0.99 1.12 1.23 1.33 1.41 ...]
[ 1. 1.11 1.23 1.35 1.47 1.57 ...]
[ 0.87 0.99 1.12 1.23 1.33 1.41 ...]
[ 0.5 0.66 0.74 0.81 0.87 0.92 ...]
[ 0. 0. 0. 0. 0. 0. ...]]