# EDP parabólicas d2u/dx2 = K du/dt# método explícito,usando diferencias divididasimport numpy as np
import matplotlib.pyplot as plt
# INGRESO# Valores de frontera
Ta = 1
Tb = 0
#T0 = 25
fx = lambda x: np.cos(3*np.pi/2*x)
# longitud en x
a = 0
b = 1
# Constante K
K = 2
# Tamaño de paso
dx = 0.2
dt = dx/100
# iteraciones en tiempo
n = 10
# PROCEDIMIENTO# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
fi = fx(xi)
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,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = 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
whilenot(j>=ultimot): # igual con lazo forfor i inrange(1,ultimox,1):
u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
j=j+1
# SALIDAprint('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)
# Gráfica
salto = int(n/10)
if (salto == 0):
salto = 1
for j inrange(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()