Publicado en por Edison Del Rosario

s2Eva_IIT2011_T3 EDP Parabólica, explícito

Ejercicio: 2Eva_IIT2011_T3 EDP Parabólica, explícito

La ecuación a resolver es:

δuδt−δ2uδx2=2 \frac{\delta u}{\delta t} - \frac{\delta^2 u}{\delta x^2} =2

Como el método requerido es explícito se usan las diferencias divididas:

d2udx2=ui+1,j−2ui,j+ui−1,j(Δx)2 \frac{d^2 u}{dx^2} = \frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2} dudt=ui,j+1−ui,jΔt \frac{du}{dt} = \frac{u_{i,j+1} - u_{i,j} }{\Delta t}

Reordenando la ecuación al modelo realizado en clase:

δ2uδx2=δuδt−2\frac{\delta^2 u}{\delta x^2} = \frac{\delta u}{\delta t} - 2 ui+1,j−2ui,j+ui−1,j(Δx)2=ui,j+1−ui,jΔt−2\frac{u_{i+1,j} - 2 u_{i,j} + u_{i-1,j}}{(\Delta x)^2} = \frac{u_{i,j+1} - u_{i,j} }{\Delta t} - 2

multiplicando cada lado por Δt

Δt(Δx)2[ui+1,j−2ui,j+ui−1,j]=ui,j+1−ui,j−2Δt \frac{\Delta t}{(\Delta x)^2} \Big[u_{i+1,j} - 2 u_{i,j} + u_{i-1,j} \Big]= u_{i,j+1} - u_{i,j} - 2\Delta t

Se establece el valor de

λ=Δt(Δx)2 \lambda = \frac{\Delta t}{(\Delta x)^2} λui+1,j−2λui,j+λui−1,j=ui,j+1−ui,j−2Δt\lambda u_{i+1,j} - 2 \lambda u_{i,j} + \lambda u_{i-1,j} = u_{i,j+1} - u_{i,j} - 2\Delta t

obteniendo la ecuación general, ordenada por índices de izquierda a derecha:

λui−1,j+(1−2λ)ui,j+λui+1,j+2Δt=ui,j+1\lambda u_{i-1,j} +(1- 2 \lambda)u_{i,j} + \lambda u_{i+1,j} + 2\Delta t = u_{i,j+1}

con valores de:

λ = Δt/(Δx)2 = (0.04)/(0.25)2 = 0.64
P = λ = 0.64
Q = (1-2λ) = -0.28
R = λ = 0.64
0.64ui−1,j−0.28ui,j+0.64ui+1,j+0.08=ui,j+10.64 u_{i-1,j} - 0.28 u_{i,j} + 0.64 u_{i+1,j} + 0.08 = u_{i,j+1}

Se realizan 3 iteraciones:

i= 1, j=0

u1,1=0.64u0,0−0.28u1,0+0.64u2,0+0.08u_{1,1} = 0.64 u_{0,0} -0.28u_{1,0}+0.64 u_{2,0}+0.08 u1,1=0.64[sin⁡(π∗0)+0∗(1−0)]− u_{1,1} = 0.64 [\sin(\pi*0)+0*(1-0)]- 0.28[sin⁡(π∗0.25)+0.25∗(1−0.25)] 0.28[\sin(\pi*0.25)+0.25*(1-0.25)] +0.64[sin⁡(π∗0.5)+0.5∗(1−0.5)]+0.08+0.64[\sin(\pi*0.5)+ 0.5*(1-0.5)]+0.08

u[1,1] =0.89

i = 2, j=0

0.64u1,0−0.28u2,0+0.64u3,0+0.08=u2,10.64 u_{1,0} - 0.28 u_{2,0} + 0.64 u_{3,0} + 0.08 = u_{2,1}

u[1,0] = 1.25

i = 3, j=0

0.64u2,0−0.28u3,0+0.64u4,0+0.08=u3,10.64 u_{2,0} - 0.28 u_{3,0} + 0.64 u_{4,0} + 0.08 = u_{3,1}

u[3,1] = 0.89

Algoritmo en Python

con los valores y ecuación del problema

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito, usando diferencias finitas
# 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 = 0
Tb = 0
T0 = lambda x: np.sin(np.pi*x)+x*(1-x)
# longitud en x
a = 0
b = 1
# Constante K
K = 1
# Tamaño de paso
dx = 0.1
dt = dx/10/2
# iteraciones en tiempo
n = 50

# 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] = T0(xi[1:ultimox])
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 = 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]+2*dt
    u[m-1,j+1] = Tb
    j=j+1

# 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, '.r')
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()