# periodic boundary conditions for advection for a<0

18 views (last 30 days)
lulu on 20 Dec 2020
Commented: lulu on 21 Dec 2020
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');

Alan Stevens on 20 Dec 2020
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 on 21 Dec 2020
Hello, turning rate which allow to turn right and left. Where the flow of right moving> the flow of left moving.

### Categories

Find more on Logical (Boolean) Operations in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!