periodic boundary conditions for advection for a<0

The following shows the codes for advection of population move to right and left . I got the correct BCs for right moving that I am able to see the solution profile leaving the domain at one end, and re-entering at the other end. Whereas, I am struggled to get a correct BCs for left moving. Could you please help.

clear all;

% numerical simulation for ur_t+a*ur_x=0 and ul_t-a*ul_x=0 where ur: population moving to right, ul: population moving to left

xl=0; xr=1; %domain[xl,xr]

J=1000; % J: number of dividion for x

dx=(xr-xl)/J; % dx: mesh size

tf=50; % final simulation time

Nt=10000; %Nt:number of time steps

dt=tf/Nt;

c=50; % paremeter for the solution

lambdaR=0.005; %turning rate

lambdaL=0.003; %turning rate

a=0.1; % speed

mu=a*dt/dx;

if mu>1.0 %make sure dt satisfy stability condition

error('mu should<1.0!')

end

%Evaluate the initial conditions

x=xl:dx:xr; %generate the grid point

f=exp(-c*(x-0.5).^2); %dimension f(1:J+1)

%store the solution at all grid points for all time steps

ur=zeros(J+1,Nt);

ul=zeros(J+1,Nt);

%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

% u(j,n)=f(j)- mu *(f(j) - f(j-1)); %upwind scheme

ur(j,n)=f(j)-0.5*mu*(f(j+1)-f(j-1))+0.5*mu^2*(f(j+1)-2*f(j)+f(j-1))+lambdaR*f(j)-lambdaL*f(j); %lax wendroff

ul(j,n)=f(j)+0.5*mu*(f(j+1)-f(j-1))+0.5*mu^2*(f(j+1)-2*f(j)+f(j-1))-lambdaR*f(j)+lambdaL*f(j);

end

ur(J+1,1)=ur(J,1); %peridic BC for right moving

ur(1,1)=ur(J+1,1); %peridic BC for right moving

ul(J+1,1)=ul(1,1); %peridic BC for left moving

ul(J,1)=ul(J-1,1); %peridic BC for left moving

else

for j=2:J % interior nodes

% u(j,n)=u(j,n-1)- mu *(u(j,n) - u(j-1,n)); %upwind scheme

ur(j,n)=ur(j,n-1)-0.5*mu*(ur(j+1,n-1)-ur(j-1,n-1))+0.5*mu^2*(ur(j+1,n-1)-2*ur(j,n-1)+ur(j-1,n-1))+lambdaR*ul(j,n-1)-lambdaL*ur(j,n-1);

ul(j,n)=ul(j,n-1)+0.5*mu*(ul(j+1,n-1)-ul(j-1,n-1))+0.5*mu^2*(ul(j+1,n-1)-2*ul(j,n-1)+ul(j-1,n-1))-lambdaR*ul(j,n-1)+lambdaL*ur(j,n-1);

end

ur(J+1,n)=ur(J,n); %peridic BC for right moving

ur(1,n)=ur(J+1,n); %peridic BC for right moving

ul(J+1,n)=ul(1,n); %peridic BC for left moving

ul(J,n)=ul(J-1,n); %peridic BC for left moving

end

%calculate the analytical solution

for j=1:J+1

xj=xl+(j-1)*dx;

u_ex(j,n)=exp(-c*(xj-a*t-0.5)^2);

end

end

%plot the analytical and numerical solution at different times

figure;

hold on;

n=1;

plot(x,ur(:,n),'r-');

plot(x,ul(:,n),'r--');

%plot(x,ur(:,n),'r-',x,ul(:,n),'k-',x ,u_ex(:,n),'r.'); %red

n=1000;

plot(x,ur(:,n),'g-');

plot(x,ul(:,n),'g--');

%plot(x,ur(:,n),'g-',x,ul(:,n),'b-',x, u_ex(:,n),'g.'); %green

n=1600;

plot(x,ur(:,n),'b-'); %blue

plot(x,ul(:,n),'b--');

n=2000;

plot(x,ur(:,n),'m-');

plot(x,ul(:,n),'m--');

n=3000;

plot(x,ur(:,n),'k-'); %black

plot(x,ul(:,n),'k--'); %black

% %%% time=n*\deta t=

legend('Numericalrght-moving t=0','Numericalleft-moving t=0','Numericalrght-moving t=5','Numericalleft-moving t=5','Numericalrght-moving t=8','Numericalleft-moving t=8','Numericalrght-moving t=10','Numericalleft-moving t=10','Numericalrght-moving t=15','Numericalleft-moving t=15');% title('Numerical and Analytical Solutions at t=0.1,0.3,0.5');

Answer 1

[Alan Stevens]

For left moving instead of

ul(J,n)=ul(J-1,n);  %peridic BC for left moving

shouldn't you have

ul(J,n)=ul(J+1,n);  %peridic BC for left moving

assuming J+1 is to the right of J.

[lulu]

Thank you for your responce , I changed but the profile not correct. I has leave fron one end and re-enter from the other end e.g for green curve at x=0 , the solution is 0.15 , but at x=1 the solution is 0.

see the attached.

[Alan Stevens]

Is there anything that means the left moving flow would be different from the right moving flow (apart from direction, of course!). If not, then you could simply do the right moving flow calculations, then at the end set the left moving flow as

ul(J+1:-1:1,:) = ur(1:J+1,:);

[lulu]

Thank you, could you please explain to me how it working with that set. The flow of right moving with turning rate lambdaR>lambdaL(turning rate of left moving)

[Alan Stevens]

I don't know what turning rate means here. However, if lambdaR being different from lambdaL means flow to the left would not be the same as flow to the right then do the left-flow calculations exactly as for the right-flow ones in the first instance, except with lambdaL and lambdaR interchanged. Call this UL. Then at the end do

ul(J+1:-1:1,:) = UL(1:J+1,:);

[lulu]

Hello, turning rate which allow to turn right and left. Where the flow of right moving> the flow of left moving.