Nonlinear and Linearised Primal and Dual Initial Boundary Value Problems: When are they Bounded? How are they Connected?

Linearisation is often used as a first step in the analysis of nonlinear initial boundary value problems. The linearisation procedure frequently results in a confusing contradiction where the nonlinear problem conserves energy and has an energy bound but the linearised version does not (or vice versa). In this paper we attempt to resolve that contradiction and relate nonlinear energy conserving and bounded initial boundary value problems to their linearised versions and the related dual problems. We start by showing that a specific skew-symmetric form of the primal nonlinear problem leads to energy conservation and a bound. Next, we show that this specific form together with a non-standard linearisation procedure preserves these properties for the new slightly modified linearised problem. We proceed to show that the corresponding linear and nonlinear dual (or self-adjoint) problems also have bounds and conserve energy due to this specific formulation. Next, the implication of the new formulation on the choice of boundary conditions is discussed. A straightforward nonlinear and linear analysis may lead to a different number and type of boundary conditions required for an energy bound. We show that the new formulation shed some light on this contradiction. We conclude by illustrating that the new continuous formulation automatically lead to energy stable and energy conserving numerical approximations for both linear and nonlinear primal and dual problems if the approximations are formulated on summation-by-parts form.


Introduction
The classical way of analysing linear initial boundary value problems (IVBPs) uses the energy method [1,2,3,4,5] as the main tool. The energy method applied to linear IBVPs lead to well posed boundary conditions and energy estimates [6,7,8,9]. This procedure can also be used directly for certain types of nonlinear problems as recently shown in [10,11,12]. Another common procedure to obtain estimates for nonlinear problems is to use the entropy stability theory [13,14,15] which originated in [16,17,18,19,20,21] and is applied for example in [22,23,24,25,26].
In this paper we will not focus on nonlinear problems per se, but rather on the relation between a nonlinear problem, and its linearised counterpart, the variable coefficient problem. This relation is interesting for many reasons, perhaps the most important one (except for being an important analysis tool) is that it constitutes the starting point for the derivation of the dual (or adjoint) problem often used in optimisation [27,28,29,30,31,32]. Although we study energy estimates for nonlinear and linearised problems which eventually lead to energy stable schemes, our analysis may also shed some light on the linear stability problems recently observed for nonlinear entropy stable schemes [33,34]. A nonlinear problem connects to its linearised problem through the linearisation procedure which frequently results in confusing contradictions. Given that the nonlinear problem conserves energy and has an energy bound, the resulting linearised variable coefficient version often does not (or vice versa). In addition to being confusing, this leads to practical problems, namely: ill-posedness of the corresponding dual problems [35,36,37,38,39].
The observed ill-posedness of the dual problem typically manifests itself as a numerical stability problem (sometimes referred to as "butterfly effects" [40,41,42,43]) and various stabilisation techniques (often referring to the famous Lorenz 63 system [44]) have been developed, see for example [42,45,46,47,48,49,50]. Although we will not deal specifically with this issue, we note that the acclaimed Lorenz 63 system deals with ordinary differential equations, where the stabilising effect of correctly imposed boundary conditions in IBVPs is missing. The analysis in this paper indicate that some of these stabilisation techniques may not be required, and that both the linear and nonlinear dual problems have energy bounds if properly posed.
Nonlinear dual problems have received considerable interest recently. The general concept of nonlinear self-adjointness [51,52], which includes strict self-adjointness [51,53], quasi self-adjointness [54] and weak self-adjointness [55], was originally introduced to construct conservation laws associated with certain symmetry properties of differential equations [55,56,57]. Nonlinear self-adjointness has also been used to derive exact transformations between linear and nonlinear formulations [58]. We will also consider nonlinear self-adjoint problems, but our focus is different. We will consider the relation between nonlinear problems and related linear ones where the linear problems are obtained using the standard linearisation procedure which neglects quadratic terms [2,4,59].
Our ambition is to derive energy conserving and bounded formulations for both the nonlinear and linearised problem, later to be transformed to energy stable numerical approximations. We focus on the formulation of the continuous equations, but also discuss boundary conditions. A nonlinear and linear analysis may lead to a different number and type of boundary conditions required for an energy bound. This was shown in [60] for the nonlinear and linearised shallow water equations. There are also cases where this discrepancy does not occur, see [10,11,12,61] for examples. We show that the new skew-symmetric formulation presented herein shed some light on this issue.
Using lifting approaches [7,62] and proper continuous boundary conditions, it is straightforward to apply the results (formulations that lead to energy conservation and bounds) from the continuous analysis and develop stable numerical schemes. This procedure enables research groups using different numerical techniques such as finite difference [63,64,65], finite volume [66,67], spectral elements [68,69], flux reconstruction [70,71,71], discontinuous Galerkin [72,73,74] and continuous Galerkin schemes [75,76] to make use of the results. The only requirement is that one can formulate the numerical procedure on summation-by-parts (SBP) form with weak boundary conditions on simultaneous approximation term (SAT) form [77,78] or equivalently through numerical flux functions [74]. We conclude the paper by exemplifying how both linear and nonlinear stability follows almost automatically using the SBP-SAT technique with proper boundary conditions. The paper is organised as follows: We start by illustrating the confusing contradiction mentioned above for the scalar Burgers' equation in Section 2. This sets the stage for the discussion in Section 3 regarding which form a general, provably bounded nonlinear or linear hyperbolic IBVP should have. By exploiting a particular skew-symmetric form, we show how to relate the linearised problem to the nonlinear one in Section 4. Section 5 presents the nonlinear and linear dual problems and show that they also have energy bounds and preserve energy. In Section 6 we return to Burgers' equation and show how the initial contradiction discussed in Section 2 is resolved. In addition we discuss more complex examples involving systems of IBVPs. In Section 7, we discuss the relation between a linear and nonlinear boundary treatment. Section 8 illustrate the close relation between the skew-symmetric continuous formulation and stability and energy conservation of the numerical scheme. A summary and conclusions are provided in Section 9.

The standard linearisation procedure
We start by considering the one-dimensional (1D) hyperbolic (inviscid) Burgers' equation. The primitive, conservative and skew-symmetric forms are respectively. Here u(x, t) is the solution where x ∈ Ω = [a, b] and t ≥ 0 are the spatial and time coordinates. We also require that no discontinuities are present (which could be accomplished by adding suitable dissipative terms). With a smooth solution, all formulations in (2.1) are mathematically equivalent.

Energy bound in the nonlinear Burgers' equation
By applying the energy method (multiplying the equation by the solution and integrating over the domain), the formulations in (2.1) all lead to i.e. the energy u 2 2 = b a u 2 dx grows or decays only through boundary effects, which can be controlled by suitable boundary conditions. Note in particular that no volume terms contribute to the growth or decay.

Growth in the linearised Burgers' equation
Next we proceed in the standard way and make the ansatz u =ū(x, t) + u ′ (x, t) whereū is O(1) while u ′ is a small pertubation (|u ′ | ≪ |ū|). By inserting the ansatz into (2.1) we find In the standard linearisation procedure one cancels the last quadratic term due to it being negligible, and assumes [2,4,59] thatū holds. This procedure leads to two confusing results: the first equation in (2.4), which states thatū satisfies the original governing equation (2.1) can, strictly speaking, only hold if u ′ = 0 (disregarding the trivial case whereū is constant and assuming uniqueness). Applying the energy method to the second equation yields The volume term on the right-hand side may cause exponential growth (or decay) of energy, even if suitable boundary conditions are supplied. The variable coefficientū x produces the volume term. No such effect exist in the nonlinear problem. This initial example illustrates that even though the nonlinear problem has an energy bound and conserves energy as in (2.2), the linearised problem may not, as seen in (2.5).

Energy bounded and energy conserving linear and nonlinear primal problems
The results above for the Burgers' equation show that a general formulation that allows for an energy estimate for both the nonlinear and linearised equations would be of interest. Throughout this paper we restrict the analysis to the hyperbolic (inviscid) part of IBVPs, where the nonlinearity normally resides. The parabolic (viscous) part could be added on and properly posed provide dissipation or damping effects (although they can complicate the boundary treatment [2,3,6,8,10,11,12]).
Consider the following general hyperbolic IBVP augmented with homogeneous boundary conditions L p U = 0 at the boundary δΩ. In (3.1), the Einstein summation convention is used and P is a symmetric positive definite (or semi-definite) time-independent matrix that defines an energy norm (or semi-norm) U 2 P = Ω U T P U dΩ. We assume that U and V are smooth. The n × n matrices A i , B i , C are smooth functions (each matrix element is smooth) of the n component long vector V , but otherwise arbitrary. Note that (3.1) encapsulates both linear (V = U ) and nonlinear (V = U ) problems. The following two concepts are essential for a proper treatment of (3.1).
holds. It is energy bounded if it is energy conserving and the boundary conditions L p U = 0 are such that Proof. The energy method applied to (3.1) yields where (n 1 , .., n k ) T is the outward pointing unit normal. The right-hand side of (3.4) is cancelled by (3.2) leading to energy conservation. If in addition (3.3) holds, an energy bound is obtained.
Remark 3.2. The relation (3.2) leads a to skew-symmetric formulation, where symmetric matrices [5,7,8,79] are allowed, but not required for the energy method to work. Proposition 3.1 shows that whatever form the original IBVP has, energy boundedness and energy conservation can be proved if it can be rewritten in the form given by (3.1)-(3.2).
For later reference we introduce the linear (V = U ) and nonlinear (V = U ) primal problem to be analysed where F is a forcing function, g is boundary data and f initial data.
Remark 3.3. The skew-symmetric form in (3.5) leads to a complete derivative in the energy method since Remark 3.4. Inhomogeneous boundary conditions L p U = g will be discussed in Section 7.

A new non-standard linearisation procedure
In the first step, we divide the solution into a non-constant mean valueŪ and small perturbation U ′ as V = U =Ū + U ′ and insert this into the homogeneous version of (3.5). Next we introduce perturbed . This partition can always be done, since the matrix entries are smooth. The result is collects the nonlinear terms. The next step is interpreted in slightly different ways by different authors [2,4,59] as pointed out above. In the standard version one neglects H and requires that holds. However, this interpretation is doubtful sinceŪ can only solve (4.2) if U ′ vanishes as seen in (4.1). By carefully considering Proposition 3.1 we instead group the terms in the following way: A more plausible interpretation is thus obtained by neglecting H and demanding that which leaves the linearised primal equation to be Equations (4.5) and (4.6) are now both in the skew-symmetric form required in Proposition 3.1 for energy boundedness. This formulation removes the common confusing growth issue of the linearised problem, when the nonlinear original problem (3.5) has a bound. Note also that (4.5) and (4.6) are coupled.

Energy bounded and energy conserving linear and nonlinear dual problems
where G is a vector weight function. To derive the dual problem to (3.5), we search for the dual solution Φ such that J(U, G) = J(Φ, F ), where F is the forcing function in (3.5). Integration by parts yields With homogeneous initial condition for the primal problem, it follows that the dual end condition is Φ( x, T ) = 0. Furthermore, the dual homogeneous boundary conditions L d Φ = 0 must be such that Φ T (n i A i )U = 0, x ∈ ∂Ω when the primal homogeneous boundary conditions are applied. Finally, the dual equation becomes The derivation of the dual problem (5.1) holds both in the linear (V = U ) and nonlinear (V = U ) cases. By setting G = 0 and letting V = Φ we find that the dual problem is identical to the primal problem (3.5) with F = 0, and hence it is nonlinearly self-adjoint, or more precisely strictly self-adjoint [51,52,53].
Remark 5.1. The concept of nonlinear self-adjointness was originally introduced to aid construction of conservation laws associated with certain symmetry properties of differential equations. One may speculate (no proof exists) that many conservation laws could possibly be transformed to the skew-symmetric form in Proposition 3.1 and proven energy bounded. This provides interesting future research as many practical problems are formulated as conservation laws, but lack a proof of boundedness.
By introducing the time transformation τ = T − t, both the linear (V = Φ) and nonlinear (V = Φ) dual problems become where G is the forcing function, q is boundary data and r initial data. We state the corresponding results for the dual problem that hold for the primal one.
Proof. The proof is identical to the one for Proposition (3.1).

Examples
We will discuss four examples of increasing complexity. We start by again considering Burgers' equation.

The 1D scalar Burgers' equation
Any of the formulations in (2.1) can be used to arrive at the skew-symmetric (or self-adjoint) forms below, but perhaps starting with the confusing formulation in (2.3) is the most instructive way. We havē where h = u ′ u ′ x . By collecting terms based on the formulation (4.3) we find After neglecting the nonlinear term h we set both groups to zero and obtain the equations for the mean and the perturbation. This leads to skew-symmetric self-adjoint formulations for both the mean and the perturbation and the preceding theory in Sections 3, 4 and 5 applies.

The 2D incompressible Euler equations
The incompressible 2D Euler equations in primitive form with velocity field (u, v) in the (x, y) direction and pressure p (divided by the constant density) are 3) In a more compact formulation on matrix-vector form [10,11,12] we havẽ To arrive at the proper skew-symmetric form, we rewrite (6.4) using the splitting technique [80] as By inserting (6.6) into (6.4) and recalling that we are dealing with an incompressible fluid, i.e., A x + B y = (u x + v y )Ĩ = 0, we obtain the final form of the incompressible 2D Euler equations To connect with the general formulation (3.5): P =Ĩ, A 1 = A/2, A 2 = B/2, C = 0. Since the matrices A, B are already symmetric in the original formulation, the formulation (6.7) is in the proper skew-symmetric form, and the preceding theory in Sections 3, 4 and 5 can be applied.
Remark 6.1. We obtain an estimate in the semi-norm U 2 I = Ω U TĨ U dxdy involving only the velocities. Remark 6.2. The divergence relation can also be used to formulate (6.4) in strictly conservative form.

The 3D incompressible Euler equations in cylindrical coordinates
A slightly more complicated example is given by the 3D incompressible Euler equations in cylindrical coordinates [81] u t + u ∂u ∂r or in a compact (almost) conservative formulation with subscripts denoting differentiatioñ In (6.8),Ĩ = diag(1, 1, 1, 0), U = (u, v, w, p) ⊤ is the solution vector, where u, v, w are the velocities in the r, θ, z direction, respectively, and p is the pressure (divided by the constant density). Furthermore, we have Striving to use Propostion 3.1, we need the skew-symmetric form of (6.8). Again using the splitting technique [80] we rewrite the conservative derivative formulation using the generic formula and insert them into (6.8). The final result is where the derivatives on the matrices and R(U ) are grouped together to form To connect with the general formulation (3.5): multiply (6.9) with r and obtain P = rĨ, A 1 = rA/2, A 2 = B/2, A 3 = rC/2 while the zero order term C in (3.5) corresponds to D in (6.10). Since the matrices in (6.9) are symmetric (except D which is skew-symmetric as required), Proposition 3.1 applies.
Remark 6.3. The estimate is obtained in a semi-norm, in this case U 2 I = Ω U TĨ U rdrdθdz.

The 2D shallow water equations
The 2D shallow water equations (SWE) in primitive form [82] excluding the influence of bottom topography are where φ = gh [5] is the geopontential, h is the water height, g is the gravitational constant and u is the fluid velocity. The Coriolis forces are included with the function f which is typically a function of latitude [82,83]. Note that h > 0 and φ > 0 from physical considerations. In both Burgers' equation and the two forms of incompressible Euler equations discussed above, transformation of variables was not required in order to obtain the proper skew-symmetric form. However, for the SWE this must be done. The total energy ǫ is the sum of the kinetic and potential energy. By multiplying the total energy with the gravitational constant g we find This is an auxiliary conserved positive quantity for smooth solutions of (6.11) from which we choose our new variables U = (U 1 , U 2 , U 3 ) T = (φ, √ φu, √ φv)) T (each component now have the same physical units). To get an energy estimate we need governing equations for U . Repeated use of the chain rule in (6.11) yields Equation (6.13) is not in the skew-symmetric form (3.5) required for an estimate. To derive the skewsymmetric form, we first note that the different coordinate directions can be treated separately. We start with the x-direction and make the ansatz 16) and α, β, γ, θ, ψ are parameters to be determined. Equality in (6.16) leads to the one-parameter solution Remark 6.4. The 1D skew-symmetric version of (6.17) was first derived in [84].
Based on the result (6.17) in the x-direction and the structure of the original matrices (6.14) and (6.15) we make a guess (and verify) that the matrix A 2 in the y-direction is The skew-symmetric matrix C now denoted C in (6.15) remain the same after the transformation. The 2D SWE are now transformed to which has the required skew-symmetric form in Proposition 3.1. The rest of the theory presented in Sections 3, 4 and 5 follows provided that the correct number and type of boundary conditions are given.
Remark 6.5. The energy rate cannot depend on the free parameters α and β in the matrices A and B since they are not present in (6.13), (6.14) and (6.15). Hence as a sanity check we compute the boundary contraction (6.20) where for compactness we introduced the normal velocity u n = n 1 u + n 2 v. The dependency on the free parameters α and β vanishes. Foregoing the analysis on boundary conditions in the next section, we remark that these parameters in the boundary contraction do not vanish in the linearised case.
Remark 6.6. The estimate (3.4) with the same energy norm U 2 P , the same boundary term (6.20) and cancelled volume terms was also obtained in [60]. In that case the original equations (6.11), the variables U = (φ, u, v) and the norm matrix P = diag(1, φ, φ) was used. The time and space dependent norm matrix was obtained from a symmetrisation requirement, which the formulation (6.19) bypasses (such that P = I 3 can be used). This simplifies the production of a stable discrete approximation as discussed in Section 8.

Analysis and discussion of boundary conditions for nonlinear and linearised problems
A nonlinear and linear analysis may lead to a different number and type of boundary conditions required for an energy bound. This was shown in [60] where the boundary coefficient matrix n i A i were found to be distinctly different in the linear and nonlinear case. There are also cases where this discrepancy do not occur, see [10,11,12,61] for examples. We will discuss this issue in light of the examples and novel formulations discussed above.
Comparing the nonlinear skew-symmetric problem with the linearised one, i.e.
we find that they differ only in the arguments of the matrices. There are no terms missing, and hence the combined boundary coefficient matrix n i A i has the same structure in both cases. This suggests that the same procedure for deriving boundary conditions can be used in both the nonlinear and linearised case. In other words, the discrepancy when it comes to boundary conditions seen in [60] seems to vanish with the new linearisation formulation. This is, in fact, also the situation for the Burgers' equation and the two forms incompressible Euler equations discussed above. However, when considering the 2D SWE, the situation becomes more complicated due to the presence of the free parameters α and β. As shown in Remark 6.5, the dependence of α and β in the nonlinear boundary term U T (n i A i (U ))U is an illusion and cancel after contraction. However, the linear case is different because two sets of different variablesŪ and U ′ are involved, and the boundary term (U ′ ) T (n i A i (Ū ))U ′ retain the dependence on α and β. From (6.20) we found the parameter independent nonlinear boundary term to be while the parameter dependent linearised version is The formulation (7.2) holds for all α, β.
To shed some light on the difference between the linear and nonlinear boundary treatment, consider the special choice α = β = 1 which produces a boundary term similar to the nonlinear one Let us now compare (7.1) and (7.3). For inflow, when both u n andū n are negative, three boundary conditions seem to be needed in both cases while at outflow when bothū n and u n are positive, none is required. However, there is a significant difference between (7.1) and (7.3) since the diagonal entries in (7.1) are functions of the solution (u n = (n 1 U 2 + n 2 U 3 )/ √ U 1 ) while they can be considered as external data in the linear case asū n is independent of U ′ . Following [60] we aim for a minimal number of inflow boundary conditions and rewrite the nonlinear boundary term as where we introduced the normal U n = n 1 U 2 + n 2 U 3 = u n √ φ and tangential U τ = −n 2 U 2 + n 1 U 3 = u τ √ φ scaled velocities. Considering non-glancing boundaries we have min |U n | ≥ δ n > 0 and since φ > 0 we also have min | √ U 1 | ≥ δ 1 > 0. This leads to a bound using only two boundary conditions (recall that U n < 0) instead of three. By specifying U 2 n + U 2 1 = g 2 2 and U n U τ = g 2 3 where g 2 , g 3 are bounded functions, we get The reformulation of the boundary term in the nonlinear case from (7.1) to (7.4) was possible only because the diagonal entries in the boundary matrix were functions of the solution. In the linearised case, no such reformulation can be done.
Remark 7.1. In the nonlinear outflow case, the rewritten boundary term in (7.4) indicates that one boundary condition is required. However, this is an illusion. Aiming for a minimal number, the relation (7.1) shows that no boundary condition is required, precisely as in the linear case.
2. An energy estimate can always be obtained by specifying too many boundary conditions [6,7]. But, in general overprescribing the boundary conditions means that existence cannot be obtained for linear problems. The precise existence requirement is not known for nonlinear problems, but it is reasonable to assume (unless proven otherwise) that similar requirements as for linear problems hold.
Remark 7.3. The free parameters α and β in the 2D SWE appeared since we had to derive new evolution equations for the transformed variables forming the energy. New equations were not required for Burgers' equation and the incompressible Euler equations since they already had the required form in the original variables. The number of boundary conditions for the linear 2D SWE depend upon the parameters α and β which may or may not match the number required in the nonlinear case. One may speculate that nonlinear equations written on skew-symmetric form with the original variables have the same number and form of boundary conditions as the linearised version. However, if a transformation is required, caution must be taken since free parameters may appear (implicitly or explicitly), rendering the linear analysis unreliable for the nonlinear problem (or vice versa). This situation may have bearing on the recent results in [33,34].

A stable and energy conserving numerical approximation
To exemplify the straightforward construction of stable and energy conserving linear and nonlinear schemes based on the skew-symmetric formulation in Proposition 3.1, we consider a summation-by-parts (SBP) approximation of the following general 2D problem The boundary conditions (imposed through SAT terms or numerical flux functions) are assumed to be dissipative and ignored (we focus on the energy conserving properties of the numerical scheme). Equation (8.1) is semi-discretised in space using summation-by-parts operators as where P = P ⊗ I x ⊗ I y and U = ( U T 1 , U T 2 , ..., U T n ) T include approximations of U = (U 1 , U 2 , ..., U n ) T in each node. The matrix elements of A 1 , A 2 are matrices with node values of the matrix elements in A 1 , A 2 injected on the diagonal as exemplified below Moreover D x = I n ⊗ D x ⊗ I y and D y = I n ⊗ I x ⊗ D y where D x,y = P −1 x,y Q x,y are 1D SBP difference operators, P x,y are positive definite diagonal quadrature matrices, Q x,y satisfies the SBP constraint Q x,y + Q T x,y = B x,y = diag[−1, 0, ..., 0, 1], ⊗ denotes the Kronecker product and I with subscripts denote identity matrices. All matrices have appropriate sizes such that the matrix-matrix and matrix-vector operations are defined. Based on the 1D SBP operators, the 2D SBP relations mimicking multi-dimensional integration by parts are given by where U T B x V and U T B y V contain numerical integration along rectangular domain boundaries. In (8.4) we have usedP = I n ⊗ P x ⊗ P y , B x = (I n ⊗ B x ⊗ P y ) and B y = (I n ⊗ P x ⊗ B y ). The discrete energy method (multiply (8.2) from the left with U TP ) yields Next we introduce the notation U 2 PP = U T (PP) U , apply the SBP relations (8.4) and arrive at which mimics (3.4) perfectly (ignoring the zero order term). By using the fact thatPA 1 = A 1P and PA 2 = A 2P (since the matrices involved consist of diagonal blocks), the right hand side of (8.6) vanishes. Hence, the energy rate depends only on boundary terms and the scheme is energy conserving in the semidiscrete setting.
Remark 8.1. An energy bound and stability can be obtained by adding proper weak dissipative boundary conditions using the SAT technique or numerical flux functions.
Remark 8.2. It is irrelevant whether the matrices A 1 and A 2 are functions of the solution or not, i.e., if the problem is linear or nonlinear. The skew-symmetric formulation, a summation-by-parts discretisation and a proper weak boundary treatment are all that matters for stability and energy conservation.
Remark 8.3. The energy conservation for the linear and nonlinear primal problem illustrated above, holds also for the corresponding linear and nonlinear dual problems.

Summary and conclusions
The standard linearisation procedure often results in a confusing contradiction where the nonlinear problem conserves energy and has an energy bound but the linearised variable coefficient version does not. We have shown that a specific skew-symmetric form of the primal nonlinear problem leads to an energy bound and energy conservation. Next, it was shown that this skew-symmetric form together with a nonstandard linearisation procedure lead to an energy bounded and energy conserving formulation also of the new slightly modified linearised variable coefficient problem. We also showed that the corresponding linear and nonlinear dual (or self-adjoint) problems retain these properties due to this specific formulation.
The scalar Burgers' equation, the incompressible 2D Euler equations, the incompressible 3D Euler equations in cylindrical coordinates and the 2D shallow water equations were analysed to illustrate the new theory. From a study of these examples, we tentatively found that nonlinear equations on skew-symmetric form in the original variables have the same number and form of boundary conditions as the linearised version. However, if a variable transformation was required, caution must be taken, and the linear analysis may not be trustworthy. We concluded the paper by connecting the continuous analysis to a semi-discrete approximation. We showed that the skew-symmetric formulation automatically produced energy conserving and energy stable and numerical schemes for both linear and nonlinear primal and corresponding dual problems if these are formulated on summation-by-parts form.
The main finding in this paper can be summarised as follows. To obtain a formulation leading to an energy bound and energy conservation for general nonlinear hyperbolic IBVP, a number of actions are in general required. Firstly, one must choose dependent variables that forms an appropriate energy norm. Secondly, evolution equations for the new variables on skew-symmetric form must be derived. Finally, one must derive boundary conditions that limit the boundary terms. Some of the steps mentioned above can often be bypassed. In the examples treated, neither the Burgers' nor the incompressible Euler equations needed new variables and the skew-symmetric form followed directly. For a general nonlinear IBVP, choosing a set of transformation variables to produce a skew-symmetric formulation is likely the most difficult task.