RBF_aula_extendida

Esta es una vista previa del archivo. Inicie sesión para ver el archivo original

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
# Define las constantes del problema
D = 5000.0 # m^2/hr
U = 100.0 # m/hr
k = 2.0 # hr^-1
L = 100.0 # m
cin = 100.0 # mol/L
# Define la solución analítica
lam1_lam2_ratio = U / (2*D) * (1 + np.sqrt(1 + 4*k*D/U**2))
lambda1 = -U / (2*D) + lam1_lam2_ratio
lambda2 = -U / (2*D) - lam1_lam2_ratio
c_analytical = cin * U / (U - D*lambda1) / lambda2 * np.exp(lambda2*L) \
 - (U - D*lambda2) / (U - D*lambda1) * lambda1 / lambda2 * np.exp(lambda1*L) \
 * (lambda2*np.exp(lambda2*L)*np.exp(lambda1*np.arange(0, L+0.1, 0.1)) 
 - lambda1*np.exp(lambda1*L)*np.exp(-lambda2*np.arange(0, L+0.1, 0.1)))
 
# Define la ecuación diferencial
def ode(x, c):
 return np.vstack((c[1], U/D*c[1] + k*c[0]))
# Define las condiciones de frontera
def bc(c0, cL):
 return np.array([c0[1] - U*c0[0]/D + U*cin/D,
 cL[0]])
# Definir la malla espacial
x = np.linspace(0, L, 101)
# Definir la solución inicial
c_guess = np.zeros((2, x.size))
# Resolver la ecuación diferencial con solve_bvp
sol = solve_bvp(ode, bc, x, c_guess, tol=1e-6)
# Graficar la solución numérica y analítica
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(sol.x, sol.y[0], 'b', label='Solución Numérica')
ax.plot(np.arange(0, L+0.1, 0.1), c_analytical, 'r--', label='Solución Analítica')
ax.set_xlabel('Posición (m)')
ax.set_ylabel('Concentración (mol/L)')
ax.set_title('Solución de la Ecuación Diferencial con Funciones de Base Radial')
ax.legend()
plt.show()