How I can implement the advection model with PCs

The model is

{u_t}^{+} +{\gamma}{u_x}^{+}= \lambda(v)(u^{-}-u^{+})-\beta{u^+}{v}

{u_t}^{-} -{\gamma}{u_x}^{-}= \lambda(v)(u^{+}-u^{-})-\beta{u^-}{v}

v_t= \beta{u v}-\alpha{v}, , \text{with}\;\;\; u=u^{+}+u^{-}.

Regarding the initial conditions we use

u^+(x, 0)=u_{*}^{+}+a_1 sin(10xk_1),\; u^-(x, 0)=u_{*}^{-}+a_2 sin(10xk_1),\nonumber\\

v(x, 0)=v_{*}+a_3 sin(10xk_1).

k_1 := 2 \Pi/L, L=1, and I want to use periodic conditions

Answer 1

[Gautam]

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');

●