Ejercicio: 2Eva2022PAOII_T3 EDP Parabólica con coseno 3/4π
\frac{\partial^2 u}{\partial x^2} = b \frac{\partial u}{\partial t}Literal a. gráfica de malla
literal b. Planteamiento
\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = b\frac{u_{i,j+1}-u_{i,j}}{\Delta t}literal c. Modelo discreto
Agrupando constantes(λ)
\frac{\Delta t}{b} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = \frac{\Delta t}{b}b\frac{u_{i,j+1}-u_{i,j}}{\Delta t} λ = \frac{\Delta t}{b(\Delta x)^2}Analiza λ para convergencia del método
λ = \frac{0.002}{2(0.2)^2} =0.025como λ<0.5 el método converge.
\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] = u[i,j+1]-u[i,j]literal d. Usando el método explícito
reagrupando nodos
u[i,j+1] =\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] + u[i,j] u[i,j+1] =\lambda u[i+1,j]+(1-2\lambda)u[i,j]+\lambda u[i-1,j]iteración i=1, j=0
u[1,1] =\lambda u[0,0]+(1-2\lambda)u[1,0]+\lambda u[2,0] u[1,1] =0.025(1) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.2\Big)+0.025 \cos \Big( \frac{3π}{2}0.4\Big)iteración i=2, j=0
u[2,1] =\lambda u[1,0]+(1-2\lambda)u[2,0]+\lambda u[3,0] u[2,1] =0.025\cos \Big( \frac{3π}{2}0.2\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.4\Big) +0.025 \cos \Big( \frac{3π}{2}0.6\Big)iteración i=3, j=0
u[3,1] =\lambda u[2,0]+(1-2\lambda)u[3,0]+\lambda u[4,0] u[3,1] =0.025\cos \Big( \frac{3π}{2}0.4\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.6\Big) +0.025 \cos \Big( \frac{3π}{2}0.8\Big)iteración i=4, j=0
u[4,1] =\lambda u[3,0]+(1-2\lambda)u[4,0]+\lambda u[5,0] u[4,1] =0.025\cos \Big( \frac{3π}{2}0.6\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.8\Big) +0.025 (0)continuar con las iteraciones en el algoritmo
literal e. Algoritmo en Python
Se observa que el método es convergente, dado que λ<0.5
Método explícito EDP Parabólica
lambda: 0.02
x: [0. 0.2 0.4 0.6 0.8 1. ]
t: [0. 0. 0. 0.01 0.01] ... 0.02
Tabla de resultados en malla EDP Parabólica
j, U[:,: 5 ], primeras iteraciones de 50
5 [ 1. 0.53 -0.28 -0.86 -0.73 0. ]
4 [ 1. 0.54 -0.28 -0.88 -0.74 0. ]
3 [ 1. 0.55 -0.29 -0.89 -0.76 0. ]
2 [ 1. 0.56 -0.3 -0.91 -0.78 0. ]
1 [ 1. 0.58 -0.3 -0.93 -0.79 0. ]
0 [ 1. 0.59 -0.31 -0.95 -0.81 0. ]gráfica de temperatura inicial de barra
Resultado luego de al menos 10 iteraciones
Gráfica en 3D
Algoritmo con Python
# 2Eva2022PAOII_T3 EDP Parabólica con coseno 3/4π
# EDP parabólicas d2u/dx2 = K du/dt
# método explícito,usando diferencias divididas
import numpy as np
# INGRESO
# Valores de frontera
Ta = 1 # izquierda de barra
Tb = 0 # derecha de barra
#Tc = 25 # estado inicial de barra en x
fxc = lambda x: np.cos(3*np.pi/2*x) # f(x) en tiempo inicial
# dimensiones de la barra
a = 0 # longitud en x
b = 1
t0 = 0 # tiempo inicial, aumenta con dt en n iteraciones
K = 2 # Constante K
dx = 0.2 # muestreo en x, tamaño de paso
dt = dx/100
n = 50 # iteraciones en tiempo
verdigitos = 2 # decimales a mostrar en tabla de resultados
# coeficientes de U[x,t]. factores P,Q,R,
lamb = dt/(K*dx**2)
P = lamb # izquierda P*U[i-1,j]
Q = 1 - 2*lamb # centro Q*U[i,j]
R = lamb # derecha R*U[i+1,j]
# PROCEDIMIENTO
# iteraciones en x, longitud
xi = np.arange(a,b+dx/2,dx)
m = len(xi)
ultimox = m-1
ultimot = n-1
# u[xi,tj], tabla de resultados
# se incluye el estado a t=0: n+1
u = np.zeros(shape=(m,n+1), dtype=float)
# u[i,j], valores iniciales
u[0,:] = Ta # Izquierda
u[ultimox,:] = Tb # derecha
# estado inicial de barra en x, Tc
fic = fxc(xi) # f(x) en tiempo inicial
u[1:ultimox,0] = fic[1:ultimox]
# Calcula U para cada tiempo + dt
tj = np.arange(t0,(n+1)*dt,dt)
for j in range(0,n,1):
for i in range(1,ultimox,1):
# ecuacion discreta, entre [1,ultimox]
u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
# SALIDA
j_mostrar = 5
np.set_printoptions(precision=verdigitos)
print('Método explícito EDP Parabólica')
print('lambda: ',np.around(lamb,verdigitos))
print('x:',xi)
print('t:',tj[0:j_mostrar],'...',tj[-1])
print('Tabla de resultados en malla EDP Parabólica')
print('j, U[:,:',j_mostrar,'], primeras iteraciones de ',n)
for j in range(j_mostrar,-1,-1):
print(j,u[:,j])
Gráficas en Python
# GRAFICA ------------
import matplotlib.pyplot as plt
tramos = 10
salto = int(n/tramos) # evita muchas líneas
if (salto == 0):
salto = 1
for j in range(0,n,salto):
vector = u[:,j]
plt.plot(xi,vector)
plt.plot(xi,vector, '.',color='red')
plt.xlabel('x[i]')
plt.ylabel('u[j]')
plt.title('Solución EDP parabólica')
plt.show()
# GRAFICA en 3D ------
tj = np.arange(0,n*dt,dt)
tk = np.zeros(tramos,dtype=float)
# Extrae parte de la matriz U,acorde a los tramos
U = np.zeros(shape=(m,tramos),dtype=float)
for k in range(0,tramos,1):
U[:,k] = u[:,k*salto]
tk[k] = tj[k*salto]
# Malla para cada eje X,Y
Xi, Yi = np.meshgrid(xi,tk)
U = np.transpose(U)
fig_3D = plt.figure()
graf_3D = fig_3D.add_subplot(111, projection='3d')
graf_3D.plot_wireframe(Xi,Yi,U, color ='blue')
graf_3D.plot(xi[0],tk[1],U[1,0],'o',color ='orange')
graf_3D.plot(xi[1],tk[1],U[1,1],'o',color ='green')
graf_3D.plot(xi[2],tk[1],U[1,2],'o',color ='green')
graf_3D.plot(xi[1],tk[2],U[2,1],'o',color ='salmon',
label='U[i,j+1]')
graf_3D.set_title('EDP Parabólica')
graf_3D.set_xlabel('x')
graf_3D.set_ylabel('t')
graf_3D.set_zlabel('U')
graf_3D.legend()
graf_3D.view_init(35, -45)
plt.show()