numerical simulation of travelling wave of epidemic model

조회 수: 9 (최근 30일)
lulu
lulu 2021년 5월 10일
댓글: lulu 2021년 5월 12일
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);

이 질문에 답변하려면 로그인하십시오.

답변 (1개)

Rafael Hernandez-Walls
Rafael Hernandez-Walls 2021년 5월 11일
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
lulu
lulu 2021년 5월 12일
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

댓글을 달려면 로그인하십시오.

카테고리

Help Center 및 File Exchange에서 Marine and Underwater Vehicles에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by