A reaction-diffusion equation with delay

The aim of this tutorial is to give a numerical method for solving a partial differential equation with a constant delay.

Introduction

We consider the following one-dimentional reaction-diffusion equation with logistic production and delayed term,

this equation was suggested in [1] as a model of viral infection spreading in tissues. For the existence of solution and Global stability of the homogeneous in space equilibrium we refer the reader to [2]. Here, $u(x,t)$ is the concentration of virus with respect to the space variable $x$ and time $t$. The parameters $D$ and $r$ are respectively the diffusion coefficient and replication rate constant. The function $f(u_{\tau})$ describes the concentration of immune cells as a function of the virus concentration at time $t-\tau, u_{\tau}(x, t) = u(x, t - \tau)$.

Numerical method

We rewrite the one-dimentional reaction-diffusion equation with a constant time delay described above with Neumann boundary condition:

where the delay $\tau$ is a positive constant. we use an implicit finite difference approximation for the diffusion term and classical approach of the resolution of delay equations. let N and M denote the number of space steps and time steps with the notations $\Delta x= L/N$ and $\Delta t= T/M$, respectively.

The discretization of space $x_i$ and time $t_n$ are given by

Let $u^n_i$ be the approximation of the function $u$ at $(x,t)=(x_i,t_n)$, moreover,

We write the scheme of the equation (\ref{equa1}) at the point $u^n_i$,

Where, $u^{n-k}_i \approx u(x_i,t_n-\tau)$ is the delay variable and $k$ is determined by the equality $\tau= k \Delta t$. For simplicity of notation, we denote by $F^{n-k}_i := f(u^{n-k}_i)$ and $G^{n,k}_i:=r u^n_i( 1-u^n_i - F^{n-k}_i)$, for $i=0,..,N-1, \quad n=1,..,k$ and $n=0,..,M-1$.

The simplest form of this expression is given as follows

where $\lambda= D \Delta t /(\Delta x)^2$.

Now, we can write the following semi-linear system with Neumann boundary condition

First, we fix the parametrs as follows $T=10$, $\tau=1$, $L=1$, $D=1$, $r=1$,

clear ;
close ;
clc,
tic
T=10;
M=100;
L=1;
epsilon=10^(-6);
D=0.01;
r=1;
to=2;
N=200;
dx=L/(N);
dt=T/M;
x= (0:N)*dx ;
lambda=dt*D/(dx^2);
t=(0:M)*dt;

The tridiagonal matrix with Neumann condition

C=eye(N+1);
D=(2*lambda+1)*C-lambda*diag(ones(1,N),1)-lambda*diag(ones(1,N),-1);
B=[[-lambda;zeros(N,1)] [zeros(N+1,N-1)] [zeros(N,1);-lambda]];
A=D+B;

The initial function $u_0$

U=zeros(N+1,1);

for i=1:N+1
    U(i,1)=u0(x(i),0);
end
figure (1)
plot(x,U)
xlabel('space')
title('Initial function u_0(x)')

The delay function $f(u_{\tau})$

f=@( t ) 2*t ;
Q=zeros(N+1,1);
for i=1:N+1
 Q(i,1)=f(x(i))   ;
end

figure (2)
plot(x,Q)
xlabel('space')
title('The function f')

Ut=zeros(N+1,1);
Fretard=zeros(N+1,1);
Uf=zeros(N+1,M+1);
Uf=U;

We introduce the following test to get $u(x, t-\tau)$,

fig =        figure;
ax  = axes('Parent',fig,'XLim',[0 L],'YLim',[0 1.1]);
ax.XLabel.String = 'space';

%%gif('pdedelay4.gif','frame',fig,'DelayTime',1)
for n=1:M
    if (((t(n)-to)<=0)||(abs(t(n)-to)<epsilon))
        for i=1:N+1
            Ut(i,1)=u0(x(i),t(n)-to);
            Fretard(i,1)=f(Ut(i,1));
        end
    else
        s=floor((t(n)-to)/dt)+1;
        for i=1:N+1
            Ut(i,1)=Uf(i,s);
            Fretard(i,1)=f(Ut(i,1));
        end
    end
    Uold=U;
    Ur=dt*r*Uold.*(1-Uold-Fretard);
    U=A\(Uold+Ur);
    %%
    %% We plot the evolution of solution $u(x,t)$ at
    %% $t= \Delta t, \tau, 2 \tau, 3\tau, ...$,
    %% when $\tau=0$, we plot the solution for $t=5, 10, 15, ...$
    %%
    if (((t(n+1)==dt)|| (mod(t(n+1),to)==0)) || ((to==0)&&(((mod(t(n+1),5)==0)||(n==M)))))

        ll =line(x,Uold,'Color','r','Parent',ax);
        ax.Title.String = ['t=',num2str(t(n+1))];
        pause(1)
        %%gif;
        delete(ll)
    end
    Uf=[Uf,U];
end

The solution $u(x,t)$ of PDE with delay

figure(4);
[X , Y] = meshgrid(x,[t]);
mesh (X , Y , Uf')
title('The evolution of solution u(x,t)')
xlabel('space')
ylabel('time')

References

[1] G.Bocharov, A. Meyerhans, N. Bessonov, S. Trofimchuk, V. Volpert. Spatiotemporal dynamics of virus infection spreading in tissues. PlosOne, December 20, 2016.

[2] T. M. Touaoula, M. N. Frioui, N. Bessonov, V. Volpert Dynamics of solutions of a reaction-diffusion equation with delayed inhibition. Manuscript submitted for publication.