import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
def solve_rcd_equation(epsilon, mu, a, b, f, x_min, x_max, N):
# Discretization
dx = (x_max - x_min) / (N - 1)
x = np.linspace(x_min, x_max, N)
# Constructing the differentiation matrix
diagonals = [[epsilon/dx**2], [-(epsilon/dx**2 + mu*a(x)/2*dx)], [epsilon/dx**2 + mu*a(x)/2*dx - b(x)]]
D = diags(diagonals, [-1, 0, 1], shape=(N, N)).toarray()
# Boundary conditions
D[0, 0] = 1
D[0, 1] = 0
D[-1, -1] = 1
D[-1, -2] = 0
# Right-hand side
rhs = f(x)
rhs[0] = 0
rhs[-1] = 0
# Solve the system
y = spsolve(D, rhs)
return x, y
# Example functions for a, b, and f
def a(x):
return 1
def b(x):
return 1
def f(x):
return np.sin(np.pi * x)
# Parameters
epsilon = 0.1
mu = 0.5
x_min = 0
x_max = 1
N = 100
# Solve the equation
x, y = solve_rcd_equation(epsilon, mu, a, b, f, x_min, x_max, N)
# Plot the solution
import matplotlib.pyplot as plt
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of Reaction-Convection-Diffusion Equation')
plt.grid(True)
plt.show()
