Write a program to compute the unsteady behavior of the following equation =v in

Question

Write a program to compute the unsteady behavior of the following equation =v in a periodic domain of length l with v= 0.01. This is usually called the nonlinear advection-diffusion equation or the Burgers equation Take the initial conditions to be f(x,t = 0)-sin(2m)+ 1.0 and approximate the equation using a forward in time approximation for the time derivative, and centered approximations for the first and second derivative. Follow the evolution up to time 1.0 using at least three different grid resolutions, the finest of which should use at least 200 grid points. in class

Explanation / Answer

% one-dimensional advection-diffusion by the FTCS scheme
n=21; nstep=100; length=2.0; h=length/(n-1); dt=0.05; D=0.05;
f=zeros(n,1); y=zeros(n,1); ex=zeros(n,1); time=0.0;
for i=1:n, f(i)=0.5*sin(2*pi*h*(i-1)); end; % initial conditions
for m=1:nstep, m, time
for i=1:n, ex(i)=exp(-4*pi*pi*D*time)*...
0.5*sin(2*pi*(h*(i-1)-time)); end; % exact solution
hold off; plot(f,'linewidt',2); axis([1 n -2.0, 2.0]); % plot solution
hold on; plot(ex,'r','linewidt',2);pause; % plot exact solution
y=f; % store the solution
for i=2:n-1,
f(i)=y(i)-0.5*(dt/h)*(y(i+1)-y(i-1))+...
D*(dt/h^2)*(y(i+1)-2*y(i)+y(i-1)); % advect by centered differences
end;
f(n)=y(n)-0.5*(dt/h)*(y(2)-y(n-1))+...
D*(dt/h^2)*(y(2)-2*y(n)+y(n-1)); % do endpoints for
f(1)=f(n); % periodic boundaries
time=time+dt;
end;

% MATLAB program to produce publication quality plots.
% The data is in file example.data with the first column containing time. Here we plot the
% data in the second to fourth column versus time.
%--------------------------------------------------------------
% Load and plot the data
data=load('example.data');
d1=plot(data(:,1),data(:,2)); set(d1,'LineWidth',2,'LineStyle','-','Color', 'k')
hold on
d2= plot(data(:,1),data(:,3)); set(d2,'LineWidth',2,'LineStyle','--','Color', 'k')
d3= plot(data(:,1),data(:,4)); set(d3,'LineWidth',2,'LineStyle','-.','Color', 'k')
% Set the identity of each line
l1=0.5;h1=0.25; c1=plot([l1 l1+0.4],[h1 h1],'-'); set(c1,'LineWidth',2,'LineStyle','-','Color',
'k'); text(l1+0.5,h1,'data1','Fontsize',18)
l2=0.5;h2=0.75; c2=plot([l2 l2+0.4],[h2 h2],'--'); set(c2,'LineWidth',2,'LineStyle','--','Color',
'k'); text(l2+0.5, h2,'data2','Fontsize',18)
l3=0.5;h3=1.25; c3=plot([l3 l3+0.4],[h3 h3],'-.'); set(c3,'LineWidth',2,'LineStyle','-.','Color',
'k'); text(l3+0.5,h3,'data3','Fontsize',18)
% Set the appearance of the plot
axis([0 5 0 5])
xlabel('Time','Fontsize',18)
ylabel('Reynolds number','Fontsize',18)
set(gca,'Box','on'); set(gca,'Fontsize',18, 'LineWidth',2)
hold off
% to make a eps file of the figure, use: print -deps exampleplot

Related Questions

Navigate

• Browse (All)
• Browse W
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.