matlab - 使用 FD 方法求解一维平流-频散方程得到的负值

Question

我正在尝试解决 1D ADE

到目前为止，这是我的代码：

clc; clear; close all
 
%Input parameters
Ao    = 1;      %Initial value
L     = 0.08;   %Column length [m]
nx    = 40;     %spatial gridpoints
dx    = L/nx;   %Length step size [m]
T     = 20/24;  %End time [days]
nt    = 100;    %temporal gridpoints
dt    = T/nt;   %Time step size [days] 
Vel   = dx/dt;  %Velocity in each cell [m/day]
alpha = 0.002;  %Dispersivity [m]
De    = alpha*Vel;   % Dispersion coeff. [m2/day]
 
%Gridblocks
x    = 0:dx:L;
t    = 0:dt:T;
 
%Initial and boundary conditions
f  = @(x) x;   % initial cond. 
% boundary conditions 
g1 = @(t) Ao;   
g2 = @(t) 0;
 
%Initialization
A      =  zeros(nx+1, nt+1);
A(:,1) =  f(x);     
A(1,:) =  g1(t);    
 
gamma = dt/(dx^2);
beta  = dt/dx;
 
% Implementation of the explicit method
for j= 1:nt-1      % Time Loop
    for i= 2:nx-1  % Space Loop
        A(i,j+1) = (A(i-1,j))*(Vel*beta + De*gamma)...
            + A(i,j)*(1-2*De*gamma-Vel*beta) + A(i+1,j)*(De*gamma);
    end
    
    % Insert boundary conditions for i = 1 and i = N
    A(2,j+1) =  A(1,j)*(Vel*beta + De*gamma) + A(2,j)*(1-2*De*gamma-Vel*beta) + A(3,j)*(De*gamma);
    A(nx,j+1) =  A(nx-1,j)*(Vel*beta + 2*De*gamma) + A(nx,j)*(1-2*De*gamma-Vel*beta)+ (2*De*gamma*dx*g2(t));
end
 
figure
plot(t, A(end,:), 'r*', 'MarkerSize', 2)
title('A Vs time profile (Using FDM)')
xlabel('t'),ylabel('A')

现在，我已经能够使用 MATLAB 的 pdepe 函数（参见图）解决问题，但我试图将结果与有限差分法（在上面的代码中实现）进行比较。然而，我得到了因变量的负值，但我不确定我到底做错了什么。因此，如果有人能在这里帮助我，我将不胜感激。非常感谢期待。

PS：如果有人想看到它，我可以发布我用于 pdepe 的代码。