how I can add this condition inside the loop if k= (1-b(g(j)+g(j)) if k>0 then k=(1-b(g(j​)+g(j))+df​(j)))else k=0

5 views (last 30 days)
clear all;
xl=0; xr=1; %domain[xl,xr]
J=200; % J: number of dividion for x
dx=(xr-xl)/J; % dx: mesh size
tf=1000; % final simulation time
Nt=100000; %Nt:number of time steps
dt=tf/Nt;
c=50; % paremeter for the solution
lambdau1=0.05; %turning rate of susp(right)
lambdau2=0.02; %turning rate of susp(left)
beta=0.1; % transimion rate
alpha=0.0; %recovery/death rate
b=110;
d=0.0;
au=0.001; % advection speed - susceptible
A=0.01;
mu1=au*dt/dx;
if mu1>1.0 %make sure dt satisfy stability condition
error('mu1 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
%gr=exp(-c*(x-0.4).^2); %initial conditions of infected
%gl=exp(-c*(x-0.4).^2); %initial conditions of infected
k=6*pi/(xr-xl);
fr=0.0+0.001*(1+sin(k*x)); %dimension f(1:J+1)initial conditions of right moving suscp
fl=0.0+0.001*(1+sin(k*x)); %initial conditions of left moving suscep
ur=zeros(J+1,Nt);
ul=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
%%%%%%%%%%%%%%%%%%%upwind of right-moving of susceptible
ur(j,n)=fr(j)-0.5*mu1*(fr(j)-fr(j-1))+dt*(lambdau2*(1-b*(gr(j)+gl(j)))*(fl(j))-lambdau1*(1-b*(gr(j)+gl(j))+d*fl(j))*(fr(j))-beta*fr(j)*(gr(j)+gl(j))+0.5*alpha*(gr(j)+gl(j)));
%%%%%%%%%%%%%%%%%%%upwind of left-moving of susceptible
ul(j,n)=fl(j)+0.5*mu1*(fl(j+1)-fl(j))+dt*(lambdau1*(1-b*(gr(j)+gl(j)))*(fr(j))-lambdau2*(1-b*(gr(j)+gl(j))+d*fr(j))*(fl(j))-beta*fl(j)*(gr(j)+gl(j))+0.5*alpha*(gr(j)+gl(j)));
%%%%%%%%%%%%%%%%%%%upwind of right-moving of infected
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(J+1)-fr(J))+dt*(lambdau2*(1-b*(gr(J+1)+gl(J+1)))*(fl(J+1))-lambdau1*(1-b*(gr(J+1)+gl(J+1)))*(fr(J+1))^2-beta*fr(J+1)*(gr(J+1)+gl(J+1))+0.5*alpha*(gr(J+1)+gl(J+1)));
%peridic BC for right moving of suscep(j=1,0=J)
ur(1,1)=fr(1)-0.5*mu1*(fr(1)-fr(J))+dt*(lambdau2*(1-b*(gr(1)+gl(1)))*(fl(1))-lambdau1*(1-b*(gr(1)+gl(1)))*(fr(1))-beta*fr(1)*(gr(1)+gl(1))+0.5*alpha*(gr(1)+gl(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(1))+dt*(lambdau1*(1-b*(gr(1)+gl(1)))*(fr(1))-lambdau2*(1-b*(gr(1)+gl(1)))*(fl(1))-beta*fl(1)*(gr(1)+gl(1))+0.5*alpha*(gr(1)+gl(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*(lambdau1*(1-b*(gr(J+1)+gl(J+1)))*(fr(J+1))-lambdau2*(1-b*(gr(J+1)+gl(J+1)))*(fl(J+1))^2-beta*fl(J+1)*(gr(J+1)+gl(J+1))+0.5*alpha*(gr(J+1)+gl(J+1))); %peridic BC for right moving(j=J+1); %peridic BC for left moving
else
for j=2:J % interior nodes
%%%%%%%%%%%%%%%%%%%upwind of right-moving of susceptible
ur(j,n)=ur(j,n-1)-0.5*mu1*(ur(j,n-1)-ur(j-1,n-1))+dt*(lambdau2*(1-b*(vr(j,n-1)+vl(j,n-1)))*(ul(j,n-1))-lambdau1*(1-b*(vr(j,n-1)+vl(j,n-1)))*(ur(j,n-1))-beta*ur(j,n-1)*(vr(j,n-1)+vl(j,n-1))+0.5*alpha*(vr(j,n-1)+vl(j,n-1)));
%%%%%%%%%%%%%%%%%%%upwind of right-moving of susceptible
ul(j,n)=ul(j,n-1)+0.5*mu1*(ul(j+1,n-1)-ul(j,n-1))+dt*(lambdau1*(1-b*(vr(j,n-1)+vl(j,n-1)))*(ur(j,n-1))-lambdau2*(1-b*(vr(j,n-1)+vl(j,n-1)))*(ul(j,n-1))-beta*ul(j,n-1)*(vr(j,n-1)+vl(j,n-1))+0.5*alpha*(vr(j,n-1)+vl(j,n-1)));
end
%peridic BC for right moving of suscep(j=J+1,j+2=1)
ur(J+1,n)=ur(J+1,n-1)-0.5*mu1*(ur(J+1,n-1)-ur(J,n-1))+dt*(lambdau2*(1-b*(vr(J+1,n-1)+vl(J+1,n-1)))*(ul(J+1,n-1))-lambdau1*(1-b*(vr(J+1,n-1)+vl(J+1,n-1)))*(ur(J+1,n-1))-beta*ur(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))+0.5*alpha*(vr(J+1,n-1)+vl(J+1,n-1)));
%peridic BC for right moving of suscep(j=1,0=J)
ur(1,n)=ur(1,n-1)-0.5*mu1*(ur(1,n-1)-ur(J,n-1))+dt*(lambdau2*(1-b*(vr(1,n-1)+vl(1,n-1)))*(ul(1,n-1))-lambdau1*(1-b*(vr(1,n-1)+vl(1,n-1)))*(ur(1,n-1))-beta*ur(1,n-1)*(vr(1,n-1)+vl(1,n-1))+0.5*alpha*(vr(1,n-1)+vl(1,n-1)));
%peridic BC for left moving of suscep (j=1)
ul(1,n)=ul(1,n-1)+0.5*mu1*(ul(2,n-1)-ul(1,n-1))+dt*(lambdau1*(1-b*(vr(1,n-1)+vl(1,n-1)))*(ur(1,n-1))-lambdau2*(1-b*(vr(1,n-1)+vl(1,n-1)))*(ul(1,n-1))-beta*ul(1,n-1)*(vr(1,n-1)+vl(1,n-1))+0.5*alpha*(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+1,n-1))+dt*(lambdau1*(1-b*(vr(J+1,n-1)+vl(J+1,n-1)))*(ur(J+1,n-1))-lambdau2*(1-b*(vr(J+1,n-1)+vl(J+1,n-1)))*(ul(J+1,n-1))^2-beta*ul(J+1,n-1)*(vr(J+1,n-1)+vl(J+1,n-1))+0.5*alpha*(vr(J+1,n-1)+vl(J+1,n-1)));
end
end
figure
ss=surf(T,X,ur+ul);
set(ss,'edgecolor','none');
xlabel('t','FontSize',28);
ylabel('x','FontSize',28);
title('u^++u^-','FontSize',28);
set(gca,'Fontsize',28);
view(2)
  4 Comments

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 14 May 2022
K = max((1-b*(gr(j)+gl(j))+d*fr(j)), 0)
No explicit "if" is needed: max() has conditional logic built in.
  3 Comments
Torsten
Torsten on 14 May 2022
Edited: Torsten on 14 May 2022
In your loops, a variable named "k" is not referenced. So setting k to anything won't influence your results.
Walter Roberson
Walter Roberson on 17 May 2022
You wrote,
"if we assume K=(1-b*(gr(j)+gl(j))+d*fr(j)) then if K>0 then K=(1-b*(gr(j)+gl(j))+d*fr(j)) else K=0 "
There is no existing K in your code, so we do not know where you were intending to insert that, or what you were intending to do with it. But where-ever you were planning to put it, you would insert
K = max((1-b*(gr(j)+gl(j))+d*fr(j)), 0);
Considering the reference to j perhaps you were thinking of inserting it into both of your for j loops.
Your current code does not define gl, gr, or fr (all of those are commented out.)

Sign in to comment.

More Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!