Numerical Methods for Nonlinear Equations

For the course notes of the Numerical Methods for Nonlinear Equations, please refer to the following links:

• Basic Theory Note, also blog;
• Numerical Methods, also this blog.

In this note, we consider the numerical methods for nonlinear equations, including the finite difference method (FDM), finite volume method (FVM), and finite element method (FEM). Here we take the periodic Burgers’ equation -periodic boundary condition

where are constants, as an example to illustrate the numerical methods.

Finite Volume Method

The finite volume method (FVM) is a numerical technique used for solving partial differential equations (PDEs) that arise in various fields such as fluid dynamics, heat transfer, and electromagnetics. The method is based on the principle of conservation, where the integral form of the governing equations is applied over control volumes.

Godunov’s method

In Godunov’s method, a piecewise constant is defined where for . The is known numerical approximation at time . For a short time interval , the exact solution can be found for as initial values since it is piecewise constant and a sequence of Riemann problems are to be solved. By conservation laws,

Using the relations , and noting that over the cell , we obtain

where the numerical flux function is given by

Because of the fact that the solution of the Riemann problem at is a similarity solution, is constant at the point over the time interval . Moreover, the constant along the line depends only on the data and for this Riemann problem. Denote this value by and the numerical flux function reduces to

and Godunov’s method becomes

For large time interval, the solution may not remain constant at because of the effect of waves arising from neighboring Riemann problems. Since the wave speeds are bounded by the eigenvalues of and the neighboringRiemann problems are distance away, will be constant over provided is sufficiently small. We require that

where .

The Riemann problem with data always has one weak solution consisting simply of the discontinuity propagating with speed . In order to determine for convex equations like Burgers’ equation, there are four cases that must be considered:

1. ,
2. ,
4. where .

The algorithm of solving the periodic Burgers’ equation with Godunov’s method is as follows:

Code: MATLAB code for solving the periodic Burgers' equation with Godunov's method.

Godunov.m

function u = Godunov(m, T)
% solve the periodic Burgers' equation by Godunov's method
% u_t + (0.5 * u^2)_x = 0,  x in [-pi, pi], u(x, 0) = 0.5 + sin(x)
% m is the number of points in space and T is the final time

% initialization
dx = 2*pi / m;
xvec = linspace(-pi, pi, m+1)' - dx/2; xvec = [xvec; pi + dx/2];
u = (0.5 * (diff(xvec)) + (-cos(xvec(2:end)) + cos(xvec(1:end-1))))/dx; % average of integral
tnow = 0;

while tnow < T
    % determine the time step
    dt = 0.5 * dx / (max(abs(u)));  % f'(u) = u for Burgers' equation
    dt = min(dt, T - tnow);  % ensure that the last step is not out of range
    
    % determine the numerical flux
    fhat = nflux_Burgers(u);
    
    % Godunov's method
    u = u - dt / dx * (fhat - [fhat(end-1); fhat(1:end-1)]);  % fhat is periodic
    
    % update
    tnow = tnow + dt;
end
end

function fhat = nflux_Burgers(u)
% compute the numerical flux of Burgers equation(i.e., solve Riemann subproblems)
% The computation is simplified for convex case(f'(u) = u for Burgers equation)

lu = length(u);
fhat = zeros(lu - 1, 1);

for i = 1:lu-1
    % determine the numerical flux by rules of convex function
    if u(i) >= 0
        if u(i+1) >= 0
            % f'(uL) >= 0, f'(uR) >= 0
            fhat(i) = u(i)^2 / 2;
        else
            % f'(uL) > 0, f'(uR) < 0, a shock
            s = (u(i) + u(i+1)) / 2;
            if s > 0
                fhat(i) = u(i)^2 / 2;
            else
                fhat(i) = u(i+1)^2 / 2;
            end
        end
    else
        if u(i+1) <= 0
            % f'(uL) <= 0, f'(uR) <= 0
            fhat(i) = u(i+1)^2 / 2;
        else
            % f'(uL) < 0, f'(uR) > 0, a rarefaction
            fhat(i) = 0;
        end
    end
end
fhat = [fhat; fhat(1)]; % because of periodicity
end

High order Godunov’s method

In high order FV method, i.e., the method of lines, the equation is integrated on and averaged to get

If the solution is smooth enough, the integration and differentiation can be switched:

Thus, we have the semi-discrete numerical scheme

where the numerical flux . We then have an a series of ODE systems

If we could solve a series of above ODE problems numerically, the fully discret numerical solution is gotted.

In order to obtain high order of accuracy, both schemes in space and time should be high order method. Runge-Kutta method is usually chosen to solve ODE problems with high order of accuracy. The 3rd order Runge-Kutta scheme applied here is

As for the space scheme, we need to reconstruct the high order left and right values at the half grids from averages of integration for solving Riemann subproblems. The 3rd order space formula applied here is a central format:

where and .

Once we have the left and right values for Riemann subproblems, we could solve for numerical flux. The Godunov flux is defined as follows for scalar conservation problems

The algorithm for computing 3rd order Godunov flux is as follows:

If is convex, above Godunov flux can be reduced to

1. , ,
2. , ,
3. , ,
4. where , .

In order to get stable results, the time step should not be too large in case that the flux is constant. Thus, we require that

where .

Therefore, the algorithm of solving the periodic Burgers’ equation with Godunov’s method can be summarized as follows:

Code: MATLAB code for solving the periodic Burgers' equation with 3rd Godunov flux.

hoGodunov.m

function u = hoGodunov(m, T)
% solve the periodic Burgers' equation by 3rd order FV scheme
% u_t + (0.5 * u^2)_x = 0,  x in [-pi, pi], u(x, 0) = 0.5 + sin(x) 
% m is the number of points in space and T is the final time

% initialization for cell averages
dx = 2*pi / m;
xvec = linspace(-pi, pi, m+1)' - dx/2; xvec = [xvec; pi + dx/2];
u = (0.5 * (diff(xvec)) + (-cos(xvec(2:end)) + cos(xvec(1:end-1))))/dx;

tnow = 0;
while tnow < T
    % determine the time step
    dt = 0.5 * dx / (max(abs(u)));  % f'(u) = u for Burgers' equation
    dt = min(dt, T - tnow);  % ensure that the last step is not out of range
    
    % 3rd order RK scheme
    fhat = hoGodunov_flux(u);  % 3rd order Godunov flux
    u1 = u - dt / dx * diff(fhat);
    fhat1 = hoGodunov_flux(u1);
    u2 = 0.75 .* u + 0.25 .* (u1 - dt / dx * diff(fhat1));
    fhat2 = hoGodunov_flux(u2);
    u = u ./ 3 + (2/3) .* (u2 - dt / dx * diff(fhat2));
    
    % update for time
    tnow = tnow + dt;
end
end

function fhat = hoGodunov_flux(u)
% compute the high order numerical Godunov  flux for Burgers equation,
% (i.e., solve Riemann subproblems)
% The computation is simplified for convex case(f'(u) = u for Burgers equation)

lu = length(u);
fhat = zeros(lu+1, 1);

% reconstruction by primitive function(k = 2 and S = {I_j-1, I_j, I_j+1})
u = [u(end-2:end-1); u; u(2:3)];
p = zeros(lu+1, 2);
for i = 1:lu+1
    % the left value of half grids
    p(i, 1) = -u(i)/6 + 5*u(i+1)/6 + u(i+2)/3;
    % the right value of half grids
    p(i, 2) = u(i+1)/3 + 5*u(i+2)/6 - u(i+3)/6;
end

for i = 1:lu+1
    % determine the numerical flux by rules of convex function
    if p(i, 1) >= 0
        if p(i, 2) >= 0
            % f'(uL) >= 0, f'(uR) >= 0
            fhat(i) = p(i, 1)^2 / 2;
        else
            % f'(uL) > 0, f'(uR) < 0, a shock
            s = (p(i, 1) + p(i, 2)) / 2;
            if s > 0
                fhat(i) = p(i, 1)^2 / 2;
            else
                fhat(i) = p(i, 2)^2 / 2;
            end
        end
    else
        if p(i, 2) <= 0
            % f'(uL) <= 0, f'(uR) <= 0
            fhat(i) = p(i, 2)^2 / 2;
        else
            % f'(uL) < 0, f'(uR) > 0, a rarefaction
            fhat(i) = 0;
        end
    end
end
end

Lax-Friedrichs flux

In this subsection, we will discuss the 3rd order finite volume scheme with Lax-Friedrichs flux with TVD and TVB limiter.

Recall that the semi-discrete numerical scheme

where the numerical flux . Also, the ODE systems

If we could solve a series of above ODE problems numerically, the fully discret numerical solution is gotted.

In order to obtain high order of accuracy, both schemes in space and time should be high order method. Runge-Kutta method is usually chosen to solve ODE problems with high order of accuracy. The 3rd order Runge-Kutta scheme applied here is

As for the space scheme, we need to reconstruct the high order left and right values at the half grids from averages of integration for solving Riemann subproblems. The 3rd order space formula applied here is a central format:

where and .

Once we have the left and right values for Riemann subproblems, we could solve for numerical flux. The Lax-Friedrichs flux is defined as follows for scalar conservation problems

In order to get stable results, the time step should not be too large in case that the flux is constant. Thus, we require that

where .

In order to eliminate the oscillations in high order finite volume methods, we would modify the left and right values for Riemann subproblems to achieve that goal. This is so called TVD limiter. Recall that the first order FV methods do not suffer from oscillations, the main idea is to combine the high order information into the first order cell averages to obtain high order accuracy while avoiding oscillations.

Denote

and define

where the minmod function is

The modification in are added to the cell averages to get the modified high order FV methods, i.e.,

The left and right values in above equation are substituted into the LF flux to get the numerical fluxes.

With TVD limiter used in and , , however, the order of accuracy would decrease. The remedy is the TVB limiter where TVD is applied if the discrepancies in and are too large. We replace the minmod function above with the new minmod function

where is a prescribed parameter. If , the TVB limiter degenerate to TVD limiter. If is too large, which means no limiter is applied.

Code: solve the periodic Burgers' equation by 3rd order FV scheme using LF flux.

LF3rd.m

function u = LF3rd(m, T)
% solve the periodic Burgers' equation by 3rd order FV scheme
% u_t + (0.5 * u^2)_x = 0,  x in [-pi, pi], u(x, 0) = 0.5 + sin(x) 
% m is the number of points in space and T is the final time

% initialization for cell averages
dx = 2*pi / m;
xvec = linspace(-pi, pi, m+1)' - dx/2; xvec = [xvec; pi + dx/2];
u = (0.5 * (diff(xvec)) + (-cos(xvec(2:end)) + cos(xvec(1:end-1))))/dx;

tnow = 0;
while tnow < T
    % determine the time step
    dt = 0.5 * dx / (max(abs(u)));  % f'(u) = u for Burgers' equation
    dt = min(dt, T - tnow);  % ensure that the last step is not out of range
    
    % 3rd order RK scheme
    fhat = LF3rd_flux(u);  % 3rd order Lax-Friedriches flux
    u1 = u - dt / dx * diff(fhat);
    fhat1 = LF3rd_flux(u1);
    u2 = 0.75 .* u + 0.25 .* (u1 - dt / dx * diff(fhat1));
    fhat2 = LF3rd_flux(u2);
    u = u ./ 3 + (2/3) .* (u2 - dt / dx * diff(fhat2));
    
    % update for time
    tnow = tnow + dt;
end
end

function fhat = LF3rd_flux(u)
% compute the high order numerical Lax-Friedriches flux for Burgers equation,

lu = length(u);

% reconstruction by primitive function(k = 2 and S = {I_j-1, I_j, I_j+1})
up = [u(end-2:end-1); u; u(2:3)];
p = zeros(lu+1, 2);
for i = 1:lu+1
    % the left value of half grids
    p(i, 1) = -up(i)/6 + 5*up(i+1)/6 + up(i+2)/3;
    % the right value of half grids
    p(i, 2) = up(i+1)/3 + 5*up(i+2)/6 - up(i+3)/6;
end

% alpha = max(|f'(u)|) where f'(u) = u
alpha = max(abs(u));
% fhat = [f(a) + f(b) - alpha * (b-a)]/2
fhat = (sum(p .^ 2, 2) ./ 2 - alpha .* diff(p, 1, 2)) ./ 2;
end

Code: solve the periodic Burgers' equation by 3rd order FV scheme using LF flux with TVD limiter.

LF3rd_TVD.m

function u = LF3rd_TVD(m, T)
% solve the periodic Burgers' equation by 3rd order FV scheme with TVD
% limiter. m is the number of points in space and T is the final time. 
% the equation is u_t + (0.5 * u^2)_x = 0,  x in [-pi, pi], u(x, 0) = 0.5 + sin(x) 

% initialization for cell averages
dx = 2*pi / m;
xvec = linspace(-pi, pi, m+1)' - dx/2; xvec = [xvec; pi + dx/2];
u = (0.5 * (diff(xvec)) + (-cos(xvec(2:end)) + cos(xvec(1:end-1))))/dx;

tnow = 0;
while tnow < T
    % determine the time step
    dt = 0.5 * dx / (max(abs(u)));  % f'(u) = u for Burgers' equation
    dt = min(dt, T - tnow);  % ensure that the last step is not out of range
    
    % 3rd order RK scheme
    fhat = LF3rd_flux_TVD(u);  % 3rd order Lax-Friedriches flux with TVD limiter
    u1 = u - dt / dx * diff(fhat);
    fhat1 = LF3rd_flux_TVD(u1);
    u2 = 0.75 .* u + 0.25 .* (u1 - dt / dx * diff(fhat1));
    fhat2 = LF3rd_flux_TVD(u2);
    u = u ./ 3 + (2/3) .* (u2 - dt / dx * diff(fhat2));
    
    % update for time
    tnow = tnow + dt;
end
end

function fhat = LF3rd_flux_TVD(u)
% compute the high order numerical Lax-Friedriches flux for Burgers equation

lu = length(u);
% reconstruction by primitive function(k = 2 and S = {I_j-1, I_j, I_j+1})
up = [u(end-2:end-1); u; u(2:3)];
p = zeros(lu+1, 2);
for i = 1:lu+1
    % the left value of half grids
    p(i, 1) = -up(i)/6 + 5*up(i+1)/6 + up(i+2)/3;
    % the right value of half grids
    p(i, 2) = up(i+1)/3 + 5*up(i+2)/6 - up(i+3)/6;
end

% TVD limiter
u_t = p(:, 1) - up(2:lu+2);
u_tt = up(3:lu+3) - p(:, 2);
diffup = diff(up);
umod_t = zeros(lu+1, 1);
umod_tt = zeros(lu+1, 1);
for i = 1:lu+1
    umod_t(i) = minmod(u_t(i), diffup(i), diffup(i+1));
    umod_tt(i) = minmod(u_tt(i), diffup(i+1), diffup(i+2));
end
% correction
p = [up(2:lu+2) up(3:lu+3)] + [umod_t -umod_tt];


% alpha = max(|f'(u)|) where f'(u) = u
alpha = max(abs(u));
% fhat = [f(a) + f(b) - alpha * (b-a)]/2
fhat = (sum(p .^ 2, 2) ./ 2 - alpha .* diff(p, 1, 2)) ./ 2;
end

Code: solve the periodic Burgers' equation by 3rd order FV scheme using LF flux with TVB limiter.

LF3rd_TVB.m

function u = LF3rd_TVB(m, T, M)
% solve the periodic Burgers' equation by 3rd order FV scheme using LF flux
% with TVB limiter. m is the number of points in space, T is the final time
% and M is TVB parameter. The equation is u_t + (0.5 * u^2)_x = 0,  x in
% [-pi, pi], u(x, 0) = 0.5 + sin(x)

% initialization for cell averages
dx = 2*pi / m;
xvec = linspace(-pi, pi, m+1)' - dx/2; xvec = [xvec; pi + dx/2];
u = (0.5 * (diff(xvec)) + (-cos(xvec(2:end)) + cos(xvec(1:end-1))))/dx;

tnow = 0;
while tnow < T
    % determine the time step
    dt = 0.5 * dx / (max(abs(u)));  % f'(u) = u for Burgers' equation
    dt = min(dt, T - tnow);  % ensure that the last step is not out of range
    
    % 3rd order RK scheme
    fhat = LF3rd_flux_TVB(u, M * dx^2);  % 3rd order Lax-Friedriches flux with TVB limiter
    u1 = u - dt / dx * diff(fhat);
    fhat1 = LF3rd_flux_TVB(u1, M * dx^2);
    u2 = 0.75 .* u + 0.25 .* (u1 - dt / dx * diff(fhat1));
    fhat2 = LF3rd_flux_TVB(u2, M * dx^2);
    u = u ./ 3 + (2/3) .* (u2 - dt / dx * diff(fhat2));
    
    % update for time
    tnow = tnow + dt;
end
end

function fhat = LF3rd_flux_TVB(u, Mbound)
% compute the high order numerical Lax-Friedriches flux for Burgers equation

lu = length(u);
% reconstruction by primitive function(k = 2 and S = {I_j-1, I_j, I_j+1})
up = [u(end-2:end-1); u; u(2:3)];
p = zeros(lu+1, 2);
for i = 1:lu+1
    % the left value of half grids
    p(i, 1) = -up(i)/6 + 5*up(i+1)/6 + up(i+2)/3;
    % the right value of half grids
    p(i, 2) = up(i+1)/3 + 5*up(i+2)/6 - up(i+3)/6;
end

% TVB limiter
u_t = p(:, 1) - up(2:lu+2);
u_tt = up(3:lu+3) - p(:, 2);
diffup = diff(up);
umod_t = zeros(lu+1, 1);
umod_tt = zeros(lu+1, 1);
for i = 1:lu+1
    if abs(u_t(i)) < Mbound
        umod_t(i) = u_t(i);
    else
        umod_t(i) = minmod(u_t(i), diffup(i), diffup(i+1));
    end
    if abs(u_tt(i)) < Mbound
        umod_tt(i) = u_tt(i);
    else
        umod_tt(i) = minmod(u_tt(i), diffup(i+1), diffup(i+2));
    end
end
% correction
p = [up(2:lu+2) up(3:lu+3)] + [umod_t -umod_tt];


% alpha = max(|f'(u)|) where f'(u) = u
alpha = max(abs(u));
% fhat = [f(a) + f(b) - alpha * (b-a)]/2
fhat = (sum(p .^ 2, 2) ./ 2 - alpha .* diff(p, 1, 2)) ./ 2;
end

Code: the minmod function used in the above code.

minmod.m

function m = minmod(a, b, c)
% the minmod function returns the minimum magnitude of elements if
% they have the same sign and 0 otherwise
vec = [a; b; c];
signs = sign(vec);
same_sign = all(signs == signs(1));
if same_sign
    if signs(1) == 1
        m = min(vec);
    else
        m = max(vec);
    end
else
    m = 0;
end
end

WENO

Recall that the semi-discrete numerical scheme

where the numerical flux . Also, the ODE systems

If we could solve a series of above ODE problems numerically, the fully discret numerical solution is gotted.

As for the space discretization, we need to reconstruct the high order left and right values at the half grids from averages of integration for solving Riemann subproblems. There are different stencils and we would like to choose the one in which less oscillations are observed. At the same time, we would like to take advantage of all the information provided by the stencils we have compared. To be precise, if the small stencils we used are

then the reconstruction from the union of these stencils is

where and are the linear weights.

The main idea of WENO scheme is to change the linear weights to nonlinear weights s.t. and several conditions are satisfied:

1. , . (This implies consistency.)
2. In smooth case, , . (The high order of accuracy is maintained.)
3. When has a discontinuity in one or more stencils, we would hope the corresponding weights to be essentially 0. (No oscillations is introduced.)

The third statement requires that smoothness is determined during the computation. The smooth indicator is introduced as follows

For the special case , we have

Therefore, the nonlinear weights can be defined as

These weights can be approximated separately in smooth and nonsmooth case:

The weighted reconstruction is

We see that is small if is not smooth in and the stencils related to that is less involved in the computation. To obtain the procedure above should be modified symmetrically w.r.t . The flux applied here is Lax-Friedriches flux, i.e.,

In order to obtain high order of accuracy, the time discretization should be high order method. Runge-Kutta method is usually chosen to solve ODE problems with high order of accuracy. The 3rd order Runge-Kutta scheme applied here is

In order to get stable results, the time step should not be too large and satisfy the CFL condition. Thus, we require that

where . In WENO scheme, the space order of accuracy is 5 if the underlying solution is smooth and we requires

where is from CFL condition. This makes the order of time scheme match the order of space scheme.

Code: solve the periodic Burgers' equation by WENO scheme

WENO.m

function u = WENO(m, T)
% solve the periodic Burgers' equation by WENO scheme
% u_t + (0.5 * u^2)_x = 0,  x in [-pi, pi], u(x, 0) = 0.5 + sin(x) 
% m is the number of points in space and T is the final time

% initialization for cell averages
dx = 2*pi / m;
xvec = linspace(-pi, pi, m+1)' - dx/2; xvec = [xvec; pi + dx/2];
u = (0.5 * (diff(xvec)) + (-cos(xvec(2:end)) + cos(xvec(1:end-1))))/dx;

tnow = 0;
texp = 5/3; % in accordance with order in space
while tnow < T
    % determine the time step
    dt = 0.5 * dx / (max(abs(u)));  % f'(u) = u
    dt = dt ^ texp;
    dt = min(dt, T - tnow);  % last step is treated here
    
    % 3rd order RK scheme
    fhat = WENO_LF(u);  % 5th order WENO scheme
    u1 = u - dt / dx * diff(fhat);
    fhat1 = WENO_LF(u1);
    u2 = 0.75 .* u + 0.25 .* (u1 - dt / dx * diff(fhat1));
    fhat2 = WENO_LF(u2);
    u = u ./ 3 + (2/3) .* (u2 - dt / dx * diff(fhat2));
    
    % update for time
    tnow = tnow + dt;
end
end

function fhat = WENO_LF(u)
% compute the high order numerical Lax-Friedriches flux using WENO scheme

lu = length(u);
epsilon = 1e-6;

% reconstruction by primitive function(k = 2 and S_0, S_1, S_2)
U = [u(end-3:end-1); u; u(2:4)]; % extended vector
uvecs = {U(1:lu+1); U(2:lu+2); U(3:lu+3); U(4:lu+4); U(5:lu+5); U(6:lu+6)};

% the left values(u-_{j+1/2})
ul = [(1/3) .* uvecs{1} + (-7/6) .* uvecs{2} + (11/6) .* uvecs{3}, ...
    (-1/6) .* uvecs{2} + (5/6) .* uvecs{3} + (1/3) .* uvecs{4}, ...
    (1/3) .* uvecs{3} + (5/6) .* uvecs{4} + (-1/6) .* uvecs{5}];

% the right values(u+_{j-1/2})
ur = [ul(:, 2), ul(:, 3), (11/6) .* uvecs{4} + (-7/6) .* uvecs{5} + (1/3) .* uvecs{6}];

% linear weights
dl = [1/10, 3/5, 3/10];
dr = fliplr(dl);

% smooth indicators
betal = [(13/12) .* (uvecs{1} - 2 .* uvecs{2} + uvecs{3}) .^ 2 + ...
    (1/4) .* (uvecs{1} - 4 .* uvecs{2} + 3 .* uvecs{3}) .^ 2, ...
    (13/12) .* (uvecs{2} - 2 .* uvecs{3} + uvecs{4}) .^ 2 + ...
    (1/4) .* (uvecs{2} - uvecs{4}) .^ 2, ...
    (13/12) .* (uvecs{3} - 2 .* uvecs{4} + uvecs{5}) .^ 2 + ...
    (1/4) .* (3 .* uvecs{3} - 4 .* uvecs{4} + uvecs{5}) .^ 2];
betar = [(13/12) .* (uvecs{2} - 2 .* uvecs{3} + uvecs{4}) .^ 2 + ...
    (1/4) .* (uvecs{2} - 4 .* uvecs{3} + 3 .* uvecs{4}) .^ 2, ...
    (13/12) .* (uvecs{3} - 2 .* uvecs{4} + uvecs{5}) .^ 2 + ...
    (1/4) .* (uvecs{3} - uvecs{5}) .^ 2, ...
    (13/12) .* (uvecs{4} - 2 .* uvecs{5} + uvecs{6}) .^ 2 + ...
    (1/4) .* (3 .* uvecs{4} - 4 .* uvecs{5} + uvecs{6}) .^ 2];

% nonlinear weights
alphal = dl ./ (betal + epsilon) .^ 2;
alphar = dr ./ (betar + epsilon) .^ 2;
wl = alphal ./ sum(alphal, 2);
wr = alphar ./ sum(alphar, 2);

% flux from reconbination
ul = sum(wl .* ul, 2);
ur = sum(wr .* ur, 2);

% LF flux = [f(a) + f(b) - alpha * (b-a)]/2
fhat = ((ul .^ 2 + ur .^ 2) ./ 2 - max(abs(u)) .* (ur - ul)) ./ 2;
end

Finite Difference Method

The finite difference method (FDM) is the most classical numerical technique for solving differential equations by approximating derivatives with finite differences. The basic idea is to replace the continuous derivatives in the differential equation with discrete approximations using values at grid points. Here we still restrict our discussion to the Burgers’ equation.

Recall that the semi-discrete numerical scheme for finite difference method is

If we want to get high order FV scheme, the derivative should be approximated by

where is the numerical flux, s.t.

Thus a high order FD scheme can be obtained. Suppose there is a function , s.t.

Then,

and we could take in order to satisfy the high order approximation requirement . The reconstruction process in finite volume method can be used to fulfill that purpose. Given the values on grids , we have

Then, should be reconstructed to get and we simply take . If we could solve a series of ODE problems numerically, the fully discret numerical solution is gotted.

Since there would be two s (namely the values on the left and right end of and , respectively), we have to choose which numerical flux should be used. We apply flux splitting method which splits with

and have as many derivatives as the order of the scheme. The value on the right end of using and that on the left end of using are added up to get the final numerical flux. The Lax-Friedriches splitting is used:

Algorithm for flux splitting

As for the WENO scheme, the space reconstruction on small stencils are obtained

and the reconstruction from the union of these stencils is

where and are the linear weights.

The main idea of WENO scheme is to change the linear weights to nonlinear weights s.t. and several conditions are satisfied:

1. , . (This implies consistency.)
2. In smooth case, , . (The high order of accuracy is maintained.)
3. When has a discontinuity in one or more stencils, we would hope the corresponding weights to be essentially 0. (No oscillations is introduced.)

The third statement requires that smoothness should be determined during the computation. The smooth indicator is introduced as follows

For the special case , we have

Therefore, the nonlinear weights can be defined as

These weights can be approximated separately in smooth and nonsmooth case:

The weighted reconstruction is

We see that is small if is not smooth in and the stencils related to that is less involved in the computation. To obtain the procedure above should be modified symmetrically w.r.t .

In order to obtain high order of accuracy, the time discretization should be high order method. Runge-Kutta method is usually chosen to solve ODE problems with high order of accuracy. The 3rd order Runge-Kutta scheme applied here is

In order to get stable results, the time step should not be too large and satisfy the CFL condition. Thus, we require that

where . In WENO scheme, the space order of accuracy is 5 if the underlying solution is smooth and we requires

where is from . This makes the order of time scheme match the order of space scheme.

Main Algorithm

Algorithm for WENO scheme

Code: solve the periodic Burgers' equation by fifth FD-WENO scheme.

FD_WENO.m

function u = FD_WENO(m, T)
% solve the periodic Burgers' equation by fifth FD WENO scheme
% u_t + (0.5 * u^2)_x = 0,  x in [-pi, pi], u(x, 0) = 0.5 + sin(x) 
% m is the number of points in space and T is the final time

% initialization for values of grids
dx = 2*pi / m;
xvec = linspace(-pi, pi, m+1)';
u = 0.5 + sin(xvec);

tnow = 0;
texp = 5/3; % in accordance with order in space
while tnow < T
    % determine the time step
    dt = 0.5 * dx / (max(abs(u)));  % f'(u) = u
    dt = dt ^ texp;
    dt = min(dt, T - tnow);  % last step is treated here
    
    % 3rd order RK scheme
    fhat = FD_WENO_LF(u);  % 5th order WENO scheme
    u1 = u - dt / dx * diff(fhat);
    fhat1 = FD_WENO_LF(u1);
    u2 = 0.75 .* u + 0.25 .* (u1 - dt / dx * diff(fhat1));
    fhat2 = FD_WENO_LF(u2);
    u = u ./ 3 + (2/3) .* (u2 - dt / dx * diff(fhat2));
    
    % update for time
    tnow = tnow + dt;
end
end

function fhat = FD_WENO_LF(u)
% compute the numerical flux for FD scheme using Lax-Friedriches splitting
% and WENO scheme

% Lax-Friedriches splitting
alpha = max(abs(u));  % max|f'(u)|
fplus = ((u .^ 2 ) ./ 2 + alpha .* u) ./ 2;
fminus = ((u .^ 2 ) ./ 2 - alpha .* u) ./ 2;

fhat = WENOreconstruction(fplus, fminus);
end

function fhat = WENOreconstruction(uplus, uminus)
% reconstruction by WENO scheme. Only values of right ends of subintervals
% of uplus and values of left ends of subintervals of uminus are taken

lu = length(uplus);
epsilon = 1e-6;

% reconstruction by primitive function(k = 2 and S_0, S_1, S_2)
Up = [uplus(end-3:end-1); uplus; uplus(2:3)]; % extended vector
Um = [uminus(end-2:end-1); uminus; uminus(2:4)];
upvecs = {Up(1:lu+1); Up(2:lu+2); Up(3:lu+3); Up(4:lu+4); Up(5:lu+5)};
umvecs = {Um(1:lu+1); Um(2:lu+2); Um(3:lu+3); Um(4:lu+4); Um(5:lu+5)};

% the left values(u-_{j+1/2})
ul = [(1/3) .* upvecs{1} + (-7/6) .* upvecs{2} + (11/6) .* upvecs{3}, ...
    (-1/6) .* upvecs{2} + (5/6) .* upvecs{3} + (1/3) .* upvecs{4}, ...
    (1/3) .* upvecs{3} + (5/6) .* upvecs{4} + (-1/6) .* upvecs{5}];

% the right values(u+_{j-1/2})
ur = [(-1/6) .* umvecs{1} + (5/6) .* umvecs{2} + (1/3) .* umvecs{3},...
    (1/3) .* umvecs{2} + (5/6) .* umvecs{3} + (-1/6) .* umvecs{4},...
    (11/6) .* umvecs{3} + (-7/6) .* umvecs{4} + (1/3) .* umvecs{5}];

% linear weights
dl = [1/10, 3/5, 3/10];
dr = fliplr(dl);

% smooth indicators
betal = smooth_indicators(upvecs);
betar = smooth_indicators(umvecs);

% nonlinear weights
alphal = dl ./ (betal + epsilon) .^ 2;
alphar = dr ./ (betar + epsilon) .^ 2;
wl = alphal ./ sum(alphal, 2);
wr = alphar ./ sum(alphar, 2);

% flux from reconbination
ul = sum(wl .* ul, 2);
ur = sum(wr .* ur, 2);

% average
fhat = ul + ur;
end

function beta = smooth_indicators(u)
% smooth indicator for k = 2
beta = [(13/12) .* (u{1} - 2 .* u{2} + u{3}) .^ 2 + ...
    (1/4) .* (u{1} - 4 .* u{2} + 3 .* u{3}) .^ 2, ...
    (13/12) .* (u{2} - 2 .* u{3} + u{4}) .^ 2 + ...
    (1/4) .* (u{2} - u{4}) .^ 2, ...
    (13/12) .* (u{3} - 2 .* u{4} + u{5}) .^ 2 + ...
    (1/4) .* (3 .* u{3} - 4 .* u{4} + u{5}) .^ 2];
end

Finite Element Method

The finite element method (FEM) is a numerical technique for finding approximate solutions to boundary value problems for partial differential equations. The basic idea of FEM is to divide the domain into smaller subdomains (elements) and use piecewise polynomial functions to approximate the solution within each element. The global solution is then obtained by assembling the local solutions from all elements.

Here we restrict our discussion to the discontinuous Galerkin (DG) method, which is a specific type of FEM that allows for discontinuities between elements. The DG method is particularly useful for hyperbolic conservation laws, such as the Burgers’ equation.

In DG method, the basis functions are similar with those of FE method except that there are no constraints upon the continuity. Therefore, the space of basis functions are

where , are subintervals over the whole domain and denotes polynomials whose degree are no more than .

The weak form of source equation is gotted by multiplying it with a test function and integrating on

The target of DG method is to find an approximation function in s.t. the weak form are satisfied for any . To be specific, that is

where the flux is replaced by the numerical flux and the limit values of in are taken since it may be discontinuous. The numerical flux can be taken as monotone flux or approximate Riemann solver. Here we choose the local Lax-Friedriches flux where , i.e., the maximum velocity over the subintervals adjacent to .

More details of DG method can be shown here. It is well known that the dimension of is and a common choice of basis is monomial . However, this basis may lead to ill-conditioned linear system and result in loss of accuracy. A better choice is Legendre polynomial which are mutually orthogonal in and , . The shifted Legendre polynomials in are

where . It’s easy to prove that the polynomials in are mutually orthogonal in and , . The first three Legendre polynomials are easy to get as following

To get the approximation of initial condition , we can do projection onto . That is, we seek s.t.

If is satisfied for every basis function in , it would be satisfied by all . Thus, we could take as the basis of and solve for the equations formed. Once we expand into in the subinterval , we get

due to the orthogonality of shift Legendre polynomials. If we compute the rhs of by high order Gauss quadrature (in order to avoid destroying the converge order of the final scheme), i.e., , the full system is

where is a diagonal matrix, and .

The same idea for obtaining initial approximation can be applied to solve . If we have in , then

As for the second term in , we compute it by high order numerical quadrature scheme, i.e.

Note that

Hence

where and are weights and points of Gauss quadrature.

As long as the second term in can be calculated accurate enough (without hurting the accuracy of final scheme), we would have

where is again a diagonal matrix, is the derivative of coefficients with respect to time, , and .

There are numerous methods to solve the ODEs and a third order Runge-Kutta method is used here.

where and is the numerical solution at time . The restrict on time step should be placed by CFL condition

where is the order in polynomial space .

Algorithm of DG method

Code: solve the periodic Burgers' equation by DG scheme.

DG.m

function u = DG(m, T)
% solve the periodic Burgers' equation by DG scheme
% u_t + (0.5 * u^2)_x = 0,  x in [-pi, pi], u(x, 0) = 0.5 + sin(x) 
% m is the number of points in space and T is the final time

% initialization for coefficients in subintervals
dx = 2*pi / m;
u = initialcoeffs(@(x) 0.5 + sin(x), -pi, pi, m);

tnow = 0;
while tnow < T
    % determine the time step
    dt = 0.2 * dx / (coeffs2maxabs(u));  % f'(u) = u
    dt = min(dt, T - tnow);  % last step is treated here
    
    % 3rd order RK scheme
    rhs = coeffs2rhs(u);  % compute the rhs
    u1 = u + dt * rhs;
    rhs1 = coeffs2rhs(u1);
    u2 = 0.75 .* u + 0.25 .* (u1 + dt * rhs1);
    rhs2 = coeffs2rhs(u2);
    u = u ./ 3 + (2/3) .* (u2 + dt * rhs2);
    
    % update for time
    tnow = tnow + dt;
end
end

% ------------ coeffs of f in m equaispaced subintervals ----------------%
function u = initialcoeffs(f, a, b, m)
% get the coeffs of f in m equaispaced subintervals of [a, b]
u = zeros(m+1, 3);  % P^2 polynomial
dx = (b - a) / m;
xvec = linspace(a, b, m+1)';

% projection on Legendre polynomials
xleg = [-sqrt(15)/5; 0; sqrt(15)/5];
alphaleg = [5, 8, 5] ./ 9;
scale = [1, 3, 5] ./ 2;
P0 = [1; 1; 1];  % P0 = 1
P1 = xleg;  % P1 = x
P2 = (3 .* xleg .^2 - 1) ./ 2;  % P2 = (3 * x^2 - 1) / 2
for i = 1:m+1
    fgauss = f(xleg .* (dx/2) + xvec(i));
    % 3 points Gauss quadrature
    a0 = alphaleg * (fgauss .* P0);
    a1 = alphaleg * (fgauss .* P1);
    a2 = alphaleg * (fgauss .* P2);
    
    % scale (P^i are orthogonal but not orthonormal)
    u(i, :) = [a0, a1, a2] .* scale;
end
end

% ----- compute the maximum of absolute of u over all subintervals ------%
function alpha = coeffs2maxabs(u)

% Note that the extreme point of a Legendre polynomial(k <= 2) is -a_1/(3 * a_2)

% compute the values at end points first
endvals = u * [1, 1; -1, 1; 1, 1];
alpha = max(abs(endvals), [], 'all');
% compare the extreme point
expoint = (u(:, 2) .* (-1/3)) ./ u(:, 3);
m = size(u, 1);
for i = 1:m
    if expoint(i) > -1 && expoint(i) < 1
        alpha = max(alpha, abs(u(i, :) * [1; expoint(i); (3 * expoint(i)^2 - 1)/2]));
    end
end
end

% --------------------- compute the rhs in DG method ---------------------%
function rhs = coeffs2rhs(u)

m = size(u, 1);
dx = 2 * pi / (m-1);

% compute the numerical LF flux
% compute the values at half grids
ur = u * [1; -1; 1];  % right values of half grids
ul = u * [1; 1; 1];   % left values of half grids
ur = [ur; ur(2)]; ul = [ul(end-1); ul];
% LF flux = [f(a) + f(b) - alpha * (b-a)]/2
fhat = ((ul .^ 2 + ur .^ 2) ./ 2 - maxabs(u) .* (ur - ul)) ./ 2;

% compute the values at 3 point Gauss points
gaussvals = [1,                1,               1; 
                     -sqrt(15)/5, 0, sqrt(15)/5;
                     2/5,      -1/2,           2/5];
ugauss = u * gaussvals;
fugauss = ugauss .^ 2 ./ 2;  %f(u) = u^2/2

% compute the integral
% (P^0)_x = 0, (P^1)_x = 1, (P^2)_x = 3*x
alphaleg = [5; 8; 5] ./ 9;
P1x = (fugauss .* [1, 1, 1]) * alphaleg;
P2x = (fugauss .* (3 .* gaussvals(2, :))) * alphaleg;

% add up
rhs = [zeros(m, 1) - fhat(2:end) + fhat(1:end-1),...
          P1x - fhat(2:end) + fhat(1:end-1) .* (-1),...
          P2x - fhat(2:end) + fhat(1:end-1)] .* ([1, 3, 5] ./ dx);
end

% ---- maximum of absolute of u over every adjacent two subintervals ---- %
function um = maxabs(u)
% compute the maximum of absolute of u over every adjacent two subintervals
% Note that the extreme point of a Legendre polynomial(k <= 2) is -a_1/(3 * a_2)

% compute the values at end points first
endvals = u * [1, 1; -1, 1; 1, 1];
maxu = max(abs(endvals), [], 2);
% compare the extreme point
expoint = (u(:, 2) .* (-1/3)) ./ u(:, 3);
m = size(u, 1);
for i = 1:m
    if expoint(i) > -1 && expoint(i) < 1
        maxu(i) = max(maxu(i), abs(u(i, :) * [1; expoint(i); (3 * expoint(i)^2 - 1)/2]));
    end
end
um = zeros(m+1, 1);
um(1) = max(maxu(1), maxu(end-1));
um(end) = um(1);
for i = 2:m-1
    um(i) = max(maxu(i), maxu(i+1));
end
end