How I write periodic BCs

2 views (last 30 days)
lulu
lulu on 1 Dec 2020
Commented: lulu on 1 Dec 2020
clear all;
xl=0; xr=1; %domain[xl,xr]
J=40; % J: number of dividion for x
dx=(xr-xl)/J; % dx: mesh size
tf=0.5; % final simulation time
Nt=50; %Nt:number of time steps
dt=tf/Nt;
c=50; % paremeter for the solution
mu=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.2).^2); %dimension f(1:J+1)
%store the solution at all grid points for all time steps
u=zeros(J+1,Nt);
%find the approximate solution at each time step
for n=1:Nt
t=n*dt; % current time
gl = exp(-c*(xl-t-0.2)^2); % BC at left side
gr = exp(-c*(xr-t-0.2)^2); % BC at right side
if n==1 % first time step
for j=2:J % interior nodes
u(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));
end
u(1,n) = gl; % the left-end point
u(J+1,n) = gr; % the right-end point
else
for j=2:J % interior nodes
u(j,n)=u(j,n-1)-0.5*mu*(u(j+1,n-1)-u(j-1,n-1))+0.5*mu^2*(u(j+1,n-1)-2*u(j,n-1)+u(j-1,n-1));
end
u(1,n)=gl; % the left-end point
u(J+1,n)=gr; % the right-end point
end
%calculate the analytical solution
for j=1:J+1
xj=xl+(j-1)*dx;
u_ex(j,n)=exp(-c*(xj-t-0.2)^2);
end
end
%plot the analytical and numerical solution at different times
figure;
hold on;
n=10;
plot(x,u(:,n),'r.',x, u_ex(:,n),'r-');
n=30;
plot(x,u(:,n),'g.',x, u_ex(:,n),'g-');
n=50;
plot(x,u(:,n),'b.',x, u_ex(:,n),'b-');
legend('Numerical t=0.1','Analytical t=0.1','Numerical t=0.3','Analytical t=0.3','Numerical t=0.5','Analytical t=0.5');
title('Numerical and Analytical Solutions at t=0.1,0.3,0.5');
  4 Comments
lulu
lulu on 1 Dec 2020
No, f is the initial condition
I want to change gl and gr to periodic boundary conditions. If the initial conditions leaves the domain at one end, and re-enters it at the other end then the Periodic BCs are correct.

Sign in to comment.

Answers (0)

Tags

Products

Community Treasure Hunt

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

Start Hunting!

Translated by