想采用差分法求解以下方程
采用了以下显式差分方法
求解结果不稳定
希望有专家能指点正确的差分格式
你好,这个差分格式可以这么写的
边界条件是要注意的,我用的格式可能跟你的不一样。
当然用显式差分格式是比较方便的
下面是我的试验代码:
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
可见当时间步离散的足够小,这个算法是收敛的,而且很稳定