numerical simulation of travelling wave of epidemic model

9 views (last 30 days)
lulu
lulu on 10 May 2021
Commented: lulu on 12 May 2021
i want to do the numerical simulation for the results of travelling wave where the infection speed c=sqrt(4av^2(alpha+2lambdav)(Abeta-alpha))/(Abeta+2lambdav) and the steady state(ur,ul,vr,vl) are (A/2,A/2,0,0) and (alpha/2beta,alpha/2beta,A/2-alpha/2beta,A/2-alpha/2beta) how I can combine these results to see the travelling wave of the system : my code is working but how do I incorporate these conditions into Matlab code?
clear all;
xl=0; xr=1; %domain[xl,xr]
J=1000; % J: number of dividion for x
dx=(xr-xl)/J; % dx: mesh size
tf=1000; % final simulation time
Nt=10000; %Nt:number of time steps
dt=tf/Nt;
c=50; % paremeter for the solution
lambdau=0.05; %turning rate of susp
lambdav=0.02; %turning rate of infect
beta=0.01; % transimion rate
alpha=0.001; %recovery/death rate
au=0.01; % advection speed - susceptible
av=0.005; % advection speed - infected
A=6;
mu1=au*dt/dx;
mu2=av*dt/dx;
if mu1>1.0 %make sure dt satisfy stability condition
error('mu1 should<1.0!')
end
if mu2>1.0 %make sure dt satisfy stability condition
error('mu2 should<1.0!')
end
%Evaluate the initial conditions
x=xl:dx:xr; %generate the grid point
fr=exp(-c*(x-0.5).^2); %dimension f(1:J+1)initial conditions of right moving suscp
fl=exp(-c*(x-0.5).^2); %initial conditions of left moving suscep
g=exp(-c*(x-0.4).^2); %initial conditions of infected
%store the solution at all grid points for all time steps
ur=zeros(J+1,Nt);
ul=zeros(J+1,Nt);
vr=zeros(J+1,Nt);
vl=zeros(J+1,Nt);
U=zeros(J+1,Nt);
V=zeros(J+1,Nt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tt=0:tf/(Nt-1):tf;
[T,X]=meshgrid(tt,x);
%find the approximate solution at each time step
for n=1:Nt
t=n*dt; % current time
if n==1 % first time step
for j=2:J % interior nodes
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of susceptible
ur(j,n)=fr(j)-0.5*mu1*(fr(j+1)-fr(j-1))*(1-dt*lambdau+0.5*dt*beta*(g(j)+g(j)))+0.5*(mu1)^2*(fr(j+1)-2*fr(j)+fr(j-1))...
+dt*lambdau*(fl(j)-fr(j))*(1-dt-0.5*dt*beta*(g(j)+g(j)))-dt*beta*fr(j)*(g(j)+g(j))*(1-0.5*dt*beta*(g(j)+g(j)))...
-0.25*mu2*dt*beta*fr(j)*(g(j+1)-g(j-1)-g(j+1)+g(j-1))-0.5*dt^2*beta^2*fr(j)*(g(j)+g(j))*(fr(j)+fl(j))...
+0.5*dt^2*beta*alpha*fr(j)*(g(j)+g(j))+dt*alpha*(g(j)+g(j))-0.125*alpha*mu2*dt*(g(j+1)-g(j-1))+0.125*mu2*alpha*dt*(g(j+1)-g(j-1))...
+0.25*alpha*dt^2*beta*(fr(j)+fl(j))*(g(j)+g(j))-0.25*alpha^2*dt^2*(g(j)+g(j));
%%%%%%%%%%%%%%%%%%%lax wendroff of left-moving of susceptible
ul(j,n)=fl(j)+0.5*mu1*(fl(j+1)-fl(j-1))*(1-dt*lambdau-0.5*dt*beta*(g(j)+g(j)))+0.5*(mu1)^2*(fl(j+1)-2*fl(j)+fl(j-1))...
+dt*lambdau*(fr(j)-fl(j))*(1-dt-0.5*dt*beta*(g(j)+g(j)))-dt*beta*fl(j)*(g(j)+g(j))*(1-0.5*dt*beta*(g(j)+g(j)))...
-0.25*mu2*dt*beta*fl(j)*(g(j+1)-g(j-1)-g(j+1)+g(j-1))-0.5*dt^2*beta^2*fl(j)*(g(j)+g(j))*(fr(j)+fl(j))...
+0.5*dt^2*beta*alpha*fl(j)*(g(j)+g(j))+dt*alpha*(g(j)+g(j))-0.125*alpha*mu2*dt*(g(j+1)-g(j-1))+0.125*mu2*alpha*dt*(g(j+1)-g(j-1))...
+0.25*alpha*dt^2*beta*(fr(j)+fl(j))*(g(j)+g(j))-0.25*alpha^2*dt^2*(g(j)+g(j));
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of infected
vr(j,n)=g(j)-0.5*mu2*(g(j+1)-g(j-1))*(1-dt*lambdav-dt*alpha+dt*beta*fr(j))+0.5*(mu2)^2*(g(j+1)-2*g(j)+g(j-1))...
+dt*(lambdav*(g(j)-g(j))+beta*fr(j)*(g(j)+g(j))-alpha*g(j))-0.25*dt*beta*(g(j)+g(j))*(fr(j+1)-fr(j-1))*(mu1+mu2)...
+dt^2*lambdav*(g(j)-g(j))*(lambdav+alpha)+0.5*dt^2*beta*(g(j)+g(j))*(fl(j)-fr(j))*(lambdau+lambdav)+0.5*dt^2*beta^2*fr(j)*(g(j)+g(j))...
*((fr(j)+fl(j))-(g(j)+g(j)))-dt^2*alpha*beta*fr(j)*(g(j)+g(j))+0.5*dt^2*alpha^2*g(j);
%%%%%%%%%%%%%%%%%%%lax wendroff of left-moving of infected
vl(j,n)=g(j)+0.5*mu2*(g(j+1)-g(j-1))*(1-dt*lambdav-dt*alpha+dt*beta*fl(j))+0.5*(mu2)^2*(g(j+1)-2*g(j)+g(j-1))...
+dt*(lambdav*(g(j)-g(j))+beta*fl(j)*(g(j)+g(j))-alpha*g(j))+0.25*dt*beta*(g(j)+g(j))*(fl(j+1)-fl(j-1))*(mu1+mu2)...
+dt^2*lambdav*(g(j)-g(j))*(lambdav+alpha)+0.5*dt^2*beta*(g(j)+g(j))*(fr(j)-fl(j))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(j)*(g(j)+g(j))...
*((fr(j)+fl(j))-(g(j)+g(j)))-dt^2*alpha*beta*fl(j)*(g(j)+g(j))+0.5*dt^2*alpha^2*g(j);
U(j,n)=fr(j)+fl(j);
V(j,n)=g(j)+g(j);
end
%peridic BC for right moving of suscep(j=J+1,j+2=1)
ur(J+1,1)=fr(J+1)-0.5*mu1*(fr(1)-fr(J))*(1-dt*lambdau+0.5*dt*beta*(g(J+1)+g(J+1)))+0.5*(mu1)^2*(fr(1)-2*fr(J+1)+fr(J))...
+dt*lambdau*(fl(J+1)-fr(J+1))*(1-dt-0.5*dt*beta*(g(J+1)+g(J+1)))-dt*beta*fr(J+1)*(g(J+1)+g(J+1))*(1-0.5*dt*beta*(g(J+1)+g(J+1)))...
-0.25*mu2*dt*beta*fr(J+1)*(g(1)-g(J)-g(1)+g(J))-0.5*dt^2*beta^2*fr(J+1)*(g(J+1)+g(J+1))*(fr(J+1)+fl(J+1))...
+0.5*dt^2*beta*alpha*fr(J+1)*(g(J+1)+g(J+1))+dt*alpha*(g(J+1)+g(J+1))-0.125*alpha*mu2*dt*(g(1)-g(J))+0.125*mu2*alpha*dt*(g(1)-g(J))...
+0.25*alpha*dt^2*beta*(fr(J+1)+fl(J+1))*(g(J+1)+g(J+1))-0.25*alpha^2*dt^2*(g(J+1)+g(J+1)); %lax wendroff ;
%peridic BC for right moving of suscep(j=1,0=J)
ur(1,1)=fr(1)-0.5*mu1*(fr(2)-fr(J+1))*(1-dt*lambdau+0.5*dt*beta*(g(1)+g(1)))+0.5*(mu1)^2*(fr(2)-2*fr(1)+fr(J+1))...
+dt*lambdau*(fl(1)-fr(1))*(1-dt-0.5*dt*beta*(g(1)+g(1)))-dt*beta*fr(j)*(g(1)+g(1))*(1-0.5*dt*beta*(g(1)+g(1)))...
-0.25*mu2*dt*beta*fr(j)*(g(2)-g(J+1)-g(2)+g(J+1))-0.5*dt^2*beta^2*fr(1)*(g(1)+g(1))*(fr(1)+fl(1))...
+0.5*dt^2*beta*alpha*fr(1)*(g(1)+g(1))+dt*alpha*(g(1)+g(1))-0.125*alpha*mu2*dt*(g(2)-g(J))+0.125*mu2*alpha*dt*(g(2)-g(J))...
+0.25*alpha*dt^2*beta*(fr(1)+fl(1))*(g(1)+g(1))-0.25*alpha^2*dt^2*(g(1)+g(1)); %peridic BC for right moving
%peridic BC for left moving of suscep(j=1)
ul(1,1)=fl(1)+0.5*mu1*(fl(2)-fl(J+1))*(1-dt*lambdau-0.5*dt*beta*(g(1)+g(1)))+0.5*(mu1)^2*(fl(2)-2*fl(1)+fl(J+1))...
+dt*lambdau*(fr(1)-fl(1))*(1-dt-0.5*dt*beta*(g(1)+g(1)))-dt*beta*fl(1)*(g(1)+g(1))*(1-0.5*dt*beta*(g(1)+g(1)))...
-0.25*mu2*dt*beta*fl(1)*(g(2)-g(J+1)-g(2)+g(J+1))-0.5*dt^2*beta^2*fl(1)*(g(1)+g(1))*(fr(1)+fl(1))+0.5*dt^2*beta*alpha*fl(1)*(g(1)+g(1)); %peridic BC for left moving; %peridic BC for left moving
%peridic BC for left moving of suscep(j=J+1)
ul(J+1,1)=fl(J+1)+0.5*mu1*(fl(1)-fl(J))*(1-dt*lambdau-0.5*dt*beta*(g(J+1)+g(J+1)))+0.5*(mu1)^2*(fl(1)-2*fl(J+1)+fl(J))...
+dt*lambdau*(fr(J+1)-fl(J+1))*(1-dt-0.5*dt*beta*(g(J+1)+g(J+1)))-dt*beta*fl(J+1)*(g(J+1)+g(J+1))*(1-0.5*dt*beta*(g(J+1)+g(J+1)))...
-0.25*mu2*dt*beta*fl(J+1)*(g(1)-g(J)-g(1)+g(J))-0.5*dt^2*beta^2*fl(J+1)*(g(J+1)+g(J+1))*(fr(J+1)+fl(J+1))+0.5*dt^2*beta*alpha*fl(J+1)*(g(J+1)+g(J+1)); %lax wendroff ; %peridic BC for right moving(j=J+1); %peridic BC for left moving
%peridic BC for right moving of infected(j=J+1)
vr(J+1,n)=g(J+1)-0.5*mu2*(g(1)-g(J))*(1-dt*lambdav-dt*alpha+dt*beta*fr(J+1))+0.5*(mu2)^2*(g(1)-2*g(J+1)+g(J))...
+dt*(lambdav*(g(J+1)-g(J+1))+beta*fr(J+1)*(g(J+1)+g(J+1))-alpha*g(J+1))-0.25*dt*beta*(g(J+1)+g(J+1))*(fr(1)-fr(J))*(mu1+mu2)...
+dt^2*lambdav*(g(J+1)-g(J+1))*(lambdav+alpha)+0.5*dt^2*beta*(g(J+1)+g(J+1))*(fr(J+1)-fl(J+1))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(J+1)*(g(J+1)+g(J+1))...
*((fr(J+1)+fl(J+1))-(g(J+1)+g(J+1)))-dt^2*alpha*beta*fl(J+1)*(g(J+1)+g(J+1))+0.5*dt^2*alpha^2*g(J+1);
%peridic BC for right moving of infected(j=1)
vr(1,n)=g(1)-0.5*mu2*(g(2)-g(J+1))*(1-dt*lambdav-dt*alpha+dt*beta*fr(1))+0.5*(mu2)^2*(g(2)-2*g(1)+g(J+1))...
+dt*(lambdav*(g(1)-g(1))+beta*fr(1)*(g(1)+g(1))-alpha*g(1))-0.25*dt*beta*(g(1)+g(1))*(fr(2)-fr(J+1))*(mu1+mu2)...
+dt^2*lambdav*(g(1)-g(1))*(lambdav+alpha)+0.5*dt^2*beta*(g(1)+g(1))*(fr(1)-fl(1))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(1)*(g(1)+g(1))...
*((fr(1)+fl(1))-(g(1)+g(1)))-dt^2*alpha*beta*fl(1)*(g(1)+g(1))+0.5*dt^2*alpha^2*g(1);
%peridic BC for left moving of infected(j=1)
vl(1,n)=g(1)+0.5*mu2*(g(2)-g(J+1))*(1-dt*lambdav-dt*alpha+dt*beta*fl(1))+0.5*(mu2)^2*(g(j+1)-2*g(j)+g(j-1))...
+dt*(lambdav*(g(1)-g(1))+beta*fl(1)*(g(1)+g(1))-alpha*g(1))+0.25*dt*beta*(g(1)+g(1))*(fl(2)-fl(J+1))*(mu1+mu2)...
+dt^2*lambdav*(g(1)-g(1))*(lambdav+alpha)+0.5*dt^2*beta*(g(1)+g(1))*(fr(1)-fl(1))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(1)*(g(1)+g(1))...
*((fr(1)+fl(1))-(g(1)+g(1)))-dt^2*alpha*beta*fl(1)*(g(1)+g(1))+0.5*dt^2*alpha^2*g(1);
%peridic BC for left moving of infected(j=J+1)
vl(J+1,n)=g(J+1)+0.5*mu2*(g(1)-g(J))*(1-dt*lambdav-dt*alpha+dt*beta*fl(J+1))+0.5*(mu2)^2*(g(1)-2*g(J+1)+g(J))...
+dt*(lambdav*(g(J+1)-g(J+1))+beta*fl(J+1)*(g(J+1)+g(J+1))-alpha*g(J+1))+0.25*dt*beta*(g(J+1)+g(J+1))*(fl(1)-fl(J))*(mu1+mu2)...
+dt^2*lambdav*(g(J+1)-g(J+1))*(lambdav+alpha)+0.5*dt^2*beta*(g(J+1)+g(J+1))*(fr(J+1)-fl(J+1))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(J+1)*(g(J+1)+g(J+1))...
*((fr(J+1)+fl(J+1))-(g(J+1)+g(J+1)))-dt^2*alpha*beta*fl(J+1)*(g(J+1)+g(J+1))+0.5*dt^2*alpha^2*g(J+1);
else
for j=2:J % interior nodes
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of susceptible
ur(j,n)=ur(j,n-1)-0.5*mu1*(ur(j+1,n-1)-ur(j-1,n-1))*(1-dt*lambdau+0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))+0.5*(mu1)^2*(ur(j+1,n-1)-2*ur(j,n-1)+ur(j-1,n-1))...
+dt*lambdau*(ul(j,n-1)-ur(j,n-1))*(1-dt-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))-dt*beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))*(1-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))...
-0.25*mu2*dt*beta*ur(j,n-1)*(vl(j+1,n-1)-vl(j-1,n-1)-vr(j+1,n-1)+vr(j-1,n-1))-0.5*dt^2*beta^2*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))*(ur(j,n-1)+ul(j,n-1))...
+0.5*dt^2*beta*alpha*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))+dt*alpha*(vr(j,n-1)+vl(j,n-1))-0.125*alpha*mu2*dt*(vr(j+1,n-1)-vr(j-1,n-1))+0.125*mu2*alpha*dt*(vl(j+1,n-1)-vl(j-1,n-1))...
+0.25*alpha*dt^2*beta*(ur(j,n-1)+ul(j,n-1))*(vr(j,n-1)+vl(j,n-1))-0.25*alpha^2*dt^2*(vr(j,n-1)+vl(j,n-1));%lax wendroff
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of susceptible
ul(j,n)=ul(j,n-1)+0.5*mu1*(ul(j+1,n-1)-ul(j-1,n-1))*(1-dt*lambdau-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))+0.5*(mu1)^2*(ul(j+1,n-1)-2*ul(j,n-1)+ul(j-1,n-1))...
+dt*lambdau*(ur(j,n-1)-ul(j,n-1))*(1-dt-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))-dt*beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))*(1-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))...
-0.25*mu2*dt*beta*ul(j,n-1)*(vl(j+1,n-1)-vl(j-1,n-1)-vr(j+1,n-1)+vr(j-1,n-1))-0.5*dt^2*beta^2*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))*(ur(j,n-1)+ul(j,n-1))...
+0.5*dt^2*beta*alpha*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))+dt*alpha*(vr(j,n-1)+vl(j,n-1))-0.125*alpha*mu2*dt*(vr(j+1,n-1)-vr(j-1,n-1))+0.125*mu2*alpha*dt*(vl(j+1,n-1)-vl(j-1,n-1))...
+0.25*alpha*dt^2*beta*(ur(j,n-1)+ul(j,n-1))*(vr(j,n-1)+vl(j,n-1))-0.25*alpha^2*dt^2*(vr(j,n-1)+vl(j,n-1));
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of infected
vr(j,n)=vr(j,n-1)-0.5*mu2*(vr(j+1,n-1)-vr(j-1,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ur(j,n-1))+0.5*(mu2)^2*(vr(j+1,n-1)-2*vr(j,n-1)+vr(j-1,n-1))...
+dt*(lambdav*(vl(j,n-1)-vr(j,n-1))+beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))-alpha*vr(j,n-1))-0.25*dt*beta*(vr(j,n-1)+vl(j,n-1))*(ur(j+1,n-1)-ur(j-1,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vr(j,n-1)-vl(j,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(j,n-1)+vl(j,n-1))*(ul(j,n-1)-ur(j,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))...
*((ur(j,n-1)+ul(j,n-1))-(vr(j,n-1)+vl(j,n-1)))-dt^2*alpha*beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))+0.5*dt^2*alpha^2*vr(j,n-1);
%%%%%%%%%%%%%%%%%%%lax wendroff of left-moving of infected
vl(j,n)=vl(j,n-1)+0.5*mu2*(vl(j+1,n-1)-vl(j-1,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ul(j,n-1))+0.5*(mu2)^2*(vl(j+1,n-1)-2*vl(j,n-1)+vl(j-1,n-1))...
+dt*(lambdav*(vr(j,n-1)-vl(j,n-1))+beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))-alpha*vl(j,n-1))+0.25*dt*beta*(vr(j,n-1)+vl(j,n-1))*(ul(j+1,n-1)-ul(j-1,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vl(j,n-1)-vr(j,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(j,n-1)+vl(j,n-1))*(ur(j,n-1)-ul(j,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))...
*((ur(j,n-1)+ul(j,n-1))-(vr(j,n-1)+vl(j,n-1)))-dt^2*alpha*beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))+0.5*dt^2*alpha^2*vl(j,n-1);
U(j,n)=ur(j,n)+ul(j,n);
V(j,n)=vr(j,n)+vl(j,n);
end
%peridic BC for right moving of suscep(j=J+1)
ur(J+1,n)=ur(J+1,n-1)-0.5*mu1*(ur(1,n-1)-ur(J,n-1))*(1-dt*lambdau+0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1)))+0.5*(mu1)^2*(ur(1,n-1)-2*ur(J+1,n-1)+ur(J,n-1))...
+dt*lambdau*(ul(J+1,n-1)-ur(J+1,n-1))*(1-dt-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1))-dt*beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))*(1-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1)))...
-0.25*mu2*dt*beta*ur(J+1,n-1)*(vl(1,n-1)-vl(J,n-1)-vr(1,n-1)+vr(J,n-1))-0.5*dt^2*beta^2*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))*(ur(J+1,n-1)+ul(J+1,n-1))...
+0.5*dt^2*beta*alpha*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1)))+dt*alpha*(vr(J+1,n-1)+vl(J+1,n-1))-0.125*alpha*mu2*dt*(vr(1,n-1)-vr(J,n-1))+0.125*mu2*alpha*dt*(vl(1,n-1)-vl(J,n-1))...
+0.25*alpha*dt^2*beta*(ur(J+1,n-1)+ul(J+1,n-1))*(vr(J+1,n-1)+vl(J+1,n-1))-0.25*alpha^2*dt^2*(vr(J+1,n-1)+vl(J+1,n-1));
%peridic BC for right moving of suscep(j=1)
ur(1,n)=ur(1,n-1)-0.5*mu1*(ur(2,n-1)-ur(J+1,n-1))*(1-dt*lambdau+0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))+0.5*(mu1)^2*(ur(2,n-1)-2*ur(1,n-1)+ur(J+1,n-1))...
+dt*lambdau*(ul(1,n-1)-ur(1,n-1))*(1-dt-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))-dt*beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))*(1-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))...
-0.25*mu2*dt*beta*ur(j,n-1)*(vl(j+1,n-1)-vl(J+1,n-1)-vr(2,n-1)+vr(J+1,n-1))-0.5*dt^2*beta^2*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))*(ur(1,n-1)+ul(1,n-1))...
+0.5*dt^2*beta*alpha*ur(1,n-1)*(vr(1,n-1)+vl(j,n-1))+dt*alpha*(vr(1,n-1)+vl(1,n-1))-0.125*alpha*mu2*dt*(vr(2,n-1)-vr(J,n-1))+0.125*mu2*alpha*dt*(vl(2,n-1)-vl(J,n-1))...
+0.25*alpha*dt^2*beta*(ur(1,n-1)+ul(1,n-1))*(vr(1,n-1)+vl(1,n-1))-0.25*alpha^2*dt^2*(vr(1,n-1)+vl(1,n-1));%lax wendroff
%peridic BC for left moving of suscep (j=1)
ul(1,n)=ul(1,n-1)+0.5*mu1*(ul(2,n-1)-ul(J+1,n-1))*(1-dt*lambdau-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))+0.5*(mu1)^2*(ul(2,n-1)-2*ul(1,n-1)+ul(J+1,n-1))...
+dt*lambdau*(ur(1,n-1)-ul(1,n-1))*(1-dt-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))-dt*beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))*(1-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))...
-0.25*mu2*dt*beta*ul(1,n-1)*(vl(2,n-1)-vl(J+1,n-1)-vr(2,n-1)+vr(J+1,n-1))-0.5*dt^2*beta^2*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))*(ur(1,n-1)+ul(1,n-1))+0.5*dt^2*beta*alpha*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1));
%peridic BC for left moving of suscep(j=J+1)
ul(J+1,n)=ul(J+1,n-1)+0.5*mu1*(ul(1,n-1)-ul(J,n-1))*(1-dt*lambdau-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1)))+0.5*(mu1)^2*(ul(1,n-1)-2*ul(J+1,n-1)+ul(J,n-1))...
+dt*lambdau*(ur(J+1,n-1)-ul(J+1,n-1))*(1-dt-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1))-dt*beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))*(1-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1)))...
-0.25*mu2*dt*beta*ul(J+1,n-1)*(vl(1,n-1)-vl(J,n-1)-vr(1,n-1)+vr(J,n-1))-0.5*dt^2*beta^2*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))*(ur(J+1,n-1)+ul(J+1,n-1))+0.5*dt^2*beta*alpha*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1)));
%peridic BC for right moving of infected(j=J+1)
vr(J+1,n)=vr(J+1,n-1)-0.5*mu2*(vr(1,n-1)-vr(J,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ur(J+1,n-1))+0.5*(mu2)^2*(vr(1,n-1)-2*vr(J+1,n-1)+vr(J,n-1))...
+dt*(lambdav*(vl(J+1,n-1)-vr(J+1,n-1))+beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))-alpha*vr(J+1,n-1))-0.25*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1))*(ur(1,n-1)-ur(J,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vr(J+1,n-1)-vl(J+1,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(J+1,n-1)+vl(J+1,n-1))*(ul(J+1,n-1)-ur(J+1,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))...
*((ur(J+1,n-1)+ul(J+1,n-1))-(vr(J+1,n-1)+vl(J+1,n-1)))-dt^2*alpha*beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))+0.5*dt^2*alpha^2*vr(J+1,n-1);
%peridic BC for right moving of infected(j=1)
vr(1,n)=vr(1,n-1)-0.5*mu2*(vr(2,n-1)-vr(J+1,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ur(1,n-1))+0.5*(mu2)^2*(vr(2,n-1)-2*vr(1,n-1)+vr(J+1,n-1))...
+dt*(lambdav*(vl(1,n-1)-vr(1,n-1))+beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))-alpha*vr(1,n-1))-0.25*dt*beta*(vr(1,n-1)+vl(1,n-1))*(ur(2,n-1)-ur(J+1,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vr(1,n-1)-vl(1,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(1,n-1)+vl(1,n-1))*(ul(1,n-1)-ur(1,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))...
*((ur(1,n-1)+ul(1,n-1))-(vr(1,n-1)+vl(1,n-1)))-dt^2*alpha*beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))+0.5*dt^2*alpha^2*vr(1,n-1);
%peridic BC for left moving of infected(j=1)
vl(1,n)=vl(1,n-1)+0.5*mu2*(vl(2,n-1)-vl(J+1,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ul(1,n-1))+0.5*(mu2)^2*(vl(2,n-1)-2*vl(1,n-1)+vl(J+1,n-1))...
+dt*(lambdav*(vr(1,n-1)-vl(1,n-1))+beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))-alpha*vl(1,n-1))+0.25*dt*beta*(vr(1,n-1)+vl(1,n-1))*(ul(2,n-1)-ul(J+1,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vl(1,n-1)-vr(1,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(1,n-1)+vl(1,n-1))*(ur(1,n-1)-ul(1,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))...
*((ur(1,n-1)+ul(1,n-1))-(vr(1,n-1)+vl(1,n-1)))-dt^2*alpha*beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))+0.5*dt^2*alpha^2*vl(1,n-1);
%peridic BC for left moving of infected(j=J+1)
vl(J+1,n)=vl(J+1,n-1)+0.5*mu2*(vl(1,n-1)-vl(J,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ul(J+1,n-1))+0.5*(mu2)^2*(vl(1,n-1)-2*vl(J+1,n-1)+vl(J,n-1))...
+dt*(lambdav*(vr(J+1,n-1)-vl(J+1,n-1))+beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))-alpha*vl(J+1,n-1))+0.25*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1))*(ul(1,n-1)-ul(J,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vl(J+1,n-1)-vr(J+1,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(J+1,n-1)+vl(J+1,n-1))*(ur(J+1,n-1)-ul(J+1,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))...
*((ur(J+1,n-1)+ul(J+1,n-1))-(vr(J+1,n-1)+vl(J+1,n-1)))-dt^2*alpha*beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))+0.5*dt^2*alpha^2*vl(J+1,n-1);
end
%calculate the analytical solution
for j=1:J+1
xj=xl+(j-1)*dx;
u_ex(j,n)=exp(-c*(xj-au*t-0.5)^2);
end
end
figure;
hold on;
n=1;
plot(x,ur(:,n),'r-','linewidth',2);
plot(x,ul(:,n),'r--','linewidth',2);
plot(x,vr(:,n),'r:','linewidth',2);
plot(x,vl(:,n),'r-.','linewidth',2);
n=100;
% n=2000;
plot(x,ur(:,n),'g-','linewidth',2);
plot(x,ul(:,n),'g--','linewidth',2);
plot(x,vr(:,n),'g:','linewidth',2);
plot(x,vl(:,n),'g-.','linewidth',2);
%%%%%%%%%%%%%%%%%
n=200;
plot(x,ur(:,n),'b-','linewidth',2); %blue
plot(x,ul(:,n),'b--','linewidth',2);
plot(x,vr(:,n),'b:','linewidth',2);
plot(x,vl(:,n),'b-.','linewidth',2);
%%%%%%%%%%%%%%%%%%%%%%% xlabel('x','Fontsize',20,'FontWeight','bold','Color','k');
% %%% time=n*\deta t=
legend('u^+ t=0','u^- t=0','v^+ t=0','v^- t=0','u^+ t=10','u^- t=10','v^+ t=10','v^- t=10','u^+ t=20','u^- t=20','v^+ t=20','v^- t=20');% title('Numerical and Analytical Solutions at t=0.1,0.3,0.5');
xlabel('x','Fontsize',20,'FontWeight','bold','Color','k');
ylabel('u(x,t),v(x,t)','Fontsize',20,'FontWeight','bold','Color','k');
set(gca,'Fontsize',18);

Sign in to answer this question.

Answers (1)

Rafael Hernandez-Walls
Rafael Hernandez-Walls on 11 May 2021
try this....
clear all;
xl=0; xr=1; %domain[xl,xr]
J=1000; % J: number of dividion for x
dx=(xr-xl)/J; % dx: mesh size
tf=1000; % final simulation time
Nt=10000; %Nt:number of time steps
dt=tf/Nt;
c=50; % paremeter for the solution
lambdau=0.05; %turning rate of susp
lambdav=0.02; %turning rate of infect
beta=0.01; % transimion rate
alpha=0.001; %recovery/death rate
au=0.01; % advection speed - susceptible
av=0.005; % advection speed - infected
A=6;
mu1=au*dt/dx;
mu2=av*dt/dx;
if mu1>1.0 %make sure dt satisfy stability condition
error('mu1 should<1.0!')
end
if mu2>1.0 %make sure dt satisfy stability condition
error('mu2 should<1.0!')
end
%Evaluate the initial conditions
x=xl:dx:xr; %generate the grid point
fr=exp(-c*(x-0.5).^2); %dimension f(1:J+1)initial conditions of right moving suscp
fl=exp(-c*(x-0.5).^2); %initial conditions of left moving suscep
g=exp(-c*(x-0.4).^2); %initial conditions of infected
%store the solution at all grid points for all time steps
ur=zeros(J+1,Nt);
ul=zeros(J+1,Nt);
vr=zeros(J+1,Nt);
vl=zeros(J+1,Nt);
U=zeros(J+1,Nt);
V=zeros(J+1,Nt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tt=0:tf/(Nt-1):tf;
[T,X]=meshgrid(tt,x);
%find the approximate solution at each time step
for n=1:Nt
t=n*dt; % current time
if n==1 % first time step
for j=2:J % interior nodes
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of susceptible
ur(j,n)=fr(j)-0.5*mu1*(fr(j+1)-fr(j-1))*(1-dt*lambdau+0.5*dt*beta*(g(j)+g(j)))+0.5*(mu1)^2*(fr(j+1)-2*fr(j)+fr(j-1))...
+dt*lambdau*(fl(j)-fr(j))*(1-dt-0.5*dt*beta*(g(j)+g(j)))-dt*beta*fr(j)*(g(j)+g(j))*(1-0.5*dt*beta*(g(j)+g(j)))...
-0.25*mu2*dt*beta*fr(j)*(g(j+1)-g(j-1)-g(j+1)+g(j-1))-0.5*dt^2*beta^2*fr(j)*(g(j)+g(j))*(fr(j)+fl(j))...
+0.5*dt^2*beta*alpha*fr(j)*(g(j)+g(j))+dt*alpha*(g(j)+g(j))-0.125*alpha*mu2*dt*(g(j+1)-g(j-1))+0.125*mu2*alpha*dt*(g(j+1)-g(j-1))...
+0.25*alpha*dt^2*beta*(fr(j)+fl(j))*(g(j)+g(j))-0.25*alpha^2*dt^2*(g(j)+g(j));
%%%%%%%%%%%%%%%%%%%lax wendroff of left-moving of susceptible
ul(j,n)=fl(j)+0.5*mu1*(fl(j+1)-fl(j-1))*(1-dt*lambdau-0.5*dt*beta*(g(j)+g(j)))+0.5*(mu1)^2*(fl(j+1)-2*fl(j)+fl(j-1))...
+dt*lambdau*(fr(j)-fl(j))*(1-dt-0.5*dt*beta*(g(j)+g(j)))-dt*beta*fl(j)*(g(j)+g(j))*(1-0.5*dt*beta*(g(j)+g(j)))...
-0.25*mu2*dt*beta*fl(j)*(g(j+1)-g(j-1)-g(j+1)+g(j-1))-0.5*dt^2*beta^2*fl(j)*(g(j)+g(j))*(fr(j)+fl(j))...
+0.5*dt^2*beta*alpha*fl(j)*(g(j)+g(j))+dt*alpha*(g(j)+g(j))-0.125*alpha*mu2*dt*(g(j+1)-g(j-1))+0.125*mu2*alpha*dt*(g(j+1)-g(j-1))...
+0.25*alpha*dt^2*beta*(fr(j)+fl(j))*(g(j)+g(j))-0.25*alpha^2*dt^2*(g(j)+g(j));
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of infected
vr(j,n)=g(j)-0.5*mu2*(g(j+1)-g(j-1))*(1-dt*lambdav-dt*alpha+dt*beta*fr(j))+0.5*(mu2)^2*(g(j+1)-2*g(j)+g(j-1))...
+dt*(lambdav*(g(j)-g(j))+beta*fr(j)*(g(j)+g(j))-alpha*g(j))-0.25*dt*beta*(g(j)+g(j))*(fr(j+1)-fr(j-1))*(mu1+mu2)...
+dt^2*lambdav*(g(j)-g(j))*(lambdav+alpha)+0.5*dt^2*beta*(g(j)+g(j))*(fl(j)-fr(j))*(lambdau+lambdav)+0.5*dt^2*beta^2*fr(j)*(g(j)+g(j))...
*((fr(j)+fl(j))-(g(j)+g(j)))-dt^2*alpha*beta*fr(j)*(g(j)+g(j))+0.5*dt^2*alpha^2*g(j);
%%%%%%%%%%%%%%%%%%%lax wendroff of left-moving of infected
vl(j,n)=g(j)+0.5*mu2*(g(j+1)-g(j-1))*(1-dt*lambdav-dt*alpha+dt*beta*fl(j))+0.5*(mu2)^2*(g(j+1)-2*g(j)+g(j-1))...
+dt*(lambdav*(g(j)-g(j))+beta*fl(j)*(g(j)+g(j))-alpha*g(j))+0.25*dt*beta*(g(j)+g(j))*(fl(j+1)-fl(j-1))*(mu1+mu2)...
+dt^2*lambdav*(g(j)-g(j))*(lambdav+alpha)+0.5*dt^2*beta*(g(j)+g(j))*(fr(j)-fl(j))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(j)*(g(j)+g(j))...
*((fr(j)+fl(j))-(g(j)+g(j)))-dt^2*alpha*beta*fl(j)*(g(j)+g(j))+0.5*dt^2*alpha^2*g(j);
U(j,n)=fr(j)+fl(j);
V(j,n)=g(j)+g(j);
end
%peridic BC for right moving of suscep(j=J+1,j+2=1)
ur(J+1,1)=fr(J+1)-0.5*mu1*(fr(1)-fr(J))*(1-dt*lambdau+0.5*dt*beta*(g(J+1)+g(J+1)))+0.5*(mu1)^2*(fr(1)-2*fr(J+1)+fr(J))...
+dt*lambdau*(fl(J+1)-fr(J+1))*(1-dt-0.5*dt*beta*(g(J+1)+g(J+1)))-dt*beta*fr(J+1)*(g(J+1)+g(J+1))*(1-0.5*dt*beta*(g(J+1)+g(J+1)))...
-0.25*mu2*dt*beta*fr(J+1)*(g(1)-g(J)-g(1)+g(J))-0.5*dt^2*beta^2*fr(J+1)*(g(J+1)+g(J+1))*(fr(J+1)+fl(J+1))...
+0.5*dt^2*beta*alpha*fr(J+1)*(g(J+1)+g(J+1))+dt*alpha*(g(J+1)+g(J+1))-0.125*alpha*mu2*dt*(g(1)-g(J))+0.125*mu2*alpha*dt*(g(1)-g(J))...
+0.25*alpha*dt^2*beta*(fr(J+1)+fl(J+1))*(g(J+1)+g(J+1))-0.25*alpha^2*dt^2*(g(J+1)+g(J+1)); %lax wendroff ;
%peridic BC for right moving of suscep(j=1,0=J)
ur(1,1)=fr(1)-0.5*mu1*(fr(2)-fr(J+1))*(1-dt*lambdau+0.5*dt*beta*(g(1)+g(1)))+0.5*(mu1)^2*(fr(2)-2*fr(1)+fr(J+1))...
+dt*lambdau*(fl(1)-fr(1))*(1-dt-0.5*dt*beta*(g(1)+g(1)))-dt*beta*fr(j)*(g(1)+g(1))*(1-0.5*dt*beta*(g(1)+g(1)))...
-0.25*mu2*dt*beta*fr(j)*(g(2)-g(J+1)-g(2)+g(J+1))-0.5*dt^2*beta^2*fr(1)*(g(1)+g(1))*(fr(1)+fl(1))...
+0.5*dt^2*beta*alpha*fr(1)*(g(1)+g(1))+dt*alpha*(g(1)+g(1))-0.125*alpha*mu2*dt*(g(2)-g(J))+0.125*mu2*alpha*dt*(g(2)-g(J))...
+0.25*alpha*dt^2*beta*(fr(1)+fl(1))*(g(1)+g(1))-0.25*alpha^2*dt^2*(g(1)+g(1)); %peridic BC for right moving
%peridic BC for left moving of suscep(j=1)
ul(1,1)=fl(1)+0.5*mu1*(fl(2)-fl(J+1))*(1-dt*lambdau-0.5*dt*beta*(g(1)+g(1)))+0.5*(mu1)^2*(fl(2)-2*fl(1)+fl(J+1))...
+dt*lambdau*(fr(1)-fl(1))*(1-dt-0.5*dt*beta*(g(1)+g(1)))-dt*beta*fl(1)*(g(1)+g(1))*(1-0.5*dt*beta*(g(1)+g(1)))...
-0.25*mu2*dt*beta*fl(1)*(g(2)-g(J+1)-g(2)+g(J+1))-0.5*dt^2*beta^2*fl(1)*(g(1)+g(1))*(fr(1)+fl(1))+0.5*dt^2*beta*alpha*fl(1)*(g(1)+g(1)); %peridic BC for left moving; %peridic BC for left moving
%peridic BC for left moving of suscep(j=J+1)
ul(J+1,1)=fl(J+1)+0.5*mu1*(fl(1)-fl(J))*(1-dt*lambdau-0.5*dt*beta*(g(J+1)+g(J+1)))+0.5*(mu1)^2*(fl(1)-2*fl(J+1)+fl(J))...
+dt*lambdau*(fr(J+1)-fl(J+1))*(1-dt-0.5*dt*beta*(g(J+1)+g(J+1)))-dt*beta*fl(J+1)*(g(J+1)+g(J+1))*(1-0.5*dt*beta*(g(J+1)+g(J+1)))...
-0.25*mu2*dt*beta*fl(J+1)*(g(1)-g(J)-g(1)+g(J))-0.5*dt^2*beta^2*fl(J+1)*(g(J+1)+g(J+1))*(fr(J+1)+fl(J+1))+0.5*dt^2*beta*alpha*fl(J+1)*(g(J+1)+g(J+1)); %lax wendroff ; %peridic BC for right moving(j=J+1); %peridic BC for left moving
%peridic BC for right moving of infected(j=J+1)
vr(J+1,n)=g(J+1)-0.5*mu2*(g(1)-g(J))*(1-dt*lambdav-dt*alpha+dt*beta*fr(J+1))+0.5*(mu2)^2*(g(1)-2*g(J+1)+g(J))...
+dt*(lambdav*(g(J+1)-g(J+1))+beta*fr(J+1)*(g(J+1)+g(J+1))-alpha*g(J+1))-0.25*dt*beta*(g(J+1)+g(J+1))*(fr(1)-fr(J))*(mu1+mu2)...
+dt^2*lambdav*(g(J+1)-g(J+1))*(lambdav+alpha)+0.5*dt^2*beta*(g(J+1)+g(J+1))*(fr(J+1)-fl(J+1))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(J+1)*(g(J+1)+g(J+1))...
*((fr(J+1)+fl(J+1))-(g(J+1)+g(J+1)))-dt^2*alpha*beta*fl(J+1)*(g(J+1)+g(J+1))+0.5*dt^2*alpha^2*g(J+1);
%peridic BC for right moving of infected(j=1)
vr(1,n)=g(1)-0.5*mu2*(g(2)-g(J+1))*(1-dt*lambdav-dt*alpha+dt*beta*fr(1))+0.5*(mu2)^2*(g(2)-2*g(1)+g(J+1))...
+dt*(lambdav*(g(1)-g(1))+beta*fr(1)*(g(1)+g(1))-alpha*g(1))-0.25*dt*beta*(g(1)+g(1))*(fr(2)-fr(J+1))*(mu1+mu2)...
+dt^2*lambdav*(g(1)-g(1))*(lambdav+alpha)+0.5*dt^2*beta*(g(1)+g(1))*(fr(1)-fl(1))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(1)*(g(1)+g(1))...
*((fr(1)+fl(1))-(g(1)+g(1)))-dt^2*alpha*beta*fl(1)*(g(1)+g(1))+0.5*dt^2*alpha^2*g(1);
%peridic BC for left moving of infected(j=1)
vl(1,n)=g(1)+0.5*mu2*(g(2)-g(J+1))*(1-dt*lambdav-dt*alpha+dt*beta*fl(1))+0.5*(mu2)^2*(g(j+1)-2*g(j)+g(j-1))...
+dt*(lambdav*(g(1)-g(1))+beta*fl(1)*(g(1)+g(1))-alpha*g(1))+0.25*dt*beta*(g(1)+g(1))*(fl(2)-fl(J+1))*(mu1+mu2)...
+dt^2*lambdav*(g(1)-g(1))*(lambdav+alpha)+0.5*dt^2*beta*(g(1)+g(1))*(fr(1)-fl(1))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(1)*(g(1)+g(1))...
*((fr(1)+fl(1))-(g(1)+g(1)))-dt^2*alpha*beta*fl(1)*(g(1)+g(1))+0.5*dt^2*alpha^2*g(1);
%peridic BC for left moving of infected(j=J+1)
vl(J+1,n)=g(J+1)+0.5*mu2*(g(1)-g(J))*(1-dt*lambdav-dt*alpha+dt*beta*fl(J+1))+0.5*(mu2)^2*(g(1)-2*g(J+1)+g(J))...
+dt*(lambdav*(g(J+1)-g(J+1))+beta*fl(J+1)*(g(J+1)+g(J+1))-alpha*g(J+1))+0.25*dt*beta*(g(J+1)+g(J+1))*(fl(1)-fl(J))*(mu1+mu2)...
+dt^2*lambdav*(g(J+1)-g(J+1))*(lambdav+alpha)+0.5*dt^2*beta*(g(J+1)+g(J+1))*(fr(J+1)-fl(J+1))*(lambdau+lambdav)+0.5*dt^2*beta^2*fl(J+1)*(g(J+1)+g(J+1))...
*((fr(J+1)+fl(J+1))-(g(J+1)+g(J+1)))-dt^2*alpha*beta*fl(J+1)*(g(J+1)+g(J+1))+0.5*dt^2*alpha^2*g(J+1);
else
for j=2:J % interior nodes
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of susceptible
ur(j,n)=ur(j,n-1)-0.5*mu1*(ur(j+1,n-1)-ur(j-1,n-1))*(1-dt*lambdau+0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))+0.5*(mu1)^2*(ur(j+1,n-1)-2*ur(j,n-1)+ur(j-1,n-1))...
+dt*lambdau*(ul(j,n-1)-ur(j,n-1))*(1-dt-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))-dt*beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))*(1-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))...
-0.25*mu2*dt*beta*ur(j,n-1)*(vl(j+1,n-1)-vl(j-1,n-1)-vr(j+1,n-1)+vr(j-1,n-1))-0.5*dt^2*beta^2*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))*(ur(j,n-1)+ul(j,n-1))...
+0.5*dt^2*beta*alpha*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))+dt*alpha*(vr(j,n-1)+vl(j,n-1))-0.125*alpha*mu2*dt*(vr(j+1,n-1)-vr(j-1,n-1))+0.125*mu2*alpha*dt*(vl(j+1,n-1)-vl(j-1,n-1))...
+0.25*alpha*dt^2*beta*(ur(j,n-1)+ul(j,n-1))*(vr(j,n-1)+vl(j,n-1))-0.25*alpha^2*dt^2*(vr(j,n-1)+vl(j,n-1));%lax wendroff
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of susceptible
ul(j,n)=ul(j,n-1)+0.5*mu1*(ul(j+1,n-1)-ul(j-1,n-1))*(1-dt*lambdau-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))+0.5*(mu1)^2*(ul(j+1,n-1)-2*ul(j,n-1)+ul(j-1,n-1))...
+dt*lambdau*(ur(j,n-1)-ul(j,n-1))*(1-dt-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))-dt*beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))*(1-0.5*dt*beta*(vr(j,n-1)+vl(j,n-1)))...
-0.25*mu2*dt*beta*ul(j,n-1)*(vl(j+1,n-1)-vl(j-1,n-1)-vr(j+1,n-1)+vr(j-1,n-1))-0.5*dt^2*beta^2*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))*(ur(j,n-1)+ul(j,n-1))...
+0.5*dt^2*beta*alpha*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))+dt*alpha*(vr(j,n-1)+vl(j,n-1))-0.125*alpha*mu2*dt*(vr(j+1,n-1)-vr(j-1,n-1))+0.125*mu2*alpha*dt*(vl(j+1,n-1)-vl(j-1,n-1))...
+0.25*alpha*dt^2*beta*(ur(j,n-1)+ul(j,n-1))*(vr(j,n-1)+vl(j,n-1))-0.25*alpha^2*dt^2*(vr(j,n-1)+vl(j,n-1));
%%%%%%%%%%%%%%%%%%%lax wendroff of right-moving of infected
vr(j,n)=vr(j,n-1)-0.5*mu2*(vr(j+1,n-1)-vr(j-1,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ur(j,n-1))+0.5*(mu2)^2*(vr(j+1,n-1)-2*vr(j,n-1)+vr(j-1,n-1))...
+dt*(lambdav*(vl(j,n-1)-vr(j,n-1))+beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))-alpha*vr(j,n-1))-0.25*dt*beta*(vr(j,n-1)+vl(j,n-1))*(ur(j+1,n-1)-ur(j-1,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vr(j,n-1)-vl(j,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(j,n-1)+vl(j,n-1))*(ul(j,n-1)-ur(j,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))...
*((ur(j,n-1)+ul(j,n-1))-(vr(j,n-1)+vl(j,n-1)))-dt^2*alpha*beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))+0.5*dt^2*alpha^2*vr(j,n-1);
%%%%%%%%%%%%%%%%%%%lax wendroff of left-moving of infected
vl(j,n)=vl(j,n-1)+0.5*mu2*(vl(j+1,n-1)-vl(j-1,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ul(j,n-1))+0.5*(mu2)^2*(vl(j+1,n-1)-2*vl(j,n-1)+vl(j-1,n-1))...
+dt*(lambdav*(vr(j,n-1)-vl(j,n-1))+beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))-alpha*vl(j,n-1))+0.25*dt*beta*(vr(j,n-1)+vl(j,n-1))*(ul(j+1,n-1)-ul(j-1,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vl(j,n-1)-vr(j,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(j,n-1)+vl(j,n-1))*(ur(j,n-1)-ul(j,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))...
*((ur(j,n-1)+ul(j,n-1))-(vr(j,n-1)+vl(j,n-1)))-dt^2*alpha*beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))+0.5*dt^2*alpha^2*vl(j,n-1);
U(j,n)=ur(j,n)+ul(j,n);
V(j,n)=vr(j,n)+vl(j,n);
end
%peridic BC for right moving of suscep(j=J+1)
ur(J+1,n)=ur(J+1,n-1)-0.5*mu1*(ur(1,n-1)-ur(J,n-1))*(1-dt*lambdau+0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1)))+0.5*(mu1)^2*(ur(1,n-1)-2*ur(J+1,n-1)+ur(J,n-1))...
+dt*lambdau*(ul(J+1,n-1)-ur(J+1,n-1))*(1-dt-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1))-dt*beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))*(1-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1)))...
-0.25*mu2*dt*beta*ur(J+1,n-1)*(vl(1,n-1)-vl(J,n-1)-vr(1,n-1)+vr(J,n-1))-0.5*dt^2*beta^2*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))*(ur(J+1,n-1)+ul(J+1,n-1))...
+0.5*dt^2*beta*alpha*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1)))+dt*alpha*(vr(J+1,n-1)+vl(J+1,n-1))-0.125*alpha*mu2*dt*(vr(1,n-1)-vr(J,n-1))+0.125*mu2*alpha*dt*(vl(1,n-1)-vl(J,n-1))...
+0.25*alpha*dt^2*beta*(ur(J+1,n-1)+ul(J+1,n-1))*(vr(J+1,n-1)+vl(J+1,n-1))-0.25*alpha^2*dt^2*(vr(J+1,n-1)+vl(J+1,n-1));
%peridic BC for right moving of suscep(j=1)
ur(1,n)=ur(1,n-1)-0.5*mu1*(ur(2,n-1)-ur(J+1,n-1))*(1-dt*lambdau+0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))+0.5*(mu1)^2*(ur(2,n-1)-2*ur(1,n-1)+ur(J+1,n-1))...
+dt*lambdau*(ul(1,n-1)-ur(1,n-1))*(1-dt-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))-dt*beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))*(1-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))...
-0.25*mu2*dt*beta*ur(j,n-1)*(vl(j+1,n-1)-vl(J+1,n-1)-vr(2,n-1)+vr(J+1,n-1))-0.5*dt^2*beta^2*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))*(ur(1,n-1)+ul(1,n-1))...
+0.5*dt^2*beta*alpha*ur(1,n-1)*(vr(1,n-1)+vl(j,n-1))+dt*alpha*(vr(1,n-1)+vl(1,n-1))-0.125*alpha*mu2*dt*(vr(2,n-1)-vr(J,n-1))+0.125*mu2*alpha*dt*(vl(2,n-1)-vl(J,n-1))...
+0.25*alpha*dt^2*beta*(ur(1,n-1)+ul(1,n-1))*(vr(1,n-1)+vl(1,n-1))-0.25*alpha^2*dt^2*(vr(1,n-1)+vl(1,n-1));%lax wendroff
%peridic BC for left moving of suscep (j=1)
ul(1,n)=ul(1,n-1)+0.5*mu1*(ul(2,n-1)-ul(J+1,n-1))*(1-dt*lambdau-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))+0.5*(mu1)^2*(ul(2,n-1)-2*ul(1,n-1)+ul(J+1,n-1))...
+dt*lambdau*(ur(1,n-1)-ul(1,n-1))*(1-dt-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))-dt*beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))*(1-0.5*dt*beta*(vr(1,n-1)+vl(1,n-1)))...
-0.25*mu2*dt*beta*ul(1,n-1)*(vl(2,n-1)-vl(J+1,n-1)-vr(2,n-1)+vr(J+1,n-1))-0.5*dt^2*beta^2*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))*(ur(1,n-1)+ul(1,n-1))+0.5*dt^2*beta*alpha*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1));
%peridic BC for left moving of suscep(j=J+1)
ul(J+1,n)=ul(J+1,n-1)+0.5*mu1*(ul(1,n-1)-ul(J,n-1))*(1-dt*lambdau-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1)))+0.5*(mu1)^2*(ul(1,n-1)-2*ul(J+1,n-1)+ul(J,n-1))...
+dt*lambdau*(ur(J+1,n-1)-ul(J+1,n-1))*(1-dt-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1))-dt*beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))*(1-0.5*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1)))...
-0.25*mu2*dt*beta*ul(J+1,n-1)*(vl(1,n-1)-vl(J,n-1)-vr(1,n-1)+vr(J,n-1))-0.5*dt^2*beta^2*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))*(ur(J+1,n-1)+ul(J+1,n-1))+0.5*dt^2*beta*alpha*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1)));
%peridic BC for right moving of infected(j=J+1)
vr(J+1,n)=vr(J+1,n-1)-0.5*mu2*(vr(1,n-1)-vr(J,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ur(J+1,n-1))+0.5*(mu2)^2*(vr(1,n-1)-2*vr(J+1,n-1)+vr(J,n-1))...
+dt*(lambdav*(vl(J+1,n-1)-vr(J+1,n-1))+beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))-alpha*vr(J+1,n-1))-0.25*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1))*(ur(1,n-1)-ur(J,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vr(J+1,n-1)-vl(J+1,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(J+1,n-1)+vl(J+1,n-1))*(ul(J+1,n-1)-ur(J+1,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))...
*((ur(J+1,n-1)+ul(J+1,n-1))-(vr(J+1,n-1)+vl(J+1,n-1)))-dt^2*alpha*beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))+0.5*dt^2*alpha^2*vr(J+1,n-1);
%peridic BC for right moving of infected(j=1)
vr(1,n)=vr(1,n-1)-0.5*mu2*(vr(2,n-1)-vr(J+1,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ur(1,n-1))+0.5*(mu2)^2*(vr(2,n-1)-2*vr(1,n-1)+vr(J+1,n-1))...
+dt*(lambdav*(vl(1,n-1)-vr(1,n-1))+beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))-alpha*vr(1,n-1))-0.25*dt*beta*(vr(1,n-1)+vl(1,n-1))*(ur(2,n-1)-ur(J+1,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vr(1,n-1)-vl(1,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(1,n-1)+vl(1,n-1))*(ul(1,n-1)-ur(1,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))...
*((ur(1,n-1)+ul(1,n-1))-(vr(1,n-1)+vl(1,n-1)))-dt^2*alpha*beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))+0.5*dt^2*alpha^2*vr(1,n-1);
%peridic BC for left moving of infected(j=1)
vl(1,n)=vl(1,n-1)+0.5*mu2*(vl(2,n-1)-vl(J+1,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ul(1,n-1))+0.5*(mu2)^2*(vl(2,n-1)-2*vl(1,n-1)+vl(J+1,n-1))...
+dt*(lambdav*(vr(1,n-1)-vl(1,n-1))+beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))-alpha*vl(1,n-1))+0.25*dt*beta*(vr(1,n-1)+vl(1,n-1))*(ul(2,n-1)-ul(J+1,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vl(1,n-1)-vr(1,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(1,n-1)+vl(1,n-1))*(ur(1,n-1)-ul(1,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))...
*((ur(1,n-1)+ul(1,n-1))-(vr(1,n-1)+vl(1,n-1)))-dt^2*alpha*beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))+0.5*dt^2*alpha^2*vl(1,n-1);
%peridic BC for left moving of infected(j=J+1)
vl(J+1,n)=vl(J+1,n-1)+0.5*mu2*(vl(1,n-1)-vl(J,n-1))*(1-dt*lambdav-dt*alpha+dt*beta*ul(J+1,n-1))+0.5*(mu2)^2*(vl(1,n-1)-2*vl(J+1,n-1)+vl(J,n-1))...
+dt*(lambdav*(vr(J+1,n-1)-vl(J+1,n-1))+beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))-alpha*vl(J+1,n-1))+0.25*dt*beta*(vr(J+1,n-1)+vl(J+1,n-1))*(ul(1,n-1)-ul(J,n-1))*(mu1+mu2)...
+dt^2*lambdav*(vl(J+1,n-1)-vr(J+1,n-1))*(lambdav+alpha)+0.5*dt^2*beta*(vr(J+1,n-1)+vl(J+1,n-1))*(ur(J+1,n-1)-ul(J+1,n-1))*(lambdau+lambdav)+0.5*dt^2*beta^2*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))...
*((ur(J+1,n-1)+ul(J+1,n-1))-(vr(J+1,n-1)+vl(J+1,n-1)))-dt^2*alpha*beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))+0.5*dt^2*alpha^2*vl(J+1,n-1);
end
%calculate the analytical solution
for j=1:J+1
xj=xl+(j-1)*dx;
u_ex(j,n)=exp(-c*(xj-au*t-0.5)^2);
end
%% like this...
plot(x,ur(:,n),'r-','linewidth',2);
title([' time = ' num2str(n)])
axis([0 1 0 1])
drawnow
pause(0.1)
%%
end
  1 Comment
lulu
lulu on 12 May 2021
Hell Rafeal,
I appreciate your answer, but how I used the speed of the infection here (c) with the boundary conditions which I got it.
Thank you

Sign in to comment.

Categories

Find more on Marine and Underwater Vehicles in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!