numerical simulation of travelling wave of epidemic model

9 visualizzazioni (ultimi 30 giorni)
lulu
lulu il 10 Mag 2021
Commentato: lulu il 12 Mag 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);

Accedi per rispondere a questa domanda.

Risposte (1)

Rafael Hernandez-Walls
Rafael Hernandez-Walls il 11 Mag 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 Commento
lulu
lulu il 12 Mag 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

Accedi per commentare.

Categorie

Scopri di più su Marine and Underwater Vehicles in Help Center e File Exchange

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by