耦合对流弥散方程差分法求解

想采用差分法求解以下方程

img

采用了以下显式差分方法

img

求解结果不稳定

希望有专家能指点正确的差分格式

你好，这个差分格式可以这么写的

img

边界条件是要注意的，我用的格式可能跟你的不一样。
当然用显式差分格式是比较方便的
下面是我的试验代码：

clc; clear
L = 1;
nx = 101;
x = linspace(0, L, nx);
tf = 1;
nt = 100000;
t = linspace(0, tf, nt);
C = zeros(nt, nx);
S = zeros(nt, nx);
theta_w = 1;
rho_b = 1;
S_max = 1;
k_att  = 1;
k_det = 1;
C0 = 1;
vp = 1;
D = 1;
dx = mean(diff(x));
dt = mean(diff(t));
for i = 1:nt-1
    C(i,1) = C0;
    dSdt = theta_w/rho_b*((S_max - S(i,:))/S_max*k_att.*C(i,:) - rho_b/theta_w*k_det*S(i,:));
    S(i+1,:) = S(i,:) + dSdt*dt;
    dCdt = zeros(size(dSdt));
    dCdt(2:end-1) = D*conv(C(i,:), [1,-2,1], 'valid')/(dx)^2 - ...
        vp * conv(C(i,:), [-1,0,1], 'valid')/2/dx - rho_b/theta_w*dSdt(2:end-1);
    C(i+1, 2:end-1) = C(i,2:end-1) + dCdt(2:end-1)*dt;
    C(i+1,1) = C0;
    C(i+1,nx) = C(i, nx);
end

可见当时间步离散的足够小，这个算法是收敛的，而且很稳定

本人也在准备数模，
差分法就是需要用到循环
然后格式就类似下面这种吧，上次数模也查过这个公式，照着资料敲了部分
C(i,1) = C0;
dSdt = theta_w/rho_b*((S_m - S(i,:))/S_mk_att.C(i,:) - rho_b/theta_wk_detS(i,:));
S(i+1,:) = S(i,:) + dSdtdt;
dCdt = zeros(size(dSdt));
dCdt(2:end-1) = Dconv(C(i,:), [1,-2,1], 'valid')/(dx)^2 - ...
vp * conv(C(i,:), [-1,0,1], 'valid')/2/dx - rho_b/theta_w*dSdt(2:end-1);
C(i+1, 2:end-1) = C(i,2:end-1) + dCdt(2:end-1)*dt;
C(i+1,1) = C0;