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

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)

[Jan]

It is not clear, where you want to include this condition and what "put zero" means. Please explain the problem with more details. Remove all parts from the code, which have no relation to the problem.

[Torsten]

@ima Kindy comment moved here

Hello thank you to contact me, I want to understand how I write loop or condition inside for n=1:Nt and j=2:J:

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

[Torsten]

@ima Kindy

What PDE are you trying to solve ?

[lulu]

hyperbolic type. This is bilogical model.

Answer 1

[Walter Roberson]

K = max((1-b*(gr(j)+gl(j))+d*fr(j)), 0)

No explicit "if" is needed: max() has conditional logic built in.

[lulu]

sorry, but where I have to keep

K = max((1-b*(gr(j)+gl(j))+d*fr(j)), 0).

[Torsten]

In your loops, a variable named "k" is not referenced. So setting k to anything won't influence your results.

[Walter Roberson]

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