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 vues (au cours des 30 derniers jours)
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 commentaires
lulu
lulu le 13 Mai 2022
hyperbolic type. This is bilogical model.

Connectez-vous pour commenter.

Réponse acceptée

Walter Roberson
Walter Roberson le 14 Mai 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 commentaires
Torsten
Torsten le 14 Mai 2022
Modifié(e) : Torsten le 14 Mai 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 le 17 Mai 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.)

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Loops and Conditional Statements dans Help Center et File Exchange

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by