Algoritmos101

s3Eva2018TI_T3 EDP Parabólica, temperatura en varilla

por

Edison Del Rosario

Ejercicio: 3Eva_IT2018_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 ecuacion.

La gráfica se realizó para 20 valores de t con tamaño de paso dt


El resultado pedido en el enunciado de distribución de temperatura para t=3k

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.    ...]]

Ejemplos