RBF_aula_extendida_3

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

import numpy as np
from scipy.optimize import curve_fit
from scipy.linalg import solve
from matplotlib import pyplot as plt
# Constantes del problema
D = 5000 # m^2/hr
U = 100 # m/hr
k = 2 # 1/hr
L = 100 # m
cin = 100 # mol/L
# Solución exacta
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))
def c_exact(x):
 return U * cin / ((U - D * lamda1) * lamda2 * np.exp(lamda2 * L)) * (lamda2 * np.exp(lamda2 * L) * np.exp(lamda1 * x) - lamda1 * np.exp(lamda1 * L) * np.exp(lamda2 * x)) 
def irbf(x, c):
 return 1 / np.sqrt(1 + c * x**2)
def solve_rbf_irbf(n, c):
 # Puntos de colocación
 x = np.linspace(0, L, n)
 
 # Matriz de coeficientes
 A = np.zeros((n, n))
 for i in range(n):
 for j in range(n):
 if i == j:
 A[i, j] = irbf(0, c) + irbf(L, c)
 else:
 A[i, j] = irbf(np.abs(x[i] - x[j]), c)
 
 # Vector de términos conocidos
 b = np.zeros(n)
 b[0] = U * cin / irbf(0, c)
 
 # Imponer la condición de frontera en x = L
 A[n-1, :] = 0
 A[n-1, n-1] = 1
 
 # Resolver el sistema de ecuaciones
 c = solve(A, b)
 
 return x, c
# Valores de c a probar
c_vals = np.linspace(0.1, 10, 100)
# Calcular el error L2 para cada valor de c
l2_norms = []
for c in c_vals:
 # Calcular la aproximación RBF
 x, c_rbf = solve_rbf_irbf(101, c)
 
 # Calcular el error L2
 error = c_exact(x) - c_rbf
 l2_norm = np.sqrt(np.sum(error**2) / len(x))
 l2_norms.append(l2_norm)
 
# Trazar la figura
fig, ax = plt.subplots()
ax.plot(c_vals, l2_norms)
ax.set_xlabel('Parámetro de forma c')
ax.set_ylabel('Error L2')
ax.set_title('Norma del error L2 en función del parámetro de forma')
plt.show()