The Riemann problem for a generalized Burgers equation with spatially decaying sound speed. I Large‐time asymptotics

In this paper, we consider the classical Riemann problem for a generalized Burgers equation, ut+hα(x)uux=uxx,$$\begin{equation*} u_t + h_{\alpha }(x) u u_x = u_{xx}, \end{equation*}$$with a spatially dependent, nonlinear sound speed, hα(x)≡(1+x2)−α$h_{\alpha }(x) \equiv (1+x^2)^{-\alpha }$ with α>0$\alpha >0$ , which decays algebraically with increasing distance from a fixed spatial origin. When α=0$\alpha =0$ , this reduces to the classical Burgers equation. In this first part of a pair of papers, we focus attention on the large‐time structure of the associated Riemann problem, and obtain its detailed structure, as t→∞$t\rightarrow \infty$ , via the method of matched asymptotic coordinate expansions (this uses the classical method of matched asymptotic expansions, with the asymptotic parameters being the independent coordinates in the evolution problem; this approach is developed in detail in the monograph of Leach and Needham, as referenced in the text), over all parameter ranges. We identify a significant bifurcation in structure at α=12$\alpha =\frac{1}{2}$ .


INTRODUCTION
In this paper and its companion, we investigate the effects of a spatially varying nonlinear sound speed in a generalized form of the classical Burgers equation. Specifically, we consider the situation when the nonlinear sound speed decays algebraically with increasing spatial distance. A discussion of generalized Burgers equations, and their significance in nonlinear acoustics, in the present context, can be found in the papers by, among many others 2-7 and, 8 and very recently, 9 together with the classical texts of Whitham, 10 and Rudenko and Soluyan. 11 A number of recent papers have addressed (using, among other approaches, Cole-Hopf type transformations, similarity and separable structures, and ad hoc approximations) the generalized Burgers equation, with time-dependent viscosity, 12 and linear damping, 13 while the inclusion of a forcing term is considered in Refs. [14] and [15]. Higher-dimensional considerations are made in Ref. [16], while the Burgers/KdV fusion is addressed in Ref. [17]. In relation to nonlinear acoustics, the spatially decaying sound speed considered here represents finite amplitude sound propagation in a spatially varying, inhomogeneous medium, in which, relative to the constant background sound speed, the propagation of finite amplitude acoustic disturbances has sound speed ( , ) = ( ) + ( 2 ) as → 0, with being the dimensionless acoustic disturbance and measuring dimensionless distance from a fixed origin. The function ( ) measures the inhomogeneous local variation of the finite amplitude sound speed in the material background. Qualitatively, this type of spatially decaying inhomogeneity in the local sound speed could be effected in, for example, a long tube of gaseous material which has a prescribed background local density which is bounded and decreases with distance, from a local maximum at the spatial origin. To model this situation, we make a specific choice of spatial variation in the local sound speed. This choice is motivated by a number of observations. First, it is a generic representative for any smoothly decreasing variation which decays algebraically with distance at long range. Second, for any more general decreasing variations which decay more rapidly with distance than algebraically, the results here obtained will provide upper bounds, via the parabolic comparison theorem. Our aim is, therefore, to develop a detailed analysis of the Riemann problem associated with this class of generalized, spatially inhomogeneous, Burgers equation.
We denote the Cauchy problem given by (2)-(5) as (BIVP). It should be noted that, via (4), at each ∈ ℝ at which 0 is continuous, then ( , ) → 0 ( ) as → 0 + . Moreover, when 0 is continuous for ∈ [ 1 , 2 ] ( 1 > 2 ), then ( , ) → 0 ( ) as → 0 + , uniformly for ∈ [ 1 , 2 ]. In relation to applications, we may consider the generalized, spatially inhomogeneous Burgers equation in (3) as representing a rational model for the unidirectional propagation of a transitional, weakly nonlinear acoustic wave through a viscous medium with a spatially varying local sound speed, which is weakly decaying with distance, from a localized maximum. In this paper, we focus attention on the Riemann problem for the Burgers equation (3), and thus we consider (BIVP) when the initial data 0 ∶ ℝ → ℝ has the specific form: with + , − ∈ ℝ and + ≠ − . On using (6), condition (4) can be written as: uniformly for ∈ ℝ. The aim of the present paper is to examine in detail the large-structure of the Riemann problem given by (BIVP) with (7). Generic qualitative results relating to the more general form of (BIVP) will be presented in the second of this pair of papers. At this stage, it is instructive to simply quote the results from this second paper, as they apply to (BIVP) with the specific initial data (7), and as they apply to the objectives of the current paper. Specifically, we can summarize these results in the following: There exists a global solution ∶ ∞ → ℝ to (BIVP). This solution is unique. In addition, for each > 0, the solution depends continuously on the parameters + and − , in . Also, for each > 0, ( , ) → ± as → ±∞, uniformly for ∈ (0, ], while (⋅, ) > 0(< 0) when − < + ( − > + ) for all ∈ (0, ∞).
A proof of this theorem follows directly from Proposition 2.8 and Proposition 2.9 in the second of this pair of papers, 18 henceforth referred to as (MNII). To reduce the region of the ( + , − ) parameter space that we need to investigate, we also have: Proof. This follows on observing the invariance of the PDE (3) under the transformation ( , ) ↦ (− , − ), together with the uniqueness property established in Theorem 1 for (BIVP). ■ As a consequence of this result, we need only to investigate (BIVP) for ( + , − ) ∈ , with In addition, we note that when ( + , − ) = ( + , − + ), it is an immediate consequence of (8) that the corresponding solution to (BIVP) is an odd function of the spatial variable . As stated earlier, the specific objective of the present paper is the determination of the detailed large-asymptotic structure of the solution to the problem (BIVP), with particular attention being paid to the key attracting features in the solution as → ∞, and how these key features depend upon the parameters > 0 and ( + , − ) ∈ . To achieve this, we adopt the method of matched asymptotic coordinate expansions, as developed by Leach and Needham, 1 which, specifically, has been successful in similar applications to Riemann-type problems for classes of temporally inhomogeneous, generalized Burgers equations (see, e.g., Refs. [7] and [ 8 ]). Briefly, this approach uses the classical method of matched asymptotic expansions, but with the independent coordinates as the asymptotic parameters. This involves first constructing the asymptotic solution as → 0 + , uniformly for ∈ ℝ; this is then used, via matching, to construct the asymptotic solution as → ±∞ with = (1) + , respectively; this then provides, again via matching, the far field information required to enable the key asymptotic structure as → ∞, uniformly for ∈ ℝ, to be fully determined. We find that there is a significant bifurcation in the large-structure for (BIVP) as passes through the critical value = = can fall into either the diffusion front case or the expansion/shock wave case depending on the magnitudes of the parameters − and + .
The paper is constructed as follows. In Section 2, we begin the analysis by developing the smallasymptotic structure to (BIVP). Section 3 then uses the details obtained in Section 2 to develop the aymptotic structure to (BIVP) as → ±∞. Finally, in Sections 4-6, the detailed asymptotic structure to (BIVP) can be developed as → ∞, guided by the details from Section 3. In each section, illustrative numerical solutions are provided, for comparison with the asymptotic results. A discussion is presented in Section 7.

ASYMPTOTIC SOLUTION AS → +
In this section, we construct the asymptotic solution to (BIVP) as → 0 + , uniformly for ∈ ℝ. When is small, the form of the initial data, and a balancing of terms in Equation (3) so as to retain the diffusion term to smooth out the initial discontinuity in the initial data, leads to a structure which has three asymptotic regions (a central smoothing region, together with left and right far field connecting regions), as follows: Region I: = ( 1 2 ) and = (1) as → 0 + In this region, we introduce the scaled coordinate = ∕ 1 2 = (1) as → 0 + , and expand in the form: with = (1). After writing (3) in terms of the scaled variables, substituting from (10) and applying condition (7), we readily obtain, after an integration, We note that the parameter does not appear at leading order in this region. We now move on to complete the asymptotic structure with the two additional outer regions when (±) ≥ (1), respectively.

ASYMPTOTIC SOLUTION AS → ±∞
This section constructs the asymptotic solution to (BIVP) as → ±∞ with = (1) + . The form of the expansions here are obtained directly from the corresponding expansions in regions ± , respectively, when considered with | | large and = (1) + . Without presenting details, which are standard, we obtain, after matching with regions ± when ≪ 1, respectively, the following expansions in regions ± .

ASYMPTOTIC SOLUTION AS → ∞ WITH ∈ ( , )
This section concentrates on constructing the detailed large-structure of the solution to (BIVP) when ∈ (0, 1 2 ). We begin in > 0. An examination of expansion (20) with (18) shows that when ≫ 1 this expansion remains uniform for ≫ , but fails when = ( ) + . Thus, to continue the asymptotic structure, we introduce region + , in which = ( ) + as → ∞. It should be noted here that the form of the expansions in each spatial region which is introduced when is large is determined by examining the form of the expansion in the neighboring region, when written in the current regions scaled coordinates.
An examination of the expansion in region − when, = ( ) − , determines that we must write, after which we obtain, following matching with region − as → −∞, An examination of expansions (30) and (33) demonstrates that they become nonuniform when = ( − 2 1+2 ) ± (which corresponds to = ( 1 1+2 ) ± ), respectively. Thus, to continue the asymptotic structure, we must introduce two further regions. We begin in region + .
We first introduce the scaled coordinate = − 1 1+2 = (1) + as → ∞ in this region. An examination of the expansion in region + then requires us to write, after which we expand in the form, with = (1). On substituting from (35) and (36) into Equation (3) and expanding, we obtain at leading order the equation, and while matching to region + requires, The problem for Ψ + 0 ( ) will be completed, interactively, at a later stage in this section. However, with this in mind, we make some fundamental observations regarding Equation (37) with (38) and (39). The following result is readily established.

( + )
Each of the following statements hold for Equation (37) and For 0 = ((1 + 2 ) + ) 1 (1+2 ) , the conditions (40) continue to hold for > 0 , with now, while condition (42) continues to hold. Also, when + = 0 and 0 = 0, we have, exactly, there exists a unique solution on [ 0 , ∞) which again satisfies conditions (40) and (42), with, for 0 ≠ 0, while, for 0 = 0, The above statements in ( + ) are established by identifying the local existence and uniqueness of a solution for each initial value problem, which follows standard ODE theory, then establishing monotonicity of solutions to the initial value problem, and then using this to establish a priori bounds for the initial value problem, which allows us to infer global existence and uniqueness for the initial value problem. The far field asymptotic forms are derived directly from the ODE, making use of monotonicity. At this stage, we must leave region + , and return for its completion, interactively, at a later stage. For the present, we move on to introduce region − .
In this region, = (1) − as → ∞, and it follows from the form of the expansion in region − that we must write, and expand as, with = (1). After substituting into Equation (3), when written in terms of and , and expanding, we obtain at leading order, and the requirement, while matching to region − demands that, Again, the leading order problem must be completed, interactively, at a later stage. For the present, we may readily establish the following useful properties of Equation (49) with (50) and (51).
Each of the following statements holds for Equation (49) with and 1+2 , the conditions (53) continue to hold for < 0 , with now, while condition (54) continues to hold. Again, when − = 0 and 0 = 0, we have, exactly, there exists a unique solution on (−∞, 0 ] which again satisfies conditions (52) and (54), with, for 0 ≠ 0, while, for 0 = 0, The above results in ( − ) are obtained in much the same way as described for the earlier corresponding initial value problems in ( + ). The next, and key, region determines the selection of 0 for regions ± above, and acts as a "wavefront" region, effecting the transition from the − regions to the + regions. We first deal with the case when + < − . In this case, as a consequence of the expansions in regions ± , the transition is expected to take place locally, and we must now introduce region as follows.
Region ( + < − ): = 0 + (1) and = (1) as → ∞ It is in this region that 0 will be determined. We first write: with = (1) as → ∞ and > 0 to be determined. On expanding, with = (1), and substituting into (3), we find that a nontrivial leading order balance (recalling 0 < < 1 2 ) requires that, We note that, in terms of the original spatial coordinate , we now have, (62) After substitution from (60) and (62) into (3), we obtain, at leading order, the following balances: where, in obtaining the factorization on the right-hand side above, we must choose, In this exceptional case, 0 = 0 in (59), and the leading order balance changes, now requiring a rescaling, with the choice which dictates, in fact, that = (1) as → ∞ in this region, in this case, with the leading order balance now giving In both of the above two cases, matching with regions ± , respectively, leads to the boundary conditions, noting that replaces as the independent variable with − = − + (and + < 0) in the above conditions, when dealing with case ( ).
We begin with case ( ). It is readily established (recalling that here + < − ) that the boundary value problem (63) with conditions (67) and (68) has a unique (up to translations in ) solution. In fact, this can be readily obtained, by quadrature, as, This completes the large-structure in this case, with the determination of 0 in (64) allowing for the interactive completion of the structure in regions ± . The principal feature is a thickening and decelerating smooth shock wave, in region SW. In terms of the original spatial variable , the location of the shock wave is, with thickness, as → ∞. It is instructive here to use a simple central finite difference method, which leads to secondorder accuracy in space, to generate numerical solutions of (BIVP) for comparison with this asymptotic structure as → ∞. We truncate the domain of solution to − ≤ ≤ and 0 ≤ ≤ , with > 0 sufficiently large that the choice of has no effect on the solution when = within a chosen accuracy. We discretize the spatial domain at the equally spaced points ≡ − + ( − 1)Δ , with Δ = 2 ∕( − 1) and = 1, 2, … , . At these points, (3) becomes where ( ) ≡ ( , ). We assume that 1 ( ) = 1 (0) = − and ( ) = (0) = + . We solve this set of evolution equations using the explicit midpoint method to ensure second-order accuracy in space. We choose a timestep smaller than 1 2 Δ 2 , which ensures that the method is numerically stable. In all the calculations presented in this paper, we used Δ = 0.05 and a timestep equal to Throughout the paper, we run all numerical solutions up to times when the large-structure has evidently emerged. In Figure 1, we graph the numerical solution to (BIVP) on the ( , ) plane, when = 0.25 with ( + , − ) = (0, 1), at times = 100 − 2900 in increments of 100. We can readily discern the emerging, thickening, and decelerating smooth shock wave, which, at the final time is seen to be approaching the "tanh" profile predicted above.
In the first window of Figure 2 , as calculated from the numerical solution, and compare this with the asymptotic location given in (70), noting that in this case, 0 ≈ 0.82548. Although the numerically calculated value appears not to be in good agreement with the asymptotic behavior given by the broken line, this is deceptive because the rate of convergence is slow, of ( In the second window of Figure 2, we graph ( 1 2 ( ), ) 1 3 as calculated from the numerical solution, which is a measure of the scaled front thickness. Similar remarks regarding the slow approach to the asymptotic value (shown with a broken line) apply to this quantity, as well as quantities that we plot in later graphs. The value of the limit predicted by the asymptotic theory is We now move on to consider case ( ). On recasting (66) as an equivalent integral equation, it is straightforward to establish that the boundary value problem (66) with (67) and (68) has a unique solution, and this has properties, In addition, with ∞ > 0 being a globally determined constant, dependent on + (< 0) and . This completes the large-structure in case ( ), with the determination of 0 = 0 completing interactively the structure in regions ± . Again, the principal structure is a smooth shock wave, in region SW. However, in this case, the shock wave is stationary, at location = ( ) = 0, and has constant  Figure 3 at times = 100 − 2400, at intervals of 100. Convergence to a stationary, antisymmetric monotone smooth shock wave, with thickness of (1), is clearly evident. Each case when + < − is now complete, and we move on to consider the complimentary cases, when + > − . In replacement of region SW, we have the following:

ASYMPTOTIC SOLUTION AS → ∞ WITH ∈ ( , ∞)
In this section, we construct the detailed large-structure of the solution to (BIVP) when ∈ ( 1 2 , ∞). We again begin in > 0. An examination of expansion (19) establishes that, when ≫ 1, this expansion remains uniform for ≫ but fails when = ( ) + . Thus, to continue the asymptotic structure, we again introduce region + , which has exactly the same details as that in Section 4, except that here, expansion (30) is now replaced by, with = (1) + . Here, + ( ) is an indeterminate function in the large-structure, being globally related, but matching to region + requires, Similarly, we next introduce region − , when = ( ) − as → ∞, which again has exactly the same details as in Section 4, but now with expansion (33) replaced by, with, again, − ( ) indeterminate, while matching to region − requiring, An examination of expansions (87) and (89) reveals that these expansions become nonuniform, in this case, when = ( − 1 2 ) ± (which corresponds to = ( 1 2 ) ± ), respectively. Therefore, to complete the asymptotic structure, we must next introduce two further asymptotic regions. We begin in region + .
with matching to region + requiring, The solution to (92) with (93) is readily obtained as, with the free constant + to be determined. We note for future reference that, while matching to region + requires, We now move on to the corresponding region in < 0.
Here, we expand as, after matching with region − . Again, − is a free constant to be determined. We note that, while matching to region − requires, This completes region − . We note that, via Theorem 1, we can deduce that both of the constants + and − must be strictly positive. We now require a localized region to connect region − to region + across = 0. A balancing of terms in Equation (3), and noting that this localized transition region, centered at = 0, should retain the full inhomogeneous factor, together with examination of the expansions in regions ± , dictates that = (1) in this transition region, and while we must have = + ( − 1 2 ) as → ∞. The condition (101) provides the first relation between the constants ± as, Thus, we introduce: Region Tr : = (1) and = + ( In this region, we must expand in the form, with = (1). On substituting into Equation (3) and expanding, we obtain at leading order, Matching to regions ± then requires, A first integral of Equation (104) gives with being a free constant. An application of the conditions (105) and (106) then requires that , + , and − must satisfy the two equations, Solution of these two equations, using (102), then gives, and We now observe that and both ± are determined by the three simultaneous, nonlinear, algebraic equations given in (101), (102), and (110), in terms of the given parameters ± and . To examine these equations, it is convenient to introduce the new constants, We then have immediately, from (102), that after which Equations (101) and (110) become, After substituting from (115) for into (114), we obtain a scalar equation for Δ, namely, It is straightforward to establish that Equation (117) has a unique solution, say Δ = Δ * , which depends on each of ± and . In particular, In terms of Δ * , we now have, via (109), (111), and (112), while is given by (115) with Δ = Δ * . An integration of (107) finally gives, with ( + , − , ) and ( + , − , ) as given in (121) and (115), respectively, with the value of̄(0) requiring the determination of, and matching to, terms of ( − 1 2 ) in regions ± , which we do not We observe that in this special case, we can also infer from Corollary 1 that ( , ) is an odd function of ∈ ℝ, and thus we may conclude that̄(0) = 0 in (123). The asymptotic structure of the solution to ( ) as → ∞ is now complete for all cases with ( + , − ) ∈  when > 1 2 . The principal feature is the diffusion scale front regions ± , when = ( 1 2 ), connected by the localized steady region, when = (1). To conclude this case, we compare the large-structure developed above with illustrative numerical solutions of (BIVP). For the numerical solution of (BIVP), we set = 1. First, we consider ( + , − ) = (1, 0), and the numerical solution is illustrated on the ( , ) plane in Figure 6, where regions ± , which have thickness ( 1 2 ), and , which has (1) thickness, are readily seen emerging for t large. As a check on numerical accuracy, we have from the asymptotic solution that (0, ) = + ( − 1 2 ), with ≈ 0.288, while the numerical solution gives ≈ 0.2852. In Figure 7, we graph (0, ) against as calculated from the numerical solution, which shows good agreement with the ( − 1 2 ) decay to as → ∞, in accord with the asymptotic theory. We next consider the case ( + , − ) = (0, 1), with the numerical solution illustrated in Figure 8. In this case, the region Tr is more pronounced, with more asymmetry between regions ± , and from the asymptotic solution, we again have (0, ) = + ( − 1 2 ), with now ≈ 0.952, while the numerical solution gives ≈ 0.9402. In Figure 9, we graph (0, ) against as obtained F I G U R E 7 (0, ) against as calculated from the numerical solution to (BIVP) with = 1 and ( + , − ) = (1, 0). The dashed line is at the limiting value given by the asymptotic theory from the numerical solution to (BIVP), which illustrates the ( − 1 2 ) approach to the limit as predicted by the asymptotic theory.
Finally, in Figures 10 and 11, we illustrate the numerical solution for the cases ( + , − ) = (1, −1) and ( + , − ) = (−1, 1), respectively. In each case, the antisymmetry of the profiles about theaxis is evident. It should also be noted that the two cases are close, but not exactly, antisymmetric images of each other. This is in line with the asymptotic theory.
We now move on to consider the final, marginal, case when = 1 2 . In particular, it is of principal interest to determine whether this marginal case follows the shock/expansion wave structures seen when 0 < < 1 2 or the diffusion dominated structures seen when > . We again begin in > 0, with an examination of expansion (21) determining that, when ≫ 1, this expansion remains uniform for ≫ , but fails when = ( ) + . Therefore, to continue the asymptotic structure, we again introduce region + , which has the same details as in Section 4, except that now expansion (30) is here replaced by, Here, + ( ) is an indeterminate function in the large-structure, but matching with region + requires, Similarly, we next introduce region − , when = ( ) − as → ∞, which has exactly the same details as in Section 4, but now expansion (33) is replaced by, with, again, − ( ) indeterminate, whist matching to region − requires, Expansions (124) and (126) become nonuniform when = ( − 1 2 ) ± , respectively. Thus, we must introduce a transition region, when, as in Section 5, = − 1 2 = (1) ± as → ∞. We label this as region Tr. There are a number of subregions to consider. We begin first in regions ± . For brevity, we give essential details only. At leading order, we have simply, with = (1) ± , respectively, and Γ > 0 to be determined. There are now two cases to consider. First, we address the expansive case, when + > − . + > − : In this case, expansions (128) become nonuniform when = ± ± (1) ( + > − ), respectively, where two passive corner regions are required to locally connect with the main expansive region, in which ∈ ( − + (1), + − (1)) and = (1) as → ∞. Here, ± both depend upon ± and are determined interactively. As the corner regions are passive, we omit their details (and record that the corner regions will fix Γ in (128)) and proceed directly into the main expansive region, which we label as region . In region , we have, for ∈ ( − + (1) + , + − (1) + ) with ∶ ℝ → ℝ satisfying, Here, is the solution of the boundary value problem (EVP), namely, It is straightforward to establish (via a shooting and monotonicity argument) that (EVP) has a solution, which is unique. Moreover, it is readily shown that the solution has the following properties, F I G U R E 1 2 Profiles of ( , ) against from the numerical solution of (BIVP) with = 1∕2 and ( + , − ) = (2, −1) for = 100 − 1000 at intervals of 100 We can now match expansion (129) to the respective edge regions, which requires that + and − are the respective, unique, solutions to the two algebraic equations, We note that in this case, + > − , via (134). This completes region in the cases where 0 ∉ [ − , + ]. In the remaining cases, we must introduce a final localized region when = ( − 1 2 ) ( = (1)), with = (1), as → ∞. For convenience in what follows, we denote ′ (0) by > 0. In this final region, which we label as region loc, we obtain the expansion (omitting straightforward details), with = (1). The structure is now complete in the expansive case. We have seen that the principal feature is that of a decelerating expansion wave, now between the locations = Σ ± ( ), where, We complete this case with some numerical illustrations when = 1 2 . In Figure 12, we show the numerical solution to (BIVP) in the case ( + , − ) = (2, −1), while the antisymmetric case ( + , − ) = (1, −1) is illustrated in Figure 13. In each case, excellent agreement is shown F I G U R E 1 3 Profiles of ( , ) against from the numerical solution of (BIVP) with = 1∕2 and ( + , − ) = (1, −1) for = 100 − 1000 at intervals of 100 with the above asymptotic structure, with the principle feature in each case being an expansion wave on an ( 1 2 ) length scale. + < − ∶ In this compressive case, for brevity, we restrict attention to the situation when − = − + ≡ ∞ > 0, with the remaining cases commented on at the end of the subsection. In this case, it follows from Theorem 1 and Corollary 1, that the solution to (BIVP) is an odd, monotone decreasing, function of ∈ ℝ at each > 0, and thus we need only construct the structure in > 0. We find that there are two subcases to consider. The first is for ∞ > 1. We recall from (128) that in region + , when = (1) + , we have the expansion, which matches with region + , at leading order, as → ∞. At next order, we obtain the problem, with the requirements, which follows from Theorem 1, together with matching to region + , which requires, It is readily established (again via a shooting approach) that this problem has a one parameter family of solutions, parameterized by 0 > 0, and which has, and, with ∞ ( 0 ) > 0, depending upon 0 , and 0 to be determined interactively. We must now move on to region , a final region containing a stationary shock wave. In this region, = (1) + and = (1) − as → ∞, and thus we expand in the form, At leading order, we obtain the boundary value problem, with the latter condition coming from matching to the expansion in region + , which also requires, We label this boundary value problem as (SBVP). We consider (SBVP) via a shooting method from = 0. Thus, we consider the associated initial value problem consisting of the ODE (146), with initial conditions, with 0 > 0 fixed. This initial value problem has a unique global solution, which is monotone decreasing in ≥ 0, and has, with , ∶ ℝ + → ℝ being both continuous and positive, and being monotone. In addition, asymptotic solutions to the initial value problem when 0 is small and large, establish that, and We can now return to (SBVP), for which it now follows that a unique solution exists for each ∞ > 0, which corresponds to the solution of the above initial value problem with 0 = * , where * is the unique positive solution to the algebraic equation, after which condition (148) determines 0 as, Finally, we can perform detailed matching between region + and region + , which is achieved with the requirement that, via (144), This completes the structure when ∞ > 1.
The final case has 0 < ∞ ≤ 1. In this case, the stationary shock region is absent, and replaced by a weak and passive adjustment region when = (1). First, in region + , we now expand as, with = (1) + . The leading order problem is then given by the boundary value problem, We refer to this boundary value problem as (DBVP). It is straightforward to establish that the solution to (DBVP) is monotone decreasing in ≥ 0., and can be obtained by shooting from = 0, witĥ′(0) < 0 as the shooting parameter. A numerical investigation determines that for each 0 < ∞ ≤ 1, the boundary value problem (DBVP) has a unique solution, and in particular, this solution has,̂( and,̂( with the two positive constants * and ∞ depending continuously on ∞ , while matching to region + requires, with = (1) + . The leading order problem is then given by, with the second condition coming from matching with region + . The solution to this boundary value problem is readily obtained as, with > 0 a constant whose determination requires the construction, and matching to, higherorder terms in region + . For brevity, we do not pursue this further. The structure is now complete in this final case. When ∞ > 1, the key feature is a stationary smooth shock wave about = 0, while, when 0 < ∞ ≤ 1, this changes to a stationary diffusion scale front about = 0. We comment that a combination of these two key features arises in the non-antisymmetric cases depending on the values of ± in relation to the critical values ±1, respectively. We finish with some illustrative numerical solutions to (BIVP) in this case. In the top window of Figure 14, we present the numerical solution to (BIVP) when ∞ = 0.5.
As predicted by the asymptotic theory, the large-structure involves the formation of a antisymmetric diffusion front, spreading on a length scale of ( 1 2 ). Finally, in the bottom window of Figure 14, we present the numerical solution to (BIVP) when ∞ = 3 2 . In this case, as again predicted by the asymptotic theory, the key large-structure is a stationary smooth shock wave, on an (1) length scale.
All cases in the marginal situation, when = 1 2 , have now been dealt with. We conclude that the large− structure depends critically on + and − , and can involve as the key feature, either a smooth shock wave, an expansion wave, or a diffusion front.

DISCUSSION
This paper has considered the Riemann problem (BIVP), for a spatially inhomogeneous generalized Burgers equation, with spatial inhomogeneity characterized by the parameter > 0, with the spatial rate of decay of the nonlinear sound speed being algebraic with increasing distance, and increasing with increasing . The Riemann data parameters are ( + , − ), and need only be considered with ( + , − ) ∈ . The principal focus has been on constructing the detailed largeasymptotic structure of the solution to (BIVP), and thereby identifying the key attracting features in each case, as → ∞. There are three distinct cases. First, for weaker spatial decay rates, with ∈ (0, 1 2 ), the principal large-structure is an algebraically decelerating and weakly thickening, smooth shock wave when + < − , and conversely, a smooth decelerating and algebraically expanding expansion wave when + > − . The precise details are given in Section 4. For stronger spatial decay rates, with ∈ ( 1 2 , ∞), the principal large-structure switches to a diffusion dominated front on the diffusion length scale of ( 1 2 ), in all cases of the Riemann parameters. This is detailed in Section 5. Finally, the marginal case, when = 1 2 , has a principal large-structure which develops overall on a length scale of ( 1 2 ), and can be either a smooth decelerating and weakly thickening shock wave, an algebraically expanding expansion wave, or a diffusion front, depending crucially on the Riemann parameters. This case is presented in Section 6. Throughout, the asymptotic theory is supported by careful numerical solutions to (BIVP). More general qualitative theory for (BIVP) is given in the second of this pair of papers (MNII). To conclude, we observe that for modifications of (BIVP) with, specifically, more general "single hump," continuous and piecewise smooth forms for the spatial inhomogeneity function ℎ ∶ ℝ → [0, ∞), comparison with the present theory indicates that the approach developed here will again be applicable, with minor modifications. In particular, we anticipate that ℎ ( ) = (| | −1 ) as | | → ∞ is the critical decay rate, in the sense that faster decay rates inhibit smooth shock wave or expansion wave formation, with error function type diffusion-fronts forming when is large, while slower decay rates lead to decelerating and thickening smooth shock waves or decelerating expansion waves forming when is large. As determined here, the marginal case may involve any of these structures, depending on the Riemann parameters.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Zenodo at https://zenodo. org/record/6817204#.YsvqgHbMKUk , reference number 10.5281/zenodo.6817204. O R C I D John C. Meyer https://orcid.org/0000-0003-0708-1088