numerical simulation of travelling wave of epidemic model

3 Ansichten (letzte 30 Tage)
lulu
lulu am 10 Mai 2021
Kommentiert: lulu am 12 Mai 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);

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (1)

Rafael Hernandez-Walls
Rafael Hernandez-Walls am 11 Mai 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 Kommentar
lulu
lulu am 12 Mai 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

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Digital Filter Analysis finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by