Project 4: Burger's Equation

by Matt Herman, Omar Safsaf, John Wilson

Background

In this project, we demonstrate the Implicit Newton solver by using the Backward Difference Equation with Newton iteration to solve Burgers' equation, defined as $u_t + uu_x = Du_{xx}$. First, in Exercise 8.4.3, we show that a function $u(x,t)$ is a solution of the Burgers' equation with Dirichlet boundary conditions. Next, in the two Computer Problems, we demonstrate usage of the solver and generate plots of it using various parameters.

Exercise 8.4.3

Computer Problem 8.4.1

Solve Burgers' equation on $[0, 1]$ with intial condition $f(x) = \sin 2\pi x$ and boundary conditions $l(t) = r(t) = 0$, using step sizes (a) $h = k = 0.1$ and (b) $h = k = 0.02$. Plot the approximate solutions for $0 \leq t \leq 1$. Which equilibrium solution does the solution approach as time increases?

Original burgers.m

1234567891011121314151617181920212223242526272829303132% Program 8.7 Implicit Newton solver for Burgers equation 
% input: space interval [xl,xr], time interval [tb,te], 
%        number of space steps M, number of time steps N 
% output: solution w 
% Example usage: w=burgers(0,1,0,2,20,40) 
function w=burgers(xl,xr,tb,te,M,N) 
alf=5;bet=4;D=.05; 
f=@(x) 2*D*bet*pi*sin(pi*x)./(alf+bet*cos(pi*x)); 
l=@(t) 0*t; 
r=@(t) 0*t; 
h=(xr-xl)/M; k=(te-tb)/N; m=M+1; n=N; 
sigma=D*k/(h*h); 
w(:,1)=f(xl+(0:M)*h)'; % initial conditions
w1=w; 
for j=1:n 
  for it=1:3 % Newton iteration 
    DF1=zeros(m,m);DF2=zeros(m,m); 
    DF1=diag(1+2*sigma*ones(m,1))+diag(-sigma*ones(m-1,1),1); 
    DF1=DF1+diag(-sigma*ones(m-1,1),-1); 
    DF2=diag([0;k*w1(3:m)/(2*h);0])-diag([0;k*w1(1:(m-2))/(2*h);0]); 
    DF2=DF2+diag([0;k*w1(2:m-1)/(2*h)],1)-diag([k*w1(2:m-1)/(2*h);0],-1); 
    DF=DF1+DF2; 
    F=-w(:,j)+(DF1+DF2/2)*w1; % Using Lemma 8.11 
    DF(1,:)=[1 zeros(1,m-1)]; % Dirichlet conditions for DF 
    DF(m,:)=[zeros(1,m-1) 1]; 
    F(1)=w1(1)-l(j);F(m)=w1(m)-r(j); % Dirichlet conditions for F 
    w1=w1-DF\F; 
  end 
  w(:,j+1)=w1; 
end 
x=xl+(0:M)*h;t=tb+(0:n)*k; 
mesh(x,t,w')

Minor Changes

123456789% Change this function:
f=@(x) 2*D*bet*pi*sin(pi*x)./(alf+bet*cos(pi*x)); 
% to this:
f=@(x) sin(2*pi*x)

% Add axis labels:
xlabel("x")
ylabel("t")
zlabel("u")

And Behold!!

w = burgers(0,1,0,1,10,10)

Space interval: $[0,1]$

Time interval: $[0,1]$

Number of space steps: 10

Number of time steps: 10

w = burgers(0,1,0,1,50,50)

Space interval: $[0,1]$

Time interval: $[0,1]$

Number of space steps: 50

Number of time steps: 50

w = burgers(0,1,0,5,50,250)

Space interval: $[0,1]$

Time interval: $[0,5]$

Number of space steps: 250

Number of time steps: 250

Summary

In each of these plots, it's clear that the solution approaches the equilibrium solution $u = 0$ as time increases. The last two plots, not requested in the original question, were generated to make this explicitly clear using a wider time interval and a much higher step count, as well as being captured at two different angles.

Computer Problem 8.4.2

Computer Problem 2 requires solving Burger's Equation using the implicit Newton solver and Homogeneous Dirichlet Boundary conditions. Alpha, beta, and delta were changed in Matlab code 8.7. M and N were varied to 100 and 32 respectively to account for the required step sizes. The code is:

123456789101112131415161718192021222324252627function w=burgers(xl,xr,tb,te,M,N) 
alf=4;bet=3;D=.2; 
f=@(x) 2*D*bet*pi*sin(pi*x)./(alf+bet*cos(pi*x)); 
l=@(t) 0*t; 
r=@(t) 0*t; 
h=(xr-xl)/M; k=(te-tb)/N; m=M+1; n=N; 
sigma=D*k/(h*h); 
w(:,1)=f(xl+(0:M)*h)'; % initial conditions
w1=w; 
for j=1:n 
  for it=1:3 % Newton iteration 
    DF1=zeros(m,m);DF2=zeros(m,m); 
    DF1=diag(1+2*sigma*ones(m,1))+diag(-sigma*ones(m-1,1),1); 
    DF1=DF1+diag(-sigma*ones(m-1,1),-1); 
    DF2=diag([0;k*w1(3:m)/(2*h);0])-diag([0;k*w1(1:(m-2))/(2*h);0]); 
    DF2=DF2+diag([0;k*w1(2:m-1)/(2*h)],1)-diag([k*w1(2:m-1)/(2*h);0],-1); 
    DF=DF1+DF2; 
    F=-w(:,j)+(DF1+DF2/2)*w1; % Using Lemma 8.11 
    DF(1,:)=[1 zeros(1,m-1)]; % Dirichlet conditions for DF 
    DF(m,:)=[zeros(1,m-1) 1]; 
    F(1)=w1(1)-l(j);F(m)=w1(m)-r(j); % Dirichlet conditions for F 
    w1=w1-DF\F; 
  end 
  w(:,j+1)=w1; 
end 
x=xl+(0:M)*h;t=tb+(0:n)*k; 
mesh(x,t,w')

The Plot of the solution is:

The second part of the problem requires finding an error. The error is found by subtracting the implicit solver solution from the solution given by the equation:

After finding the solution, a function was written to compute the error as a function of k step size. The code is:

123456789function s = errorp2 (n)
  u=.130921              %using the equation...
  errarray={}
  for j=4:n %running burgers equation multiple times
    z=burgers(0,1,0,2,100,(2*2^(j)))
    err= z(50,(2^(j)))-u
    errarray=[errarray,err]
    
  end

The output of the function is an array of errors. This will be plotted as a function of K. The plots:

Note that the slope of the logarithmic function is around 1. This shows that there was an expected first-order decrease in error.