Hello lulu,
Here's a basic framework for numerically solving a coupled system of advection-reaction PDEs
clear; clc;
% Parameters
L = 1;
xmin = 0; xmax = L;
N = 100; % Number of spatial points
dx = (xmax - xmin)/N;
x = linspace(xmin, xmax-dx, N); % N points, periodic
dt = 0.0005; % Time step
tmax = 0.1; % Final time
nt = round(tmax/dt);
gamma = 0.5; % Advection speed
beta = 1.0; % Reaction parameter
alpha = 0.5; % Decay rate
u_star_p = 1.0; % Steady-state values
u_star_m = 1.0;
v_star = 1.0;
a1 = 0.1; a2 = 0.1; a3 = 0.1;
k1 = 2*pi/L;
lambda = @(v) 1.0; % Example: constant, or try @(v) 1+0.5*v
% Initial conditions
u_p = u_star_p + a1*sin(10*x*k1);
u_m = u_star_m + a2*sin(10*x*k1);
v = v_star + a3*sin(10*x*k1);
% Preallocate for speed
u_p_new = zeros(1,N);
u_m_new = zeros(1,N);
v_new = zeros(1,N);
for n = 1:nt
% Periodic boundary indices
ip = [2:N, 1]; % i+1 with wrap
im = [N, 1:N-1]; % i-1 with wrap
% Advection (upwind)
% For u^+ : upwind left (since +gamma)
u_p_x = (u_p - u_p(im))/dx;
% For u^- : upwind right (since -gamma)
u_m_x = (u_m(ip) - u_m)/dx;
% Coupling/reaction terms
lam = lambda(v);
u = u_p + u_m;
% Update equations (explicit Euler)
u_p_new = u_p + dt * ( ...
- gamma * u_p_x ...
+ lam.*(u_m - u_p) ...
- beta * u_p .* v );
u_m_new = u_m + dt * ( ...
+ gamma * u_m_x ...
+ lam.*(u_p - u_m) ...
- beta * u_m .* v );
v_new = v + dt * ( ...
+ beta * u .* v ...
- alpha * v );
% Advance
u_p = u_p_new;
u_m = u_m_new;
v = v_new;
end
% Plot results
figure;
plot(x, u_p, 'r', x, u_m, 'b', x, v, 'k', 'LineWidth', 2);
legend('u^+','u^-','v');
xlabel('x'); ylabel('Value');
title('Advection-Reaction System with Periodic BC');
