Piecewise Self-Similar Solutions and a Numerical Scheme for Scalar Conservation Laws.
The solution of the Riemann problem was a building block for general Cauchy problems in conservation laws. A Cauchy problem is approximated by a series of Riemann problems in many numerical schemes.  But, since the structure of the Riemann solution holds locally in time only, and, furthermore, a Riemann solution is not piecewise constant in general, there are several fundamental issues in this approach such as the stability and the complexity of computation.In this article we introduce a new approach which is based on piecewise self-similar solutions.  The scheme proposed in this article solves the problem without the time marching process.  The approximation error enters in the step for the initial discretization only, which is given as a similarity summation of base functions.  The complexity of the scheme is linear. Convergence to the entropy solution and the error estimate are shown.  The mechanism of the scheme is introduced in detail together with several interesting properties of the scheme.
Introduction
. Self-similarity of the Cauchy problem for one dimensional conservation
laws,
with Riemann initial data
has been the basis of various schemes devised for general initial value problems, Glimm
[9] and Godunov [10], for example. The self-similarity of the Riemann problem is the
property that the solution is a function of the self-similarity variable x=t. In other
words the solution is constant along the self-similarity lines
x
The basic idea of the Godunov scheme for a general initial value problem is to approximate
the initial data by a piecewise constant function and then apply the self-similarity
structure to the series of Riemann problems.
There are two basic issues we have to consider immediately in this kind of ap-
proach. First, since the self-similarity for a piecewise constant solution holds locally
in time only, the structure of the Riemann problem can be applied for a small time
period. In other words the scheme is not free from the CFL condition and, hence, the
scheme can march just a little amount of time every time step and it costs computation
time. Furthermore, since rarefaction waves appear immediately, the solution
is not piecewise constant anymore. So a numerical scheme contains a process which
rearranges the rarefaction wave into a piecewise constant function every time step.
The numerical viscosity enters in this process and tracking down the behavior of the
scheme becomes extremely hard.
IMA, University of Minnesota, Minneapolis, MN 55455-0436 (yjkim@ima.umn.edu).
LeVeque [15] considers a large time step technique based on the Godunov method
for the genuinely nonlinear problem. In the scheme the CFL number may go beyond
1, and it is even possible to solve the propagation of a simple wave in a single step,
for the given final time T ? 0. However the scheme handles interactions
between waves incorrectly if the CFL number is so large.
One way to avoid the rearranging process is to consider a modified equation,
where h and u 0 approximate f and v 0 respectively. Dafermos [6] considers a polygonal
approximation h  f , i.e., h is a continuous, piecewise linear function. In the case
the exact solution of (1.4) is piecewise constant. So the method does not require
a rearranging process and, hence, it does not introduce numerical viscosity and the
error is controlled by taking the polygonal approximation h. In this approach the
exact behavior of the solution can be monitored more closely and we may get a more
detailed understanding of the numerical scheme based on this approach. This idea
has been developed in Holden and Holden [11], and it has been extended to multi-dimensional
problems in Holden and Risebro [13] and to systems of conservation laws
in Holden, Lie and Risebro [12]. In particular we refer Bressan [2], [3] for systems.
This front tracking method has been developed, especially by the Norwegian School,
as a computational tool.
Lucier [17] approximates the actual flux f by a piecewise parabolic function h and
achieves a second order scheme. In the case the initial data v 0 (x) is approximated by a
piecewise linear function u 0 and the solution remains piecewise linear. The difference
between the solutions of the original problem (1.1) and the modified problem (1.4) is
estimated by
Since the linear approximation is of second order, he achieves a second order scheme
for a fixed time t ? 0.
If we want to design a numerical scheme which represents the exact solution, we
have to find a way to choose the grid points correctly. If they are simply fixed, it
is clear that the scheme can not represent the exact solution and, hence, we need to
rearrange the solution to fit the solution to the fixed grid points. So it is natural
to consider moving mesh method, see Miller [18]. In Lucier [17] the moving mesh
method is used to find the exact solution of (1.4), where mesh points move along
characteristics. Another option is not to use any grid point. In numerical schemes
based on the front tracking method we mentioned earlier grid points are used just for
the initial discretization. The scheme we consider in this article does not use any grid
point, neither.
This article has two goals. The first one is to introduce the mathematical idea
which is behind the piecewise self-similar solutions. The second one is to demonstrate
how to implement the idea into a numerical scheme and show properties of the scheme.
From the study of the Burgers equation [14] it is easily observed that the primary
structure which dominates the evolution is a saw-tooth profile. In fact it is a series of
N-waves and eventually the solution evolves to a single N-wave, see Liu and Pierre [16].
The starting point of our scheme is to use this structure as the unit of the scheme.
This scheme has several unique properties that other schemes based on piecewise
constant functions do not have.
Piecewise Self-similar Solutions 3
Suppose that u(x; t) is a special solution of (1.1) which is a function of the self-similarity
variable x=t. Then the self-similarity profile (or the rarefaction wave),
easily derived from (1.1). It is natural to expect that characteristic
lines pass through the origin, i.e., they are compatible with self-similarity lines (1.3).
The piecewise self-similarity initial profile is considered in the sense that
Note that the time index t k can be a negative number. In this article we show that
the solution of (1.1) with piecewise self-similarity initial profile has such a structure
for all t ? 0, i.e.,
and give the explicit formula for this kind of solutions under several situations. First
we consider a convex flux with positive wave speed,
where f is locally Lipschitz continuous. The convexity of the flux f 00 (u)  0 is to get
the explicit formula 'g' of the self-similarity profile such that f 0 and the
self-similarity profile (1.7) can be written as
Note that the equality is included for the second derivative of the flux in (H) and,
hence, the monotonicity of f 0 is not strict and g is not exactly the inverse function of
f 0 and g(f 0 (u)) 6= u in the case. In this approach we may include a piecewise linear
flux of the front tracking method, see Remark 6.4.
In section 3 we consider a piecewise self-similar solution which can be written as
a self-similarity summation (or simply S-summation),
of finite number of base functions. We give definitions for the S-summation and base
functions in the section and show that u(x;
(x) is the solution
of (1.1) with initial data u
Theorem 3.6. We consider u
as an approximation of the solution v with the original initial data v 0 . Then the L 1
contraction theory of conservation laws implies
It is the estimate corresponding to the error estimate (1.5), which does not have the
time dependent term anymore. It is natural to expect that the error increases in time
if the flux is changed. In our approach we use the original flux and the error decreases
in time. The convergence of the scheme is now clear (see Theorem 3.6, Corollary
3.7). Note that the self-similarity summation (1.9) represents only special kind of
piecewise self-similar profiles (1.6), which have positive indexes t k ? 0 and are ordered
appropriately, i.e., c if an  b n  :::  a 2  b 2  :::  a 1  b 1 .
4 Y.-J. Kim
The self-similarity summation is coded for a numerical scheme successfully in Section
4. This scheme has several unique properties. First it does not require a time
marching procedure. So the complexity of the scheme is of order O(N ), not O(N 2 ).
Second it captures the shock place very well even if small number of base functions
(or mesh points) are used, Figure 4.3. In Figure 4.5 it is clearly observed that the solution
with finer mesh always passes through bigger artificial shocks and this property
provides a uniform a posteriori error estimate of the numerical approximation. Since
it does not introduce numerical viscosity at all, we may get very good resolution for
an inviscid problem. Our scheme also distinguishes physical shocks and artificial ones
clearly. Table 4.1 shows the time when the physical shock appears.
In Section 5 we generalize the method. For a general convex flux case,
the method is applied through the transformations (5.1) and (5.3). If the flux has
inflection points, then the scheme becomes considerably complicate and it is beyond
the purpose of this article. But, if the flux has only one inflection point, for example,
then we can easily apply the scheme through a similar transformation (5.4). Dafermos
considers a flux with a single inflection point through generalized characteristics.
The flux of the Buckley-Leverett equation satisfies this condition. The flux
which appears in thin film flows (see Bertozzi, Munch and Shearer [1]), also
belongs to this category. Figure 5.3 shows the strength of our scheme over the upwind
scheme in this case.
The scheme is not good enough for a short time behavior t !! 1 since the initial
error controlled efficiently. To resolve the situation we add an
extra structure to base functions in Section 6. Using these base functions we can
approximate the initial data with second order accuracy and still solve the exact
solution for the modified initial datum without the time discretization. Furthermore,
a general piecewise self-similarity profile (1.6) can be written in terms of self-similarity
summation of these modified base functions.
2. Self-similarity of conservation laws. Consider one dimensional scalar conservation
laws,
where the flux f is locally Lipschitz continuous. For a nonlinear flux f(u) the solution
may have a singularity and hence the solution is considered in the weak sense with
the entropy admissibility condition :
~
for any number ~ u lying between
a conservation law is from the fact that a rescaled function,
Piecewise Self-similar Solutions 5
is also the solution of (2.1) if and only if the initial profile u 0 (x) satisfies
It is clear that the Riemann initial condition,
ae
satisfies (2.4) and, hence, u(x; is a function of
the self-similarity variable,
The structure of a Riemann solution is given in Figure 2.1 together with characteristic
lines. Note that a self-similarity line R is not a characteristic line
and the solution is constant along it. This is a special property of Riemann problem
and it is not expected in a general situation. If the solution is constant along a line,
it is natural to assume that the line is a characteristic line and it is the starting point
of our scheme.
x
st
x
(a) Characteristic lines (b) Self-similarity lines
Fig. 2.1. Let f 0 (u+ lines are different from characteristic
lines. Even though the solution is constant along self-similarity lines.
If the total mass of the initial data u 0 (x) is finite,
Z
then the relation (2.4) cannot be satisfied since the transformation u
does not preserve the total mass. So the solution cannot be a function of self-similarity
variable x=t. In the following we consider techniques to achieve the Riemann
solution like self-similarity for general Cauchy problems.
In the Godunov method the space is discretized into small intervals and the initial
function u 0 is approximated by a step function which takes the cell average over those
intervals. Then the problem can be considered as a sequence of Riemann problems
and the structure of a Riemann solution holds locally in time and space. The scheme
finds the cell average of the solution after a small amount of time using the self-similarity
structure of Riemann solutions. It is fair to say that this method is more
6 Y.-J. Kim
focused on the structure of the Riemann initial data which makes the self-similarity
of the problem rather than the self-similarity itself. As a result the method takes
cell-averages every time step and loses the accessibility to the exact solution.
In the front tracking method the nonlinear flux f(u) is approximated by a continuous
function hm (u) which is linear between points f k
example, and
then the initial datum is approximated by a piecewise constant function by taking
these values, not cell averages. Then every discontinuity propagates as an admissible
shock of the modified problem,
in the sense of entropy condition (2.2) until it may possibly collide to other shocks.
We may say that the self-similarity of the original problem (2.1) has been modified
to get it fit to the piecewise constant functions. In this approach the exact solution
of the modified problem is accessible and, hence, the method can be employed as an
analytical tool as well as a computational one.
Now we suggest a new approach which keeps the self-similarity globally in time.
Suppose that characteristic lines of the solution u(x; t) pass through the origin. Then
we have
Since the right hand side diverges as t ! 0, we consider the initial datum as the
profile at a given time simplest case of L 1 initial datum of the kind is
Characteristic lines of this initial profile are given in Figure 2.2. Non-vertical characteristics
pass through the point (0; \Gammat 0 ) and there is a region in which characteristic
lines overlap with each other. The solution is given by finding the shock characteristic
correctly. In this case the shock characteristic is not a straight line
and the solution is not a function of x=(t+ t 0 ). Even though the solution is a function
of in the region
Since the shock speed s 0 (t) satisfies the Rankine-Hugoniot jump condition, the
shock place s(t) can be decided by its integral form. On the other hand, if the
convexity of the flux f is assumed,
we may consider the self-similarity profile g such that f 0 x. In the case
on (0; s(t)) and we can find the shock place s(t) easily from
the equal area rule,
(2.
Piecewise Self-similar Solutions 7
x
Fig. 2.2. Characteristic lines of a self-similarity solution are similar to self-similarity lines.
The main difference is that the shock characteristic is not a straight line anymore.
Since the conservation law (2.1) does not depend on the x variable explicitly, we
may translate the initial data (2.9) in the x direction. We can also consider initial
data which consist of finite number of these structures. A simple case is
where centers c k and shock places s k satisfy
The time index t k in (2.12) decides the slope of the initial profile and they can be
chosen differently. Condition (2.13) implies that all profiles in (2.12) are separated.
If not, the simple summation in (2.12) breaks down the self-similarity structure we
want to keep. In Section 3 we consider a self-similarity summation which preserves it.

Figure

2.3 shows characteristic lines for initial data (2.12) with 4. In this case
tracking down a shock is more complicate and (2.11) is not enough for the purpose
since waves interact to each other.
x
Fig. 2.3. Shock characteristics of (2.12) interact together and make a bigger shock.
3. Piecewise self-similar solutions. In this section we give the definition of
the self-similarity summation and show that the exact solution of (2.1) is given as
an S-summation. Notations in this section are directly converted into a numerical
scheme in Section 4. In this section we consider a flux under the hypothesis,
8 Y.-J. Kim
In this case the self-similarity profile 'g' is the profile which satisfies f 0 As
it is mentioned earlier g is not exactly the inverse function of f 0 since the monotonicity
of f 0 is not strict. We also assume f 0 this section for our convenience, and
it implies that the solution is actually assumed to be positive under (H). The result
of this section are generalized in Section 5.
3.1. Base functions. As it is mentioned earlier, the self-similarity profile
represents the asymptotic behavior of the conservation law (1.1). The function,
ae g
serves as a base function in this article. A base function has the self-similarity profile
over the interval between the center 'c' and the shock place `s'. The area (or the mass)
enclosed by the x-axis and the base function is given by
c
Z s\Gammacg(x=t)dx =: m(t; c; s):
It is convenient to consider the mass m as the fourth index of the base function, say
m;t;c;s (x), or any three of them as an index set. In any case we consider it under the
assumption that indexes m; t; c; s satisfy the relation (3.3). So if any three of them
are given, the fourth one is decided by the relation.
Consider a Cauchy problem,
It is clear from (2.10) that the solution u(\Delta; t) has the self-similarity profile with time
index t+t 0 between the original center c 0 and a new shock place s(t). Since the initial
total mass m 0 should be preserved, the solution of (3.4) is
where the shock place decided by the relation (3.3).
Remark 3.1. If we take a ffi-function as the initial datum, for example
the solution is given by u(x; t 0
(x). So the slope of the
base function represents the time of the evolution starting from the ffi-function like
initial data, and that is why we take index t for the base function.
Remark 3.2. For the Burgers case, the self-similarity profile is
given as the identity function, x. In the case (3.3) gives following relations,
Remark 3.3. The rescaling (2.3) does not preserve the total mass. So it can
not measure the invariance property for L 1 solutions of conservation laws. For the
Burgers equation we consider
Piecewise Self-similar Solutions 9
where the rescaling preserves the total mass. We can easily check that variables
are invariant under the rescaling after the translation t. These
variables are called self-similarity variables for L 1 Cauchy problems, and the Burgers
equation is transformed to
We can easily check that Bm0 ;t 0 =1;c0=0 (i) is an admissible steady state of the equation
and hence w(i; is the solution of (3.9). If we transform the
variables back to u; t; x, then we get u(x; This is another way
to show (3.5). In this example we can see that the approach with piecewise self-similar
solutions captures the self-similarity of the general Cauchy problems exactly. For a
detailed study for the transformed problem (3.9) we refer [14].
3.2. Self-similarity Summation. Since the solution of (3.4) is given by (3.5),
we can easily guess that
is the solution of the conservation law with initial data
if all the supports of the base functions in (3.10) are disjoint. But it is not usually
the case since the support of a base function expands in time. The self-similarity
summation(or simply S-summation),
is to handle the case that supports of base functions overlap with each other. The
definition is given inductively in the following.
(x). Suppose that
(x) is well
defined and supp(B
c
. Suppose there
exists a point  j ? c j such that
Under the assumption of (3.13), the left hand side of (3.14) is monotone in  j and,
hence, such a point is unique. If there is no such a point we say that the S-summation
(3.12) is not defined. If there exists such a point
ae g
Base functions are ordered by centers c k and then the S-summation is given from
the right hand side. It is because of the positiveness assumption for the wave speed,
(H). If the order of the summation is changed, the result is different. So
the S-summation is not associative.
Remark 3.4. If the time indexes are identical, t then we can
show the S-summation (3.12) is well defined. If then, since the self-similarity profile g
is an increasing function, we have g
has values of g
(3.15), the inequality (3.13) is
satisfied for all  j 2 R. Furthermore the left hand side of (3.14) has value
for and diverges to 1 as  j !1. So there exists a point  j satisfying (3.14)
and the S-summation is well defined.
Remark 3.5. We may consider  j as the j-th shock generated by the base
function Bm j
. Suppose that  j i.e., the j-th shock caught up the (j-1)-th
shock. The definition (3.15) implies that the self-similarity profile g
disappears. We can easily check that we will get the same S-summation (3.15) if we
remove the (j-1)-th base function and modify m j by adding . This property
represents the irreversibility of the conservation laws.
Theorem 3.6. Suppose that the flux f(u) satisfies Hypothesis (H). If the self-similarity
summation
(x) is well defined, then u(x;
+t;ck (x) is also well defined and it is the solution of (1.1) with initial
data u 0 . If v(x; t) is the entropy solution of (1.1) with initial data
Proof. The proof is completed through inductive arguments. Suppose that
+t;ck (x) is the solution with the initial condition
(x). It is assumed that u j (x;
(x) is well defined
and we let u j (x; t) be the solution of (1.1) with this initial data. Let
the shock characteristic given by the j-th base function, i.e.,  j for the  j in
(3.13,3.14). If x ?  j (t), then it is clear that u j (x; since characteristics
on the right hand side of do not interact with it because f 0 (u)  0.
since the vertical characteristics starting in the region x
do not touch shock characteristics moving to the right hand side. The characteristic
passing through a point (x; is a straight line connecting
and, hence, u
. Since the total mass is preserved, the shock
place should satisfy
+t;ck (x) from the definition of the S-summation and the first
part of the proof is completed. The second part (3.16) is simply the L 1 contraction
theory for conservation laws.
In the proof we employ the theory of characteristics (see [8], ch. 11). The error
estimate (3.16) implies that the initial error decreases in time, and the solution with
modified initial data is obtained exactly in a single step for any given time t ? 0. The
scheme has ideal properties for the study of asymptotic behavior.
Now we consider u
(x) as an approximation of L 1 initial
data v 0 . Let a partition be the set of centers. Its norm is defined
Piecewise Self-similar Solutions 11
by j. There can be many ways to discretize the initial data. To
guarantee the convergence of the scheme, we need the existence of ffi; L ? 0 such that
where a constant " ? 0 is given. An example of such a discretization is given in
Section 4.1. The convergence of the scheme satisfying (3.17) is clear from (3.16).
Corollary 3.7. (Convergence) The scheme of the self-similarity summation
(x) with initial discretization u
satisfying (3.17) converges to the entropy solution v(x; t) with initial data
as
Remark 3.8. Now we consider the S-summation between two base functions,

Figure

3.1. It gives a good example to figure out
the meaning of the S-summation. Furthermore, in the numerical computation, we
can possibly compare only two base functions each time and, hence, it is worth to
consider it in detail. If these two base functions are separated, s then the shock
place  of the definition (3.15) is simply
Z
implies that two base functions are merged, i.e.,K
Suppose that s 1 is far away and the shock place  is guaranteed to be between s 2
and s 1 . Then (3.18) can be written as
Z c1
Z
The solution  of (3.19) has a special meaning in the coding. We define  as an
operator between two base functions, Bm2
:= . Note that in the
definition of the operator we do not use the information m 1 at all. We just assume it
is big enough and, hence,  ! s 1 . This operator is used in Section 4 to check if two
adjacent base functions are merged or not.
For the Burgers case, implies that, in Figure 3.1, Trapezoid
has the same area as Triangle Ac 1 , and the relation (3.19) can be written as
an algebraic relation,
The operator ' ' between two base functions is now given explicitly,
c 2+s 2\Gamma2c 2 s2
which is the solution of the algebraic relation (3.20) with
We consider this operator in Section 4.1 again.
A
ff
Fig. 3.1. The equal area rule gives the shock place when two base functions interact together.
4. Coding Strategy. In this section we show how the self-similarity summation
can be implemented into a numerical scheme. To see what is really happening in each
step it is helpful to consider a specific example. For that purpose we consider the
Burgers equation,
The result of the scheme is compared with the Godunov scheme.
4.1. Implementation. Here we introduce a grid-less scheme based on the self-similarity
(a) The Equal Area Rule (b) Base Functions with overlaps
Fig. 4.1. Initial data are approximated by a piecewise self-similarity profile. It turns out to be
an S-summation of base functions.
Step 1. (Initial discretization) The first step is to design a method to approximate
the initial datum v 0 (x) by a self-similarity summation u 0 (x) which satisfies
(3.17). Consider n base functions B[k]; Each element B[k] consists
of two members B[k]:m; B[k]:c, which represent the area and the center of the base
function. We use identical time index t hence, we do not need a
member for the time index.
Piecewise Self-similar Solutions 13
0 be a cell average approximation of v 0 with steps of mesh size
profiles with time index t 0 ? 0 which pass through
the left end points of the constant parts of the step function v "

Figure

4.1 (a). Let
B[k]:c be the x-intercept of the k-th self-similarity profile from the right hand side
and B[k]:m be the area enclosed by x-axis,
, and the k-th and the
profiles. This discretization is well defined only if B[n]:c ! ::: ! B[1]:c. To
achieve it a small initial time index t 0 should be chosen depending on the initial data.
Since the initial self-similarity profile of the example (4.1) is a line with the slope
1=t 0 and the slope of the initial data is bounded by v x (x; we have to take t 0 ! 1.
In

Figure

4.1 the initial data in (4.1) has been discretized using 10 base functions,
base functions (b) have
some overlaps and the self-similarity summation (a) has a saw-tooth profile. The size
of the triangle like areas added and subtracted by the approximation is proportional
to " 2 and the total number them is proportional to 1=". So we have jjv
Theorem 3.6 says u(x;
is the solution with the
modified initial data u 0 . So the rest of the scheme is focused on how to display the
given solution. Even if it is possible to follow the inductive arguments of the definition,
we will get serious complexity in the coding if behind shocks capture the front ones,
i.e., . In the case the S-summation is not changed even if two base functions
are merged, Remark 3.5, and hence we do the merging process first. From now on
the corresponding time index is t for each k.
Step 2. (Merging) The operator ' ' between two base functions defined by (3.19)
for the general case or by (3.21) for the Burgers case plays the key role here. Suppose
that there is no contact between shocks for k
Then we can easily check that the k-th shock in (3.15) is given by
1. Suppose that
In the case  j 6= B[j]   B[j \Gamma 1] in general. Even though it implies  j ?  j \Gamma1 and the
self-similarity profile of the (j-1)-th base function B[j \Gamma 1] disappears, Remark 3.5. In
the case these two base functions B[j] and B[j \Gamma 1] can be combined, i.e., put
remove then rearrange the array B[\Delta] from
is the number of base functions left after the previous step. Since the combined base
function may take over another one again, we decrease the index j if j 6= 2. If (4.2)
does not hold, we increase index j. We continue this procedure from
Note that there is no base function B[0] and we use B[1]   B[0] := B[1]:s in (4.2) for
given by the relation (3.3).(Step 2 is complete.)
If there are no base functions merged together, there will be
of (4.2). If m base functions are merged, then base functions are left and the
maximum number of the comparison (4.2) is n In Figure 4.2 (b) base
functions at time are displayed after the merging process. There were 50 base
functions initially (a) and 38 of them are left after the merging step. It means that
small base functions has been merged together and made a big base function. The
big base function can be considered as an accumulation of small artificial shocks in
some sense and it represents the physical shock.
14 Y.-J. Kim0.050.150.25
(a) 50 Initial Base functions (b) 38 Base functions left at
Fig. 4.2. The initial base function with slope 1=t 0 has slope 1=(t0 + t) at time t without area
change. After the merging process, Step 2, some of the base functions are merged together and make
a big base function which represents a physical shock.
Now we are ready to display the solution. Suppose that base functions B[j];
are left after the merging step. Let  Then the right and
the left hand side limits are given by,
So to display the solution it is enough to plot the points
Between these point the solution has the self-similarity profile. So if
we connect these points with self-similarity profile with time index
B[j]:c, we get the solution. In Figure 4.3 solutions are displayed using different number
of base functions. We can clearly see that the solution converges as the number of
base functions is increased.0.050.150.250.350.45
(a) Initial Discretization (b) Solutions at
Fig. 4.3. Three S-summations using 10,40 and 160 base functions. The solution finds the
shock correctly even if a very rough initial discretization is given. A solution with finer mesh passes
through artificial shocks.
4.2. Comparison with Godunov. A typical way to discretize the initial data is
taking the cell average, Figure 4.4 (a). The Godunov scheme solves Riemann problems
between each cells for a short amount of time \Deltat and then repeat the process until it
Piecewise Self-similar Solutions 15
reaches a given time t ? 0. In Figure 4.4 (b) we can see that the numerical solution
converges to the same limit as the S-summation, Figure 4.3 (b), as \Deltax ! 0.0.050.150.250.350
(a) Data Discretization (b) Solutions at
Fig. 4.4. Three approximations by Godunov using 1=160. The scheme is
convergent to the same limit of the S-summation. We can observe that numerical solutions are
separated near the shock and it is hard to guess where the limit is from a single computation.
Remark 4.1. (Computation time) Let N be the number of mesh points. Then
the number of operations for the S-summation is of order N since the time marching
process is not required, Theorem 3.6. The number of operations is almost independent
from the final time t ? 0. On the other hand the Godunov scheme has operations of
and the situation becomes worse if the final time t is increased.0.220.260.30.34
Fig. 4.5. A magnification of Figure 4.3 (b) near the physical shock shows that self-similarity
solutions with finer mesh passes through the middle of artificial shocks.
Remark 4.2. (Error estimate) We can clearly see that the exact solution v of
limit of the S-summation) always passes though the artificial shocks
of self-similarity solutions, Figure 4.3 (b). This property makes it possible to get a
uniform a posteriori error estimate. Figure 4.5 is a magnification of Figure 4.3 (b).
There are couple of other things we can observe here. First the sizes of artificial shocks
decrease in time with order of O(1=(t We can also observe that, even if we use
small number of base functions, we can get the physical shock very closely.
Remark 4.3. (Shock Appearance Time) In a numerical scheme the solution
is approximated by piecewise continuous functions and hence it is hard to see if a

Table
Shock Appearance time. The exact solution with initial data (4.1) blows up at 1. The time
of shock appearance can be measured by counting the base functions after the merging step.
initial number of base functions the time when the number is decreased by 2
200 T=1.0015
discontinuity represents the physical shock or not. In our scheme, as we can see from

Figure

4.2, the accumulation of base functions represents the physical shock. So if a
base function is merged by its behind one in the sense (4.2), we may conclude that a
physical shock has appeared. The physical shock appears at time in the example
since min(@ x v 0 \Gamma1. We can easily check that if (4.2) happens around the
time.

Table

4.1 shows the time when the number of initial base functions decreases.
5. General cases. The self-similarity summation has been considered under
Hypothesis (H). In this section we generalize it under Hypothesis (H1) and (H2).
5.1. General convex flux. We consider L 1 initial function u 0 which is uniformly
bounded, say \GammaA  u 0 (x)  B. Then the solution of (2.1) is always bounded,
Consider a convex flux,
If the flux satisfies f 00 (u)  0, we may change the variable get an equation
f satisfies (H1). Note that we include the
equality in (H1) and a piecewise linear flux can be considered.
We can easily check that a new flux,
satisfies the hypothesis (H) and h 0 be the solution of
We can easily check that
is the solution with the original flux f and initial data u 0 . Since u  \GammaA, the solution
v(x; t) is positive. Now we are in the exactly same situation as in the previous sections
except the structure of the initial data. The initial data v(\Delta; 0) is not L 1 anymore. To
handle the situation we consider two special base functions with infinite mass,
Piecewise Self-similar Solutions 17
These base functions handles the transformation u A. Note that the
speed of the shock connecting the state our case.
The Self-similarity summation including these two base functions can be defined in
a similar way. We omit the detail. Figure 5.1 shows how the self-similar solution
evolves for the Burgers case. In the figure even the solution with very rough initial
discretization with only 16 base functions represents the asymptotic behavior very
correctly.
-0.4
(a) Data Discretization (b) Solutions at
Fig. 5.1. Three S-summations are displayed using 16,64 and 256 base functions. It handles
sign changing solutions correctly. This figure shows the time convergence to an inviscid N-wave.
5.2. Flux without convexity. Consider a flux with a single inflection point,
Then, under the change of variables,
the problem (2.1) is transformed to
Then the new flux h satisfies
and A is not the lower bound of the solution u(x; t) in
general, we can not expect v  0. So in this case we have to consider positive part and
negative part together. It is possible since h 0 (u) is monotone on (\Gamma1; 0) and (0; 1)
respectively. All we have to do is to consider negative base functions together with the
positive ones. Since the wave speed h 0 (u) is positive, the self-similarity summation is
defined from the right hand side as in the previous cases.
Example 5.1. Consider an inviscid thin film flow in [1],
where the initial datum is compactly supported supp(u 0
has a single inflection point under the transformation (5.4), we
get the flux It satisfies
which is not exactly same as (5.5) but has the opposite direction in the inequalities.
We do the self-similarity summation from the left hand side instead of changing the
space variable using \Gammax. Now the original problem (5.6) is transformed into
In this case the self-similarity profile (2.8) is given by,
and the corresponding base functions are,
ae
The initial data v 0 (x) converges to \GammaA as x ! \Sigma1 and we need to consider two base
functions with infinite mass,
Note again that, in our example (5.6), the infinite state is \Gamma1=3 and the shock
speed is
Numerical solutions of (5.6) with initial data,
are in

Figure

5.2. The first picture shows the initial data and the self-similarity
summation using 200 base functions. A part of it has been magnified with numerical
approximations of upwind scheme in the second picture. We can clearly see that the
solution of upwind scheme converges to the self-similarity summation. This example
shows that the self-similarity summation gives a very accurate resolution using small
number of mesh points. Furthermore, since it gives the solution without time marching
procedure, computational time is a lot smaller.
5.3. Flux with the space dependence. Since the self-similarity of the problem
(2.1) depends on the fact that the flux depends on the solution only
we have no clue how to generalize our scheme to a problem with a general space
dependent though, if the space dependence is given by
the equation is transformed to
under the change of variable
R x1=a(s)ds and our scheme can be applied.
Piecewise Self-similar Solutions 190.10.30.50.7
(a) initial data and S-summation at t=6 (b) comparison with upwind
Fig. 5.2. Flux is Picture (a) shows the initial data and the self-similarity
summation at shows that upwind converges to the self-similarity summation. 200
base functions are used in the S-summation and 800 and 4,000 meshes are used in upwind scheme.
Since the self-similarity of hyperbolic conservation laws is the one-dimensional
property, it should be possible to expand the scheme to multi-dimension problems.
Consider a 2-dimensional problem,
with a velocity vector field satisfying
Cvetkovic and Dagans [5] suggest space variables y 1
dy 1
dj
which transform (5.14) to
Problem (5.16) can be considered as a set of one-dimensional problems and, hence,
the complexity of the scheme for it is of order O(N 2 ). Since the transformation (5.15)
also has the complexity of O(N 2 ), we eventually get a scheme of O(N 2 ) for a two-dimensional
problem. In this approach each channel of the velocity vector field is
considered separately and, hence, it seems useful to channel problems.
6. Second order approximation. The scheme introduced in the previous sections
solves the problem exactly with modified initial data, and the size of the initial
error decreases in time. Even though the scheme is not good enough for the short
time behavior since the error generated by the initial discretization can be huge. Here
we add an extra structure to base functions and make the initial data discretization
to be second order. In this way we can handle general self-similarity solutions (1.7).
6.1. Modified base functions. The base function considered in the previous
sections has three indexes, say m; t; c. In this section we introduce two more indexes,
h and  t. Note that there are two time indexes t and
which play different roles. We
assume 1. For the simplicity we consider under the
hypothesis (H). It can be easily generalized as we did in Section 5.
To figure out the structure of the new base function B h;  t
m;t;c (x), we consider
and
Let g be the self-similarity profile, f 0 As an intermediate step we define
t;c (x) first. For
defined by
and, for it is defined by
The constant  c is the center of the top self-similarity profile with time index  t
and the constant x   is the x-coordinate of the intersection point between two self-similarity
profiles with index t an  t. We can easily see from (6.2) that  c ! x   for  t ? 0
and
t;c (x) is well defined for since the
corresponding domain is empty. For
Now we introduce the index m ? 0 which decides the support of the base function.
c be the solution of,
Z
c
For it always has a solution. For  t  0 it has a solution only
R  c
t;c (x)dx.
The base function is now defined by
t;c
The self-similarity summation among these base functions can be similarly defined
using the profile g
\Delta in the domain c ! x ! x   and the profile g
for x   ! x. We omit the detail. We may consider the base function (3.2) as a special
case of (6.7) with
6.2. Initial discretization and the exact solution. Suppose the initial function
be a partition of the interval [A; B]. We can approximate v 0 with
self-similarity profiles over interval
which is second
order. For the Burgers case it is simply a piecewise linear approximation. The
approximation u 0 can be written as
Piecewise Self-similar Solutions 21
R
Initially the supports of base functions
are disjoint and, hence, the self-similarity summation is the usual summation. The
exact solution of the conservation law with initial data (6.8) is
We still consider the exact solution and the contraction theory implies
Remark 6.1. The initial discretization (6.8) is trivial in comparison with Step
1 in Section 4.1. It is an additional advantage when the modified base function is
used in a numerical scheme. Even though this additional structure may cause extra
complexity when it is used as an analytical tool.
Remark 6.2. (Piecewise Constant Data) In many cases initial data are given
as piecewise constant functions from the beginning. In the case an initial datum can
be considered as a summation of base functions with In Figure 6.1 we
consider the Burgers case (4.1) using base functions B h;1
m;t;c (x). We can clearly see
that these approximations represent the shock place very well. Unlike the previous
case, the solution with finer mesh always passes though the constant parts of coarse
(a) Data Discretization (b) Solutions at
Fig. 6.1. The S-summation for the modified base functions (6.7) with
constant, piecewise self-similar solution. In the figure 3 summations are displayed together using
base functions. We can observe that finer one always passes the constant parts.
Remark 6.3. (Singular Initial Data) If singular initial data are given, then extra
mesh points are usually introduced to capture the effect of the singularity of the data.
But, since our method handles initial data individually, extra mesh points are not
needed. In Figure 6.2 the Burgers equation is solved with singular initial data (a).
We use 6 modified base functions with
Remark 6.4. (Front Tracking) It is possible to consider the front tracking
method in terms of the self-similarity summation. Consider an L 1 solution of the
Burgers equation bounded by 0  u(x; t)  1. Let h(u) be the polygonal approximation
of the flux with the partition f0; 1=n; :::; 1g. So h 0 (u) is a
step function,
22 Y.-J. Kim0.050.150.250.350
(a) Singular initial data (b) Solutions at
Fig. 6.2. The scheme does not require extra meshes to handle singular initial data (a). In the
S-summation every datum is handled exactly by a base function. Only 6 base functions solves this
example.
and the self-similarity profile g(x) is also a step function,
So the values of g(x) are the breaking points of the flux h(u). We can approximate
the given initial data v 0 by taking a cell average, not just breaking points. Then the
initial discretization u 0 can be written in a from of (6.8) with
1. This is a
simplified version of the front tracking method under Hypothesis (H).
7. Conclusion. The basic idea of the method introduced in this article is to
approximate the solution of a conservation law by a self-similarity summation of
base functions. In that approach we get the exact solution in the class of functions.
This method can be easily converted into a numerical scheme and the complexity
of the scheme is of order N , not N 2 since no time marching procedure is needed.
Convergence of the scheme is now a trivial matter, Theorem 3.6 and Corollary 3.7.
The method can be used as an analytical tool. In fact the author is preparing
an article on asymptotic behavior of scalar conservation laws through this method.
Various issues appear when we apply this idea to other cases, systems or convection-diffusion
equations. The author does not have a good understanding for these cases
yet.

Acknowledgement

The author would like to thank Professor A. E. Tzavaras.
He gave the author the motivation and valuable remarks for this work. The author
also would like to thank people in IMA for all the discussions and supports.



--R




On the partial difference equations of mathematical physics

Polygonal approximations of solutions of the initial value problem for a conservation law
Regularity and large time behaviour of solutions of a conservation law without convexity
Hyperbolic conservation laws in continuum physics
Solutions in the large for nonlinear hyperbolic systems of equations
A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics
On scalar conservation laws in one dimension.
An unconditionally stabl method for the Euler equations
A method of fractional steps for scalar conservation laws without the CFL condition
Diffusive N-waves and Metastability in Burgers equation
Large time step shock-capturing techniques for scalar conservation laws

A moving mesh numerical method for hyperbolic conservation laws

--TR
characteristics;gridless scheme;front tracking;self-similarity