Solve 1D Burgers Equation via Semi-Implicit Adams-Bash​forth/Cran​k-Nicolson Scheme

Hello world,

I'm trying to solve the 1D Burgers equation (nonlinear convection-diffusion equation) by applying the explicit Adams-Bashforth scheme to the nonlinear convective term and the implicit Crank-Nicolson scheme to the diffusive term. The equation in dimensionless form is

where the second term is nonlinear. Dirichlet boundary conditions on are , while those Neumann are .

The exact solution is

so that the initial condition and steady state solutions respectively are

and

Note that there is a discontinuity in the exact solution at . To avoid this, we can start the solution at where is a small parameter on the order of or . Thus, the initial condition becomes

I wish to construct a solution using 100 spatial steps in so that , and 278 time steps in so that . This makes the Courant number (convective stability) where since is the maximum velocity at found by a perturbation of . Note that there is no restriction on viscous stability here since the plan is to use an implicit method on the diffusion term.

I have solved the 1D diffusion equation using the implicit Crank-Nicolson method, I have solved the 1D Burgers equation using the explicit Euler method, and I have solved several ODEs using the explicit two-step Adams-Bashforth method. I think I understand the theory of this problem, and on paper I have written out the time discretized equations, but my skill in MATLAB is novice and I have no idea how to write code that combines these methods to solve a single PDE. Here, since I'm trying to use the Crank-Nicolson method on the diffusive term in the Burgers equation, I suspect that I can just apply my code from the diffusion equation (in conjunction with some unwritten code for the Adams-Bashforth method) to these initial and boundary conditions.

The MATLAB code from my Crank-Nicolson solution to the 1D diffusion equation applied directly to the initial and boundary conditions of the 1D Burgers equation is:

clear all; close all; clc;

% Construct spatial mesh

Nx = 100;           % Number of spatial steps
xl = -6;            % Left x boundary
xr = 6;             % Right x boundary
dx = (xr-xl)/Nx;    % Spatial step
x = xl:dx:xr;       % x

% Construct temporal mesh

tf = 2;             % Final time
dt = 0.0072;        % Timestep
Nt = round(tf/dt);  % Number of timesteps
t = 0:dt:tf;        % t

% Parameters

umax = 15;                      % Found by a perturbation of t=10^-2
C = umax*dt/dx;                 % Convective Stability / Courant Number
disp("Courant Number: "+C)
V = dt/(dx*dx);                 % Viscous Stability / Diffusion Number
disp("Diffusion Number: "+V)

% Define Boundary & Initial Conditions

% Boundary conditions (Dirichlet)

for j = 1:Nt+1                 % Loop through all t
    u(1,j) = 2;                % u(-6,t) = 2
    u(Nx+1,j) = -2;            % u(6,t) = -2 
end

% Initial condition (t = t0, NOT t = 0)

for i = 1:Nx+1                % Loop through all x                 
    u(i,1) = (-2*sinh(x(i)))/(cosh(x(i))-exp(0.0001));   % u(x,t0)
end

% Define the right (MMr) and left (MMl) matrices in the Crank-Nicolson method

aal(1:Nx-2) = -V;                   % Below the main diagonal in MMl
bbl(1:Nx-1) = 2+2*V;                % The main diagonal in MMl
ccl(1:Nx-2) = -V;                   % Above the main diagonal in MMl
MMl = diag(bbl,0)+diag(aal,-1)+diag(ccl,1);     % Building the MMl matrix

aar(1:Nx-2) = V;                    % Below the main diagonal in MMr
bbr(1:Nx-1) = 2-2*V;                % The main diagonal in MMr
ccr(1:Nx-2) = V;                    % Above the main diagonal in MMr
MMr = diag(bbr,0)+diag(aar,-1)+diag(ccr,1);     % Building the MMr matrix

% Implementation of the Crank-Nicholson method

for j = 1:Nt                        
    u(2:Nx,j+1) = inv(MMl)*MMr*u(2:Nx,j);   
end

Plotting this solution:

Clearly this is problematic, but hopefully a good start. In addition to incorporating the Adams-Bashforth scheme, I think I also need to incorporate the Neumann boundary conditions. If anyone out there is familiar with this semi-implicit AB/CN method, I would be very appreciative of some guidance. I know this is an elaborate questions, so please don't hesitate to write comments and communicate with me about anything.

Thank you!

Answer 1

[Mat Hunt]

Make sure your vaiue of D*dt/dx^2 isn't too large, ideally, it should be around 1.