periodic boundary conditions for advection for a<0

12 次查看（过去 30 天）
lulu
lulu 2020-12-20
评论： lulu 2020-12-21
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');

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

回答（1 个）

Alan Stevens
Alan Stevens 2020-12-20
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.
  5 个评论
显示 3更早的评论
Alan Stevens
Alan Stevens 2020-12-20
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
lulu 2020-12-21
Hello, turning rate which allow to turn right and left. Where the flow of right moving> the flow of left moving.

请先登录，再进行评论。

类别

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

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by