# 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

1 view (last 30 days)
lulu on 13 May 2022
Commented: Walter Roberson on 17 May 2022
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)
lulu on 13 May 2022
hyperbolic type. This is bilogical model.

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.
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.)

### Categories

Find more on Eigenvalue Problems in Help Center and File Exchange

R2021b

### Community Treasure Hunt

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

Start Hunting!