Plotting Solution for Different Times: Hyperbolic PDE

조회 수: 4 (최근 30일)
Matlab12345
Matlab12345 2020년 4월 19일
댓글: darova 2020년 4월 19일
I am trying to solve the Advection-Diffusion PDE delu/delt - 2*t*delu/delx = 0 using the Maccormack method. The domain is 0 <= x <= 1, t > 0, and I am given the initial conditions u(x,0) = sin(2*pi*x) for 0 <=x<=1 and the periodic boundary condition u(0, t) = u(1,t). The length step is deltax = 0.001, and the time step is deltat = 0.0005. So far, I have the following:
clear;
clc;
%% problem definition and discretization
xstep = 0.001;
tstep = 0.0005;
xdomain= [0 1];
tdomain= [0 1];
nx = round((xdomain(2)-xdomain(1))/xstep);
nt = round((tdomain(2)-tdomain(1))/tstep);
alpha=zeros((nx+1),1);
xt0 = zeros((nx+1),1); % initial condition
xold = zeros((nx+1),1); % solution at timestep k
xnew = zeros((nx+1),1);% solution at timestep k+1
xnew1 = xnew; % intermediate solution
for h = 1:nx+1
alpha(h) =-2*h;
end
%x0 = 0.0; % left boundary condition
%xl = 0.0; % right boundary condition
damping = -0.00; %-0.001;
% initial condition
for i=1:nx+1
xi = (i-1)*xstep;
xt0(i) = 0.0;
if(xi>0 && xi<=1)
xt0(i)=sin(2*3.1416*xi);
end
end
CFL = alpha*tstep/xstep;
xold = xt0;
%xold(1) = x0;
%xold(nx+1) = xl;
for k=1:nt
i = 0;
time = k*tstep;
for i=1:nx+1
% Use periodic boundary condition, u(nx+1)=u(1)
if(i==nx+1)
dudx = xold(1)-xold(nx+1);
du2dx2 = xold(i-1)-2*xold(i)+xold(1);
elseif(i==1)
dudx = xold(i+1)-xold(i);
du2dx2 = xold(nx+1)-2*xold(i)+xold(i+1);
else
dudx = xold(i+1)-xold(i); % C/2*(U(i+1,k)-U(i-1,k)
du2dx2 = xold(i+1)-2*xold(i)+xold(i-1);
end
% Predictor step
xnew1(i) = xold(i) - CFL(i)*dudx + tstep*damping*xold(i); % U(i,k+1)= U(i,k)-C/2*(U(i+1,k)-U(i-1,k)
% correction step
if(i==1)
dudx1 = xnew1(1) - xnew1(nx+1);
else
dudx1 = xnew1(i) - xnew1(i-1);
end
xnew(i) = 0.5*(xold(i) + xnew1(i) - CFL(i)*dudx1) + 0.5*tstep*damping*xnew1(i);
end
% xold= xnew;
end
x=linspace(xdomain(1),xdomain(2),nx+1);
t=linspace(tdomain(1),tdomain(2),(nx+1));
close all
hold on;
h1=plot(x,xt0,'r','linewidth',1);
hold on;
h2=plot(x,xnew,'g','linewidth',1);
hold on;
legend('Exact','Initital','t=0', 't=0.5','t=1')
xlabel('x [m]');
ylabel('Displacement, u(x,t) [m]');
ax=gca;
ax.FontSize= 12;
ax.XLim= [0 1];
ax.YLim= [-20 20];
I am trying to obtain solutions at time = 0, 0.5, and 1. However, I am not sure how to change the value of the time so that the value of the analytical solution changes, and so that the solutions can be stored and plotted. Any help would be appreciated.
  댓글 수: 2
darova
darova 2020년 4월 19일
I only see first derivatives from here
Where did you get second one?
du2dx2 = xold(i-1)-2*xold(i)+xold(1);
Matlab12345
Matlab12345 2020년 4월 19일
편집: Matlab12345 2020년 4월 19일
That line could be used if there was a second derivative, but since there isn't it's presence doesn't seem to make a difference.

댓글을 달려면 로그인하십시오.

답변 (1개)

darova
darova 2020년 4월 19일
편집: darova 2020년 4월 19일
I used this difference scheme
results
Do you have a picture of expected result?

카테고리

Help CenterFile Exchange에서 Eigenvalue Problems에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by