RBF_aula_extendida_2

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
# Definir parámetros
D = 5000 # m^2/hr
U = 100 # m/hr
k = 2 # 1/hr
L = 100 # m
cin = 100 # mol/L
lamda1 = U/(2*D)*(1 - np.sqrt(1 + 4*k*D/U**2))
lamda2 = U/(2*D)*(1 + np.sqrt(1 + 4*k*D/U**2))
# Definir la solución exacta
def exact_sol(x):
 return (U*cin)/((U - D*lamda1)*lamda2*np.exp(lamda2*L) - (U - D*lamda2)*lamda1*np.exp(lamda1*L))*(lamda2*np.exp(lamda2*L)*np.exp(lamda1*x) - lamda1*np.exp(lamda1*L)*np.exp(lamda2*x))
# Definir la ecuación diferencial y las condiciones de frontera
def ode(x, y):
 return np.vstack((y[1], U*y[1]/D - k*y[0]))
def bc(ya, yb):
 return np.array([ya[0] - (U/D - lamda1)*ya[1], yb[1]])
# Definir la función de base radial gaussiana
def phi(r, eps):
 return np.exp(-(eps*r)**2)
# Definir la función de aproximación
def approx_func(x, c, eps):
 f = np.zeros_like(x)
 for i in range(len(c)):
 f += c[i]*phi(np.abs(x - x[i]), eps)
 return f
# Definir la función para calcular los coeficientes de la solución aproximada
def solve_coefficients(x, y, eps):
 N = len(x)
 A = np.zeros((N, N))
 for i in range(N):
 for j in range(N):
 A[i,j] = phi(np.abs(x[i] - x[j]), eps)
 b = y
 c = np.linalg.solve(A, b)
 return c
# Definir el número de nodos
n = 200
# Definir los valores de eps para la RBF gaussiana
eps_list = [0.01, 0.05, 0.1, 0.5, 1.0]
# Inicializar la matriz de errores locales
error_matrix = np.zeros((n, len(eps_list)))
# Resolver la ecuación diferencial y calcular los errores locales para cada valor de eps
x = np.linspace(0, L, n)
y = np.zeros((2, n))
y[0, 0] = cin
y = solve_bvp(ode, bc, x, y)
y_exact = exact_sol(x)
for j in range(len(eps_list)):
 eps = eps_list[j]
 c = solve_coefficients(x, y.y[0], eps)
 y_approx = approx_func(x, c, eps)
 error_matrix[:,j] = np.abs(y_approx - y_exact)
# Graficar los errores locales para cada valor de eps
 plt.plot(x, error_matrix[:,j], label=f"eps = {eps_list[j]}")
 plt.legend()
 plt.title("Error local de la solución numérica para diferentes valores de eps")
 plt.xlabel("x")
 plt.ylabel("Error")
 plt.show()