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

4 次查看（过去 30 天）
lulu
lulu 2022-5-13
评论： Walter Roberson 2022-5-17
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 个评论
显示 2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2022-5-14
K = max((1-b*(gr(j)+gl(j))+d*fr(j)), 0)
No explicit "if" is needed: max() has conditional logic built in.
  3 个评论
显示 1更早的评论
Torsten
Torsten 2022-5-14
编辑：Torsten 2022-5-14
In your loops, a variable named "k" is not referenced. So setting k to anything won't influence your results.
Walter Roberson
Walter Roberson 2022-5-17
@ima Kindy
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.)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 General PDEs 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by