我的这段代码的时间不多了，所以任何帮助都将不胜感激。我目前正在 matlab 中编写 2D 热/扩散方程，但在添加源项时遇到问题。我的等式如下所示：

$$\frac{\partial u}{\partial t} = \alpha_x \frac{\partial^2 u}{\partial x^2} + \alpha_y \frac{\partial^2 u}{\partial y^2} + Au \left( 1 - \frac{u}{K} \right)$$

我到处都在使用带有通量 = 0 的 Neumann bc。u 是人口数，A 是内在增长率，K 是承载能力。我正在使用的代码是我在此处找到的代码，我已对其进行编辑以使用隐式/BTCS 方法。

太感谢了。

%%
%Specifying parameters
nx=40;                           %Number of steps in space(x)
ny=50;                           %Number of steps in space(y)       
nt=50;                           %Number of time steps 
dt=0.01;                         %Width of each time step
dx=2/(nx-1);                     %Width of space step(x)
dy=2/(ny-1);                     %Width of space step(y)
x=0:dx:2;                        %Range of x(0,2) and specifying the grid points
y=0:dy:2;                        %Range of y(0,2) and specifying the grid points
u=zeros(nx,ny);                  %Preallocating u
un=zeros(nx,ny);                 %Preallocating un
visx = 2;                        %Diffusion coefficient/viscocity
visy = 2;                        %Diffusion coefficient/viscocity
% RepRateA = 0.3;                %Intrinsic Growth Rate of Population
% CarCapK = 1.7;                 %Carrying Capacity of Environment 

UnW=0;                           %x=0 Neumann B.C (du/dn=UnW)
UnE=0;                           %x=L Neumann B.C (du/dn=UnE)
UnS=0;                           %y=0 Neumann B.C (du/dn=UnS)
UnN=0;                           %y=L Neumann B.C (du/dn=UnN)

%%
%Initial Conditions
for i=1:nx
    for j=1:ny
        if ((0.5<=y(j))&&(y(j)<=1.5)&&(0.5<=x(i))&&(x(i)<=1.5))
            u(i,j)=2;
        else
            u(i,j)=0;
        end
    end
end

%%
%B.C vector
bc=zeros(nx-2,ny-2);

bc(1,:)=-UnW/dx;       %Neumann B.Cs
bc(nx-2,:)=UnE/dx;     %Neumann B.Cs
bc(:,1)=-UnS/dy;       %Neumann B.Cs
bc(:,nx-2)=UnN/dy;     %Neumann B.Cs

%B.Cs at the corners:
bc=visx*dt*bc+visy*dt*bc;

%Calculating the coefficient matrix for the implicit scheme
Ex=sparse(2:nx-2,1:nx-3,1,nx-2,nx-2);
% Ax=Ex+Ex'-2*speye(nx-2);
Ax(1,1)=-1; 
Ax(nx-2,nx-2)=-1;                 %Neumann B.Cs
Ey=sparse(2:ny-2,1:ny-3,1,ny-2,ny-2);
% Ay=Ey+Ey'-2*speye(ny-2);
Ay(1,1)=-1; 
Ay(ny-2,ny-2)=-1;                 %Neumann B.Cs
A=kron(Ay/dy^2,speye(nx-2))+kron(speye(ny-2),Ax/dx^2);
D=speye((nx-2)*(ny-2))-visx*dt*A-visy*dt*A;
%%
%Calculating the field variable for each time step
i=2:nx-1;
j=2:ny-1;

for it=0:nt
    un=u;
    h=surf(x,y,u','EdgeColor','none');       %plotting the field variable
    shading interp
    axis ([0 2 0 2 0 2])
    title({['2-D Diffusion with alphax = ',num2str(visx)];['2-D Diffusion with alphay = ',num2str(visy)];['time (\itt) = ',num2str(it*dt)]})
    xlabel('Spatial co-ordinate (x) \rightarrow')
    ylabel('{\leftarrow} Spatial co-ordinate (y)')
    zlabel('Transport property profile (u) \rightarrow')
    drawnow; 
    refreshdata(h)

    %Implicit method:
    U=un;
    U(1,:)=[];
    U(end,:)=[];
    U(:,1)=[];
    U(:,end)=[];
    U=reshape(U+bc,[],1);
    U=D\U;
    U=reshape(U,nx-2,ny-2);
    u(2:nx-1,2:ny-1)=U;
%Neumann:
    u(1,:)=RepRateA*u(2,:)*(1-u(2,:).'/CarCapK)-UnW*dx;
    u(nx,:)=RepRateA*u(nx-1,:)*(1-u(nx-1,:).'/CarCapK)+UnE*dx;
    u(:,1)=RepRateA*u(:,2).'*(1-u(:,2)/CarCapK)-UnS*dy;
    u(:,ny)=RepRateA*u(:,ny-1).'*(1-u(:,ny-1)/CarCapK)+UnN*dy;
end