0

我想在一维反应-扩散-ODE系统的显式有限差分数值方案中实现周期性边界条件。但是，我不是这方面的专家，因此我们将非常感谢任何帮助。代码（我目前有 Neumann 边界条件）如下：

clear all
close all
tic
Fig=figure('Position',[100 100 1200 600]); % Figure Settings [left, bottom, width, height]

% Spatial discretisation
DeltaX=0.05;

% Parameter values
a=0.2;
b=0.15;
h=0.2;
d=10;
epsilon=0.01;
k=0.3;

% Number of grid cells
m=400;
NX=m;

% Timesteps
dt=0.0002; 
Time=0;
EndTime=100000;
PlotStep=10;
PlotTime=PlotStep;

% Initialisation
V = zeros(1,m);
U = zeros(1,m);
S = zeros(1,m);
dV=zeros(1,m);
dU=zeros(1,m);
dS= zeros(1,m);
NetV=zeros(1,m);
NetU=zeros(1,m);

% Boundary conditions
% Neumann
FXV = zeros(1,NX+1);        % bound.con. no flow in/out to X-direction
FXU = zeros(1,NX+1);        % bound.con. no flow in/out to X-direction

% Initial state for V
    for i=1:NX
        if (i <= NX/2)
            V(i)=0.6; 
        else
            V(i)=0;
        end
    end
% Initial state for U
U(:)=61.5;

% figtitle = sgtitle('Time: 0');
figtitle=sgtitle(['A = ' num2str(a) ', B=' num2str(b) ', H=' num2str(h) ', D=' num2str(d) ', eps=' num2str(epsilon) ', k=' num2str(k)]);

index=1;

% Timesteps
while Time<=EndTime
     
    % Reaction
    dV= U.*V.^2 -k*U.*V.^3 -b*V-h*S.*V;
    
    dU= a*(1-U)-U.*V.^2; 
    
    dS= (1/d)*(b*V+h*S.*V-S);
    
    VU=max(0,U); 
    
    % Diffusion
    FXV(2:NX)=(V(2:NX)-V(1:NX-1))/DeltaX; % flow in x-direction : dV/dx;        
    NetV=(epsilon^2)*(FXV(2:NX+1)-FXV(1:NX))/DeltaX; % eps^2 * vector for Laplacian differences
    
    FXU(2:NX)=(VU(2:NX)-VU(1:NX-1))/DeltaX; % flow in x-direction : dU/dx;
    NetU=(FXU(2:NX+1)-FXU(1:NX))/DeltaX; % vector for Laplacian differences 
    
    % Update
    V=V+(NetV+dV)*dt; 
    U=U+(NetU+dU)*dt; 
    S=S+dS*dt;
    
    Time=Time+dt; 

    PlotTime=PlotTime-dt;
  if PlotTime<=0.00001
     
    figure(1)
    % V plot    
    subplot(1,3,1)
    plot(V)
    xlabel('x');
    ylabel('V');
    ylim([0 3]);
    % U plot    
    subplot(1,3,2)
    plot(U)
    xlabel('x');
    ylabel('U');
    ylim([0 1]);
    % S plot 
    subplot(1,3,3)
    plot(S)
    xlabel('x');
    ylabel('S');
    ylim([0 2]);
    
    PlotTime=PlotStep; 
    
  end
end
4

0 回答 0