Theory of Nonlinear Equations

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

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

Introduction

To show the complexity of the nonlinear equations, we first consider the classical shock tube problem. The structure of this flow turns out to be very interesting, involving three distinct waves separating regions in which the state variables are constant.

Shock Tube Problem
The physical set-up of this problem is a tube filled with gas, initially divided by a membrane into two sections. The gas has a higher density and pressure in one half of the tube than in the other half, with zero velocity everywhere. At time , the membrane is suddenly removed or broken, and the gas allowed to flow. We expect a net motion in the direction of lower pressure. Assuming the flow is uniform across the tube, there is variation in only one direction and the one-dimensional Euler equations apply, aka

where is the density, is the velocity, is the total energy, and is the pressure. The equations are coupled through the equation of state, which relates the pressure to the density and internal energy. And the boundary is given by the initial conditions:

where , and problem with this kind boundary condition is also called Riemann problem.

Shock Tube Problem

The solution of this problem with initial condition

after some time is displayed in the following figure. It’s observed that the solution can be divided into five regions, separated by four points :

• Region I: , where the solution is constant and equal to .
• Region II: , where the solution is a rarefaction wave, which is a smooth transition from state in region I to state in region II.
• Region III: , where the solution is some constant.
• At point , there is a contact discontinuity, since the density is discontinuous but the velocity and pressure are continuous.
• Region IV: , where the solution is some constant.
• At point , there is a shock wave, since the density, velocity, and pressure are all discontinuous.
• Region V: , where the solution is constant and equal to .

Solution to a shock tube problem for the one-dimensional Euler equations

In mathematics, the difficulties of such problem include:

1. Even if the initial data is smooth, the solution may develop discontinuities in finite time, which results in the loss of regularity of the solution, i.e, the strong solution may not exist.
2. If we look for weak solutions, the uniqueness of the solution is not guaranteed, and we need to impose some additional conditions on the weak solutions.

Correspondingly, a good numerical method should satisfy the following two properties:

1. Shock tracking: track the discontinuities in the solution.
2. Shock capturing: sharp resolution of the discontinuities without excesive smearing or spurious oscillations.

Nonlinear Hyperbolic Conservation Laws

The nonlinear hyperbolic conservation laws are a class of partial differential equations that describe the evolution of a system over time. They are typically written in the form:

where is the vector value function of conserved variables and , is the flux function. The system is hyperbolic if the eigenvalues of the Jacobian matrix of the flux function are real and distinct. We first consider the one-dimensional case, where the system is deduced to a scalar conservation law:

The characteristic curves of this system are given by the equation:

where is the characteristic speed, and on this curve we have

therefore, the solution is constant along the characteristic curves, i.e.

Burger’s Equation
The most simple example of nonlinear hyperbolic conservation law is the Burger’s equation, which is given by

where is the velocity of the fluid. The flux function is given by , and the characteristic speed is given by . The solution to this equation can develop shock waves, which are discontinuities in the solution, even with smooth initial data.

For short time, the solution of Burger’s equation can be still smooth, which is shown in the following figure. Note that the characteristic curves of Burger’s equation are all straight lines, since , while we proved that the solution is constant along the characteristic curves, hence . In fact, the solution is given by the following equation:

where is the characteristic curve passing through the points and .

Characteristics and solution for Burgers' equation (small )

However, as time goes on, the solution transforms in different speeds, regions with higher value moves faster than regions with lower value, and the solution becomes steeper and steeper. Eventually at some time where the characteristics first cross, the function has an infinite slope – the wave “breaks” and a shock forms. Beyond this point there is no classical solution of the PDE, and the weak solution we wish to determine becomes discontinuous. In fact, we can show that if the initial data is smooth, then the wave will break at

Shock formation in Burgers' equation.

For times some of the characteristics have crossed and so there are points where there are three characteristics leading back to . One can view the “solution” at such a time as a triple-valued function.

Triple-valued solution to Burgers' equation at time .

Remark. Actually the Burger’s equation given above should be called the “inviscid Burgers’ equation”, since the equation studied by Burgers also includes a viscous term:

and the solution of inviscid Burgers’ equation is only a approximation of the solution of viscous Burgers’ equation if is small.

Sovle Burger’s Equation by Characteristics Method

In this subsection, we consider using the characteristics method to solve the periodic Burger’s equation by characteristic method at 3rd order Gauss points. The problem is given by

with periodic boundary condition and initial condition . Take , for example.

First do a variable substitution:

We have

and

We know that satisfies the equation

Thus, we can solve for to get .

Recall for Burgers’ equation, the characteristics are

Since the function values are the same along characteristics, we only have to determine the right initial value to get the value of . That is

For solving above nonlinear equation, we could apply Newton’s method. The initial iterate of Newton’s method is taken as the value of the first characteristic on the left of among the characteristics emitted from 100 equispaced points in . The algorithm and corresponding MATLAB code are given below.

Code: MATLAB code for solving the nonlinear equation by Newton's method.

nonlinearsolve.m

function u = nonlinearsolve(x, t)
% solve the nonlinear equation u = sin(x - u*t) by Newton's method

% convert to a reference interval
x = mod(x + pi, 2*pi) - pi;

% determine the initial value
xvec = linspace(-pi, pi, 100).';
uvec = sin(xvec);
u = uvec(find(xvec + uvec .* t >= x, 1));

for i = 1:50
    fu = u - sin(x - u * t); % f(u)
    if abs(fu) < 1e-13
        % solution found
        break
    end

    dfu = 1 + cos(x - u * t) * t; % f'(u)
    du = - dfu \ fu; % Newton update
    u = u + du;
end
end

Now the main algorithm and corresponding MATLAB code can be given as follows:

Code: MATLAB code for solving the periodic Burger's equation by characteristic method.

character.m

function uvec = character(n, t)
% solve the periodic Burger's equation by characteristic method at 3rd
% order Gauss points

% constants
alpha = 0.5; beta = 1;

% xvec = linspace(-pi, pi, n+1) + [-sqrt(15)/5; 0; sqrt(15)/5] .* (pi/n);
% xvec = xvec(:);
xvec = linspace(-pi, pi, n+1)';
uvec = zeros(length(xvec), 1);
for i = 1:length(xvec)
    xi = xvec(i) - alpha * t;
    tildeui = nonlinearsolve(xi, beta * t);
    uvec(i) = alpha + beta * tildeui;
end

And the numerical solutions with at and are given in the following figures. The solution at is smooth, while the solution at has a shock wave.

Numerical solutions of Burger's equation within the interval at .

Weak Solution

A natural way to define a generalized solution of the inviscid equation that does not require differentiability is to go back to the integral form of the conservation law. Following the analysis in note of real analysis, we define the weak solution via the integral form of the conservation law and the integration by parts. Note that if is a solution of the conservation law, then it satisfies the integral form of the conservation law:

for any and , where is any test function with enough regularity. We can integrate by parts to get

Let , and , , then we have

for all , by using for all . Especially, if we take

then we also have

for all and .

It’s easy to prove that and are equivalent, and we call the function a weak solution of the conservation law if it satisfies for all , or for all and .

Weak Solution of Burger’s equation
Now we consider the weak solution of Burger’s equation. To simplify the analysis, we restrict ourselves to the Riemann problem of Burger’s equation, which is given by the piecewise constant initial condition:

the form of the solution depends on the relation between and .

Case 1: .
In this case there is a unique weak solution, i.e. a shock wave, which is a discontinuous solution of the form

where

is the shock speed, the speed at which the discontinuity travels. Note that characteristics in each of the regions where is constant go into the shock (see the following figure) as time advances.

Shock wave

Case 2: .
In this case there are infinitely many weak solutions. One of these is again the solution in case 1, in which the discontinuity propagates with speed . Note that characteristics now go out of the shock (see the following figure) and that this solution is not stable to perturbations . If the data is smeared out slightly, or if a small amount of viscosity is added to the equation, the solution changes completely.

Entropy-violating shock

Another weak solution is the rarefaction wave, which is a smooth transition from to , given by

This solution is stable to perturbations and is in fact the vanishing viscosity generalized solution (See the following figure), that is, is the limit of the solution of the viscous equation

as . We can determine the correct physical behavior by adopting the vanishing viscosity approach, since the Burger’s equation is a model of above viscous equation valid only for small and smooth . When it breaks down, we must return to original viscous equation. If the initial data is smooth and very small, then before the wave begins to break the term is negligible compared to the other terms and the solutions to both PDEs look nearly identical. The solution would be essentially unchanged if we solved viscous equation with small rather than the Burger’s equation. However, as the wave begins to break, the second derivative term grows much faster than , and at some point the eurr term is comparable to the other terms and begins to play a role. This term keeps the solution smooth for all time, preventing the breakdown of solutions that occurs for the hyperbolic problem.

Rarefaction wave

Shock Speed

The shock speed is defined to be

where is the curve along which the solution is discontinuous, i.e. the shock front. The shock speed is the speed at which the discontinuity propagates through the medium.

The propagating shock solution is a weak solution to Burgers’ equation only if the speed of propagation is given by . The same discontinuity propagating at a different speed would not be a weak solution.

In general, the speed of propagation can be determined by conservation: First, a curve of discontinuity (shock front) is inserted, separating two regions in the plane where the solution is smooth. These discontinuities are moving internal boundaries. The basic rule governing their propagation may be obtained in an informal way, by appealing to the fundamental principle of conservation. We draw a small box of dimension by in the plane, small enough so that the solution to either side of the discontinuity can be thought of as approximately constant, denoted by , and the speed of the discontinuity is approximately constant, denoted by , with the box dimensions chosen so that if and are small enough. Requiring that the total ‘mass’ at the new time equals the old mass plus net mass flow through the boundaries gives

which is an approximation of the integral form of the conservation law

then we have

assuming , this gives the shock speed as

which is the Rankine-Hugoniot condition. This condition is a necessary condition for the shock to be a weak solution of the conservation law.

Shock relations arise naturally by integrating the solution over a small control volume.

The equal area rule. One technique that is useful for determining weak solutions by hand is to start with the “solution” constructed using characteristics (which may be multi-valued if characteristics cross) and then eliminate the multi-valued parts by inserting shocks. The shock location can be determined by the “equal area rule”, which is best understood by looking at the following figure. The shock is located such that the shaded regions cut off on either side have equal areas. This is a consequence of conservation – the integral of the discontinuous weak solution must be the same as the area “under” the multi-valued solution, since both “solve” the same conservation law.

Equal area rule for shock location

Entropy Solution

The vanishing viscosity generalized solution we defined before (for Burger’s equation) is a weak solution in the sense of and , and so this definition includes the solution we are looking for. Unfortunately, weak solutions are often not unique, and so an additional problem is often to identify which weak solution is the physically correct vanishing viscosity solution. Again, one would like to avoid working with the viscous equation directly, but it turns out that there are other conditions one can impose on weak solutions that are easier to check and will also pick out the correct solution, which are usually called entropy conditions by analogy with the gas dynamics case. The vanishing viscosity solution is also called the entropy solution because of this.

For scalar equations there is an obvious condition suggested by the figure presented before. A shock should have characteristics going into the shock, as time advances. A propagating discontinuity with characteristics coming out of it, is unstable to perturbations. Either smearing out the initial profile a little, or adding some viscosity to the system, will cause this to be replaced by a rarefaction fan of characteristics. This gives our first version of the entropy condition:

Entropy condition (Version 1).
A discontinuity propagating with speeds given by RH condition satisfies the entropy condition if

Note that is the characteristic speed. For convex , the Rankine-Hugoniot speed must lie between and , so reduces to simply the requirement that f’(ul) > f’(u,), which again by convexity requires .

A more general form of this condition, due to Oleinik, applies also to nonconvex scalar flux functions :

Entropy condition (Version 2).
is the entropy solution if all discontinuities have the property that

For convex , this requirement reduces to the first version. Another form of the entropy condition is based on the spreading of characteristics in a rarefaction fan. If is an increasing function of in some region, then characteristics spread out if . The rate of spreading can be quantified, and gives the following condition, also due to Oleinik.

Entropy condition (Version 3).
is the entropy solution if there is a constant such that for all , and ,

Entropy Function

Yet another approach to the entropy condition is to define an entropy function for which an additional conservation law holds for smooth solutions that becomes an inequality for discontinuous solutions. In gas dynamics, there exists a physical quantity called the entropy that is known to be constant along particle paths in smooth flow and to jump to a higher value as the gas crosses a shock. It can never jump to a lower value, and this gives the physical entropy condition that picks out the correct weak solution in gas dynamics.

Suppose some function satisfies a conservation law of the form

for some entropy flux . Then we can obtain from this, for smooth ,

Since is the solution of the conservation law, we have , and so we can write the above equation as

For a scalar conservation law this equation admits many solutions . For a system of equations and are still scalar functions, but now the equation above reads, which is a system of equations for the two variables and . If , this may have no solutions.

An additional condition we place on the entropy function is that it be convex, , for reasons that will be seen below.

The entropy is conserved for smooth flows by its definition. For discontinuous solutions, however, the manipulations performed above are not valid. Since we are particularly interested in how the entropy behaves for the vanishing viscosity weak solution, we look at the related viscous problem and will then let the viscosity tend to zero. The viscous equation is

Since solutions to this equation are always smooth, we can derive the corresponding evolution equation for the entropy following the same manipulations we used for smooth solutions of the inviscid equation, multiplying above equation by to obtain

We can now rewrite the right hand side to obtain

Integrating this equation over the rectangle gives

As , the first term on the right hand side vanishes. (This is clearly true if is smooth at and , and can be shown more generally.) The other term, however, involves integrating . over the If the limiting weak solution is discontinuous along a curve in this rectangle, then this term will not vanish in the limit. However, since , and (by our convexity assumption), we can conclude that the right hand side is nonpositive in the limit and hence the vanishing viscosity weak solution satisfies

for all and . Alternatively, in integral form,

Consequently, the total integral of is not necessarily conserved, but can only decrease. (Note that our mathematical assumption of convexity leads to an “entropy function” that decreases, whereas the physical entropy in gas dynamics increases.) The fact that above inequality holds for all and is summarized by saying that in the weak sense. This gives our final form of the entropy condition, called the entropy inequality.

Entropy Condition (Version 4).
The function is the entropy solution of if, for all convex entropy functions and corresponding entropy fluxes, the inequality

is satisfied in the weak sense.

This formulation is also useful in analyzing numerical methods. If a discrete form of this entropy inequality is known to hold for some numerical method, then it can be shown that the method converges to the entropy solution.

Just as for the conservation law, an alternative weak form of the entropy condition can be formulated by integrating against smooth test functions 0, now required to be nonnegative since the entropy condition involves an inequality. The weak form of the entropy inequality is