November 2011 

NASA/TM–2011-217307 

Discretely Conservative Finite-Difference 

Formulations for Nonlinear Conservation Laws 

in Split Form: Theory and Boundary Conditions 

Travis C. Fisher and Mark H. Carpenter 

Langley Research Center, Hampton, Virginia 

Jan Nordström 

Linköping University, Linköping, Sweden 

Nail Yamaleev 

North Carolina A&T State University, Greensboro, North Carolina 

R. Charles Swanson 

Langley Research Center, Hampton, Virginia 
NASA STI Program . . . in Profile

 Since its founding, NASA has been dedicated to 

the advancement of aeronautics and space science. 

The NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

 The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. It 

collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI program 

provides access to the NASA Aeronautics and Space 

Database and its public interface, the NASA Technical 

Report Server, thus providing one of the largest 

collections of aeronautical and space science STI in 

the world. Results are published in both non-NASA 

channels and by NASA in the NASA STI Report 

Series, which includes the following report types: 

 TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but having 

less stringent limitations on manuscript length 

and extent of graphic presentations. 

 TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

 CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

 CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA. 

 SPECIAL PUBLICATION. Scientific, 

technical, or historical information from NASA 

programs, projects, and missions, often 

concerned with subjects having substantial 

public interest. 

 TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission. 

 Specialized services also include creating custom 

thesauri, building customized databases, and 

organizing and publishing research results. 

 For more information about the NASA STI 

program, see the following: 

 Access the NASA STI program home page at 

http://www.sti.nasa.gov

 E-mail your question via the Internet to 

help@sti.nasa.gov

 Fax your question to the NASA STI Help Desk 

at 443-757-5803 

 Phone the NASA STI Help Desk at 

443-757-5802 

 Write to: 

 NASA STI Help Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

November 2011 

NASA/TM–2011-217307 

Discretely Conservative Finite-Difference 

Formulations for Nonlinear Conservation Laws 

in Split Form: Theory and Boundary Conditions 

Travis C. Fisher and Mark H. Carpenter 

Langley Research Center, Hampton, Virginia 

Jan Nordström 

Linköping University, Linköping, Sweden 

Nail K. Yamaleev 

North Carolina A&T State University, Greensboro, North Carolina 

R. Charles Swanson 

Langley Research Center, Hampton, Virginia 
Acknowledgments

Prof. Nordstr ̈om contributed to this work while in residence at NASA Langley Research Center,

Technical monitor Joseph Morrison. The work of Prof. Yamaleev has been supported in part by

NASA under Grant NNX09AV08A and the Army Research Laboratory under Grant

W911NF-06-R-006. Dr. Swanson is a Distinguished Research Associate (DRA) at Langley.

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not

constitute an offical endorsement, either expressed or implied, of such products or manufacturers by the

National Aeronautics and Space Administration.

Available from:

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076-1320

443-757-5802
Abstract

Simulations of nonlinear conservation laws that admit discontinuous solutions are

typically restricted to discretizations of equations that are explicitly written in divergence form. This restriction is, however, unnecessary. Herein, linear combinations

of divergence and product rule forms that have been discretized using diagonalnorm skew-symmetric summation-by-parts (SBP) operators, are shown to satisfy

the sufficient conditions of the Lax-Wendroff theorem and thus are appropriate for

simulations of discontinuous physical phenomena. Furthermore, special treatments

are not required at the points that are near physical boundaries (i.e., discrete conservation is achieved throughout the entire computational domain, including the

boundaries). Examples are presented of a fourth-order, SBP finite-difference operator with second-order boundary closures. Sixth- and eighth-order constructions

are derived, and included in E. Narrow-stencil difference operators for linear viscous terms are also derived; these guarantee the conservative form of the combined

operator.

Contents

1 Introduction 2

2 Definitions 4

2.1 Summation-By-Parts Operators . . . . . . . . . . . . . . . . . . . . . 4

2.2 Complementary Grids . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.3 The Lax-Wendroff Theorem . . . . . . . . . . . . . . . . . . . . . . . 5

2.4 Split-Operator Consistency Requirements for Lax-Wendroff . . . . . 6

3 A Discretely Conservative, Finite-Domain, Split Form 6

3.1 A Fourth-Order Example . . . . . . . . . . . . . . . . . . . . . . . . 6

3.2 General SBP Operators . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.3 A Multidomain Approach for Alternative Operators . . . . . . . . . 10

3.3.1 Multidomain Operators . . . . . . . . . . . . . . . . . . . . . 10

3.3.2 Curvilinear Coordinates . . . . . . . . . . . . . . . . . . . . . 10

4 Examples using Burgers equation 11

4.1 Energy Analysis of the Continuous Problem . . . . . . . . . . . . . . 11

4.2 Energy Analysis of the Semidiscrete Problem: Single Domain . . . . 13

4.3 Energy Analysis of the Semidiscrete Problem: Multidomain . . . . . 15

5 Numerical Tests 17

5.1 Accuracy Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

5.2 Multidomain Validation . . . . . . . . . . . . . . . . . . . . . . . . . 20

6 Conclusion 21

A Discrete Conservation 26

1
B Split Continuous Equations 27

C High-Order Skew-Symmetric Fluxes 28

C.1 The Fourth-Order Operator: (3-6-3) . . . . . . . . . . . . . . . . . . 28

C.2 The Complete Fifth-Order Operator: (4-8-4) . . . . . . . . . . . . . 28

D Conservative Finite-Difference Viscous Terms 29

E Stencil Coefficients 31

E.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

E.2 Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

E.2.1 (2-4-2) Scheme Skew-Symmetric . . . . . . . . . . . . . . . . 32

E.2.2 (3-6-3) Scheme Skew-Symmetric . . . . . . . . . . . . . . . . 36

E.2.3 (4-8-4) Scheme Skew-Symmetric . . . . . . . . . . . . . . . . 47

1 Introduction

The nonlinear product terms in the continuous Euler equations can be expressed

in many forms: 1) conservative, ∂

∂xj

(uiuj ); 2) primitive, uj

∂ui

∂xj

; 3) skew-symmetric,

1

2

∂

∂xj

(uiuj ) + 1

2

uj

∂ui

∂xj

; 4) rotational, uj (

∂ui

∂xj

−

∂uj

∂xi

) −

∂

∂xj

(

1

2

uiuj ); as well as others. Although they are equivalent for smooth flows at the continuous level, each

form can exhibit profoundly different discrete accuracy, conservation, and nonlinear stability properties. Growing evidence suggests that the nonconservative forms,

particularly the skew-symmetric form, preserve important invariants of the continuous equations (e.g., kinetic energy or entropy) and deliver enhanced accuracy and

robustness. This makes the use of split forms an attractive alternative to conventional approaches based on the conservation form of the equations. Evidence that

supports this assertion is available in the following references [1–16].

Despite many favorable properties, alternative, nonconservative forms of the

Euler equations have historically been viewed as inappropriate for simulations that

involve discontinuities, such as supersonic and hypersonic flows. The basis of this

misconception likely stems from the Lax-Wendroff theorem [17] which stipulates

that the divergence form of the governing equations must be discretized with a

conservative operator, a condition that alternative discrete forms apparently lack.1

While alternative forms are generally not explicitly conservative, many can be

manipulated into an equivalent, consistent, and conservative form, thus allowing

the Lax-Wendroff theorem to be used to guarantee that the convergent captured

discontinuities are weak solutions of the governing equations. For example, Ducros

et al. [9] demonstrated that the skew-symmetric form of the convective terms in the

Euler equations are discretely conservative (at least for periodic fourth- and sixthorder centered operators), while Jameson [18] used the skew-symmetric form of the

Burgers equation and developed discretely conservative second order operators that

1Note that Gerritsen and Olssen [6] and, later, Yee and Vinokur [10] showed that correct shock

simulations could be achieved using the skew-symmetric form, even before a proof of discrete

conservation was identified.

2
are capable of capturing shocks. Recently Pirozzoli [15, 16] showed that all periodic

centered finite-difference operators can be expressed in a conservative form (even on

curvilinear grids) if an appropriate splitting of the equations is used.

Several open issues must be resolved before high-fidelity simulations that use

alternative forms of the equations become routine. One critical requirement is the

need for high-order boundary closures and boundary conditions that maintain the

desirable properties of the interior operators. Another is the need for a systematic

methodology to extend alternative single-domain operators (e.g., a skew-symmetric

discretization) to multiple domains, thereby extending the generality of any alternative approach.

A principle contribution of this work is the recognition that the boundary closures for alternative operators follow immediately if finite-domain summation-byparts (SBP) discrete operators are used to derive near-wall boundary closures. For

example, boundary closures that preserve the skew-symmetric form of the Euler

equations up to and including the boundary conditions are achieved by appropriately

splitting the equations and using SBP operators for all derivatives. We demonstrate

this point by developing three sets of fully conservative, alternative finite-difference

operators in this study: (2-4-2), (3-6-3) and (4-8-4).2

Another contribution of this work follows from addressing the broader question,

“Which discretizations of general conservation laws are discretely conservative?”

To this end, we conjecture and support with many examples that any combination of the conservative and primitive conservation law operators, αD(vw) + (1 −

α) (vDw + wDv), when discretized using any, diagonal-norm, skew-symmetric, SBP

discretization operator, can be expressed in a conservative form after the appropriate

discrete manipulations have been made. The generality of this conjecture allows for

quite arbitrary splittings of the equations while maintaining discrete conservation.

A final contribution of this work is the recognition that all single-domain alternative operators satisfy all of the conditions that are necessary for extension

to general, multiblock discretizations, provided that interface coupling terms that

maintain stability, accuracy, and conservation are used [21].

Another critical research element is the need for dissipation operators that are

complimentary with the skew-symmetric form, but maintain the design-order accuracy of the finite-domain operator. Rather than include this topic herein, we address

it in a companion paper. Because all of the discrete operators that are derived in this

paper are nondissipative and have no inherent mechanism to suppress oscillations,

we refrain from testing problems that admit strong shocks. As such, we validate the

accuracy of the new class of operators by using the viscous Burgers equation.

The layout of the paper is as follows. Section 2 defines SBP operators, briefly

summarizes the Lax-Wendroff theorem, and presents a derivation of the discretely

conservative split-form finite-difference approximations for conservation laws. Examples of a complete fourth-order operator are presented in section 3. Section 3 also

2The (p-2p-p) nomenclature denotes that boundary-interior-boundary stencils are p-, 2p-, porder accurate, respectively. The resulting operators are (p+1)-order accurate for hyperbolic equations (Euler) and incompletely parabolic equations (Navier-Stokes) and (p+2)-order accurate for

parabolic equations (viscous Burgers) [19, 20].

3
contains a discussion on the relationship between the SBP operators and split-form

conservation, as well as a discussion of how alternative operators can be extended

to complex geometries. In section 4, several split-form operators are applied to the

viscous Burgers equation. Numerical validation of accuracy of an alternative form

is given in section 5, with the use of both single and multiple domains. For brevity,

the (3-6-3) and (4-8-4) operators appended verbatim in E.

2 Definitions

2.1 Summation-By-Parts Operators

Discretize the physical domain by using N equidistant points that are distributed

over the interval [0, 1] as

x = (x1, x2, . . . , xN−1, xN ), xi = (i − 1) 1

N − 1

, i = 1, 2, . . . , N

and define a general class of SBP derivative operators Dsbp that satisfy the following

matrix properties:

Dsbp = P

−1Q, Q + QT = B = diag(−1, 0, · · · , 0, 1)

P = diag(p(1,1), · · · , p(N,N)), ζ

TPζ > 0, ζ 6= 0

(2.1)

where p(1,1) · · · p(N,N) are the diagonal elements of the matrix and ζ is an arbitrary

vector. The SBP operators defined in eq. 2.1 utilize a diagonal norm P and are,

therefore, a subset of the more general SBP form (see ref. [21] or ref. [22]). All are

telescoping operators and are conservative in the P norm. (See A for a proof that

all operators that can be expressed in the form of eq. 2.1 are conservative.)

2.2 Complementary Grids

The discrete derivative operator approximating hx (h sufficiently smooth) is represented herein as a combination of matrix operations that transfer data between two

complementary meshes. This nomenclature automatically leads to discrete conservation. The two complementary sets of grid points (or meshes) are defined on the

finite interval 0 ≤ x ≤ 1, differ in dimension by one and are expressed by using the

vectors

x = [x1, x2, · · · , xN−1, xN ]

T

, x ̄ = [ ̄x0, x ̄1, · · · , x ̄N−1, x ̄N ]

T

with x1 =  ̄x0 = 0 and xN =  ̄xN = 1. Herein, the points x are referred to as

the “solution points,” while the points x ̄ are referred to as the “flux points.” (The

overbar nomenclature is reserved herein for those quantities that are defined at the

flux points x ̄.)

Define the finite-domain conservative difference operator δh

δx in matrix notation

as

δh

δx = P

−1∆h ̄

4
where

∆ =





−1 1 0 0 0 0

0 −1 1 0 0 0

0 0 .

.

.

.

.

. 0 0

0 0 0 −1 1 0

0 0 0 0 −1 1





(2.2)

is a rectangular (N +1×N) matrix. The flux vector h ̄ constructed using appropriate

interpolation operators on the flux points x ̄, is differenced back onto the solution

points by ∆. The matrix P

−1 accounts for the appropriate grid spacings throughout the domain. Note that discrete conservation follows immediately because the

columns of ∆ sum to [−1, 0, · · · , 0, 1]. Further details on the definitions of the dual

grids, as well as implementation details can be found elsewhere [24].

2.3 The Lax-Wendroff Theorem

Consider a system of conservation laws of the form

∂v

∂t +

∂f(v)

∂x = 0, t ≥ 0, 0 ≤ x ≤ 1,

v(x, 0) = v0(x) , v(0, t) = bc(t)

(2.3)

where v = v(x, t) is the continuous solution vector, f = f(v) is a nonlinear flux vector, and v0(x) is a bounded piecewise continuous function in L2. Although smooth

solutions to eq. (2.3) do exist, no smoothness can be guaranteed due to the nonlinearity of f. Thus, we seek a weak solution to eq. (2.3) such that the integral

relation

Z

T

0

Z

1

0

[ψtv + ψxf(v)] dxdt =





Z

1

0

ψvdx





T

0

+





Z

T

0

ψf(v)dt





1

0

(2.4)

is satisfied for all smooth test functions ψ(x, t).

Choose a vector-valued function g(u) that is both a function of 2l compactly

supported vector arguments and related to the flux f with the sole requirement

that g(v, · · · , v) = f(v). Furthermore, define the discrete fluxes g to be of the form

g(x + δx/2) = g(u−l+1, · · · , ul), g(x − δx/2) = g(u−l, · · · , ul−1)

with suitably modified one-sided expressions near the physical boundaries.

With these definitions, we state here without proof the Lax-Wendroff theorem

[17] (originally derived for an infinite domain):

Theorem 2.1. Consider the conservative discrete analogue of eq. (2.3):

u(x, t + δt) − u(x, t)

δt +

g(x +

δx

2

) − g(x −

δx

2

)

δx = 0 (2.5)

If the discrete solution u(x, t) converges boundedly almost everywhere to u ̄(x, t), then

u ̄(x, t) is a weak solution to the continuous equation.

Finite-domain proofs of the Lax-Wendroff theorem are given in the context of

SBP operators in refs. [21, 23].

5
2.4 Split-Operator Consistency Requirements for Lax-Wendroff

Consider the flux f that is defined in eq. (2.3), which satisfies the property f(u) =

v(u)w(u). The general split flux which is composed of a linear combination of the

divergence and the product-rule components, can be written as

ut + αf(u)x + (1 − α) [v(u)w(u)x + w(u)v(u)x] = 0 (2.6)

Discretizing eq. (2.6) with an SBP spatial operator results in the following semidiscrete equation:

ut + αP

−1QWv + (1 − α)



V P

−1Qw + WP−1Qv



= 0

v = [v(u1), v(u2), . . . , v(uN )]T

, V = diag(v)

w = [w(u1), w(u2), . . . , w(uN )]T, W = diag(w)

(2.7)

Equation 2.7 does not explicitly lead to a discretely conservative flux form δg

δx .

Furthermore, if the resulting discrete flux can be manipulated into a conservation

flux that satisfies g(u, · · · , u) = f(u), then g is not immediately consistent with the

physical flux f = vw that satisfies the Rankine-Hugoniot relations. Thus, if the

Lax-Wendroff theorem is to be used to guarantee convergence to the weak solution

of the continuous equations, then the following conditions must be met:

1. The discrete spatial operator αP

−1QWv + (1 − α)



V P

−1Qw + WP−1Qv



must be equivalent to a telescoping form δg

δx .

2. The discrete flux g is Lax-consistent with the physical flux f, where f is the

flux that satisfies the Rankine-Hugoniot relations S[u]−[f] = 0. (The bracket

notation [] denotes the jump in the function.)

3 A Discretely Conservative, Finite-Domain, Split Form

3.1 A Fourth-Order Example

We begin by demonstrating that the fourth-order finite-domain split-form operator

is conservative for any value of the parameter α if a diagonal-norm, finite-domain

SBP operator is used for discrete differentiation.

The conservation part P

−1Q(Wv) at an interior point i (s < i < N + 1 − s)

is obviously conservative, and therefore does not need to be further manipulated.

(The number of closure terms that are needed at each boundary is given by the

constant s.) The product-rule terms V P

−1Qw + WP−1Qv, which are discretized

by using the conventional centered finite-difference at an interior point i, are

(V P

−1Qw + WP−1Qv)i =

vi

12

(wi−2 − 8wi−1 + 8wi+1 − wi+2)

+

wi

12

(vi−2 − 8vi−1 + 8vi+1 − vi+2)

(3.1)

6
and can be expressed in the following conservative form:

 ̄fi −  ̄fi−1 =

1

12





vi−2

vi−1

vi

vi+1

vi+2





T 



0 0 1 0 0

0 0 −8 0 0

1 −8 0 8 −1

0 0 8 0 0

0 0 −1 0 0









wi−2

wi−1

wi

wi+1

wi+2





(3.2)

The individual fluxes in eq. 3.2 are given as

 ̄fi =

1

12





vi−2

vi−1

vi

vi+1

vi+2





T 



0 0 0 0 0

0 0 0 −1 0

0 0 0 8 −1

0 −1 8 0 0

0 0 −1 0 0









wi−2

wi−1

wi

wi+1

wi+2





 ̄fi−1 =

1

12





vi−2

vi−1

vi

vi+1

vi+2





T 



0 0 −1 0 0

0 0 8 −1 0

−1 8 0 0 0

0 −1 0 0 0

0 0 0 0 0









wi−2

wi−1

wi

wi+1

wi+2





(3.3)

Clearly, the coefficients in the flux matrix simply shift down and to the right for

each subsequent flux  ̄fi (s < i < N + 1 − s). The boundary fluxes are more

elaborate but telescope in the same manner. The boundary fluxes are completely

determined by the form of Q near the boundaries for both the conservation and

product-rule forms.

The discrete split operator that is defined in eq. 2.7 can be written in conservative form by using a linear combination of the conservative divergence and

product-rule fluxes as

ut + P

−1∆ ̄f = 0,  ̄f = α ̄fc + (1 − α) ̄fe (3.4)

The divergence form flux is given by

 ̄fc =





w1v1

1

96 (48w1v1 + 59w2v2 − 8w3v3 − 3w4v4)

1

96 (−11w1v1 + 59w2v2 + 51w3v3 − 3w4v4)

1

96 (−3w1v1 + 51w3v3 + 56w4v4 − 8w5v5)

1

12 (−w3v3 + 7w4v4 + 7w5v5 − w6v6)

.

.

.

1

12 (−wi−1vi−1 + 7wivi + 7wi+1vi+1 − wi+2vi+2)

.

.

.

1

12 (−wN−5vN−5 + 7wN−4vN−4 + 7wN−3vN−3 − wN−2vN−2)

1

96 (−8uN−4vN−4 + 56uN−3vN−3 + 51uN−2vN−2 − 3uN vN )

1

96 (−3wN−3vN−3 + 51wN−2vN−2 + 59wN−1vN−1 − 11wN vN )

1

96 (−3wN−1vN−1 − 8wN−2vN−2 + 59wN−1vN−1 + 48wN vN )

wN vN





(3.5)

7
while the product-rule flux is

 ̄fe =





w1v1

1

96

(59u2v1 − 8u3v1 − 3u4v1 + 59u1v2 − 8u1v3 − 3u1v4)

1

96

(−8u3v1 − 3u4v1 + 59u3v2 − 8u1v3 + 59u2v3 − 3u1v4)

1

96

(−3u4v1 + 59u4v3 − 8u5v3 − 3u1v4 + 59u3v4 − 8u3v5)

1

12

(−w5v3 − w3v5 + 8w5v4 + 8w4v5 − w6v4 − w6v4)

.

.

.

1

12

( − wi−1vi+1 − wi+1vi−1 + 8wivi+1

+ 8wi+1vi − wivi+2 − wi+2vi )

.

.

.

1

12

( − wN−5vN−5 − wN−3vN−5 + 8wN−4vN−3

+ 8wN−3vN−4 − wN−3vN−1 − wN−1vN−3 )

1

96

( − 8uN−2vN−4 + 59uN−2vN−3 − 3uN vN−3

− 8uN−4vN−2 + 59uN−3vN−2 − 3uN−3vN )

1

96

( − 3uN vN−3 + 59uN−1vN−2 − 8uN vN−2

+ 59uN−2vN−1 − 3uN−3vN − 8uN−2vN )

1

96

( − 3uN vN−3 − 8uN vN−2 + 59uN vN−1

− 3uN−3vN − 8uN−2vN + 59uN−1vN )

wN vN





(3.6)

We also have constructed operators for 2p = 6 and 2p = 8 with boundary closure

blocks of dimension s = 7 and s = 9, respectively. The interior flux forms for the

sixth- and eight-order schemes are presented in C. The full flux form operators are

appended verbatim in E.

3.2 General SBP Operators

The previous fourth-order finite-domain example (as well as the sixth- and eighthorder examples, which are given in C) suggests that conventional SBP operators

lead to conservative finite-domain split-form operators. To this end, we make the

following conjecture:

Conjecture 3.1. The discrete split-form of eq. (2.7) can be manipulated into the

conservative form

ut + P

−1∆ ̄f = 0 (3.7)

for any diagonal-norm SBP operator that can be expressed in the form of eq. (2.1)

and for any value of the parameter α. Furthermore, the resulting local fluxes  ̄fj

have compact support and are consistent with the original conservative flux f(u).

8
The Lax-Wendroff theorem guarantees that discretely captured discontinuities are

accurate.

To prove this conjecture for a general banded matrix Q of halfwidth r, one would

need to first prove that the product rule term is conservative. That is,

(V Qw + WQv) = ∆ ̄fe

Because the term P

−1Qf in eq. (2.7) is already discretely conservative by eq. (A2)

(i.e., can be manipulated into the form P

−1∆ ̄fc), the split spatial operator is discretely conservative for any value of the parameter α.

Evidence: Although a formal proof is not forthcoming at this time, considerable

evidence suggests that the conjecture is true. Specifically, all of the banded highorder SBP finite-difference operators that are derived for this work (up to order

eight) are consistent with this conjecture.

Furthermore, consider a banded matrix Q that satisfies eq. (2.1), with elements

[Q](i1,i2) = q(i1,i2) and a halfwidth r ≥ 1. Solving for the conservative flux  ̄fe in the

relation (V Qw + WQv) = ∆ ̄fe [25] yields the closed-form expression

 ̄fj =

Xr

k=1

X

k

l=1

(wj+lvj+l−k + wj+l−kvj+l)q(j+l−k,j+l)

1 ≤ j + l, j + l − k ≤ N 1 ≤ j ≤ N − 1

 ̄f0 = w1v1, ̄fN = wN vN

(3.8)

Manipulating the q(j+l−k,j+l)terms produces the consistency condition

1

2

=

Xr

k=1

X

k

l=1

q(j+l−k,j+l)

1 ≤ j + l, j + l − k ≤ N 1 ≤ j ≤ N − 1

(3.9)

Substituting wj = w(u) and vj = v(u) into eq. (3.8) and using eq. (3.9) yields

 ̄fj = 2wvXr

k=1

X

k

l=1

q(j+l−k,j+l) = wv

1 ≤ j + l, j + l − k ≤ N 1 ≤ j ≤ N − 1,

 ̄f0 = wv,  ̄fN = wv

(3.10)

Thus, the term αQWv+ (1−α) (V Qw + WQv) can be expressed in terms of a telescoping conservative flux, has a compact support of halfwidth r, and is consistent

with the original conservative flux in the conservation law given by eq. (2.3). All

of the sufficient conditions of the Lax-Wendroff theorem are met, so converged solutions using the above split operators are weak solutions to the conservation law. 

Remark: This conjecture has been verified symbolically for matrices with a halfwidth

1 ≤ r ≤ 50.

9
3.3 A Multidomain Approach for Alternative Operators

While great flexibility is afforded by the potential use of alternative forms of the

equations (e.g., skew-symmetric or canonical splittings), a systematic approach is

needed to extend alternative uniform-grid single-domain, tensor-product operators

to complex geometries. A simple approach is to adopt a general multidomain finitedifference discretization that incorporates a generalized curvilinear mesh in each

domain.

3.3.1 Multidomain Operators

A major obstacle in the application of high-order methods to realistic problems is

the development a suitable grid around complex geometric features or in regions

of strong gradients. (Constraining the grid to design-order smoothness severely

complicates grid generation.) s equation is Multidomain techniques greatly simplify

the grid generation process for complex configurations by breaking the geometry into

the union of piecewise smooth quadrilateral (hexahedral) domains in two (three)

dimensions.

Conventional high-order SBP finite-difference techniques naturally extend to

multidomain discretizations [21, 23, 26, 27]. Each domain is discretized with a stable tensor-product formulation and then connected to its adjoining neighbors using

interface conditions that maintain the stability, accuracy, and conservation of the

interior operators. Penalty type interface treatments, which are closely related to

those that are used in discontinuous Galerkin and internal penalty finite-element

approaches, are most frequently used to join the domains. Domain interfaces need

only be C0 smooth to maintain the stability, conservation and design accuracy of

the single-domain operators.

Existing derivations for multidomain high-order finite-difference schemes primarily focus on the divergence form of the equations. Indeed, stability and conservation

proofs follow naturally with the use of this form. However, nothing precludes the

use of alternative operators in each subdomain. The only additional requirement is

the need for domain-interface coupling conditions that retain the desirable interior

properties of the alternative operator. These interface operators must be developed

on a case-by-case basis; however, as the examples in the next section demonstrate,

these operators generally are simple extensions of the original coupling conditions

introduced and extended in refs. [21, 23].

3.3.2 Curvilinear Coordinates

Generalized curvilinear coordinate formulations are commonly used in high-order

formulations to facilitate the use of nonuniform grid distributions or the definition

of “complex” geometries. The likelihood exists that the desirable attributes of a

Cartesian grid, alternative formulation will not survive the curvilinear mappings

between the physical (x, y, z) and the computational (ξ, η, ζ) space. However, the

curvilinear alternative formulation can inherit the desirable Cartesian grid properties. For example, Pirozzoli [16] recently demonstrated that pseudo-kinetic energy

10
can be preserved on curvilinear grids that use an alternative (skew-symmetric) formulation, (i.e., one motivated by formulations that preserve kinetic energy on a

Cartesian grid).

4 Examples using Burgers equation

We use the viscous Burgers equation to validate the numerical accuracy of the alternative split operators that are developed in the previous sections. In conservation

form, the Burgers equation is

ut + f(u)x = uxx, f(u) = 

u

2

2



, x ∈ [0, 1], t ∈ [0,∞)

a ̃0u(0, t) − ux(0, t) − g0(t) = 0, a ̃1u(1, t) − ux(1, t) − g1(t) = 0

a ̃0 =

u(0, t) + |u(0, t)|

3

, a ̃1 =

u(1, t) − |u(1, t)|

3

(4.1)

where  represents the diffusion coefficient. Although the problem is fundamentally

different with  = 0, we form the discrete approximation such that when  = 0

the fully inviscid Burgers equation is recovered and only the appropriate boundary

conditions are imposed. The boundary conditions that are given in eq. 4.1 ensure

that the formulation is strongly well-posed; this is demonstrated further in section

4.1.

The Burgers equation must be split to find continuous and discrete energy estimates. In the next section, we examine numerically the difference between the full

conservation form and the split-form by using the discretely conservative method

that is developed herein. While our formulation allows for any value of α, only one

α yields the energy estimate. This value of α is used in the next section to evaluate

advantages, if any, that are gained from using the discrete formulation that satisfies

the energy estimate.

4.1 Energy Analysis of the Continuous Problem

The properties and behavior of the inviscid portion of eq. (4.1) are covered in an

energy analysis that is presented by Jameson [18]; we reproduce and supplement the

essential results here for completeness. The energy that is given in Burgers equation

can be analyzed directly without a linearization. If we assume a smooth solution,

then can we split the conservation form (eq. (4.1)) by the parameter α to obtain

ut + α



u

2

2



x

+ (1 − α)uux = uxx, x ∈ [0, L], t ∈ [0,∞) (4.2)

11
This equation is multiplied by u and integrated over the interval [0, L] to obtain the

energy

1

2

d

dtkuk

2 = −

Z

L

0



αu 

u

2

2



x

+ (1 − α)u

2ux



dx + 

Z

L

0

uuxxdx

= −

Z

L

0



α

2



u

3



x

+



1 −

3

2

α



u

2ux



dx + 

Z

L

0

[(uux)x − uxux] dx

=

α

2

u

3







0

L

+ uux|

L

0 −

Z

L

0



1 −

3

2

α



u

2uxdx −  kuxk

2

(4.3)

Equation 4.3 clearly shows that for α = 2/3 the integral term disappears and the

energy is dependent on the boundary data and the viscous dissipation. This is referred to as the canonical splitting or “skew-symmetric form” of the Burgers equation [18, 28]. For the viscous problem ( > 0), note that any shock that develops is

resolvable at the continuous level and that the energy analysis remains valid. For the

inviscid problem ( = 0), the integration bounds can be split, and the same canonical splitting parameter is found. See reference [18] for details. The well posedness

of the boundary data is evaluated first by substituting α = 2/3 into eq. (4.3) and

then by analyzing each boundary individually. The left boundary term (scaled by

2 for convenience) is

BT0 =

2

3

u(0, t)

3 − 2u(0, t)ux(0, t) (4.4)

Completing the squares in eq. 4.4 leads to boundary terms that can be written as

BT0 =

1

2

3

u(0, t)

(

2

3

u(0, t)

2 − ux(0, t)

2

− (ux(0, t))2

)

(4.5)

We define

a(u) = 2

3

u(x, t) (4.6)

as a sensor for the hyperbolic part of the equation;  ̃a0 and  ̃a1 in the boundary

conditions of eq. (4.1) serve as switches to eliminate the hyperbolic part of the

boundary condition if waves are not propagating toward the interior of the domain

from the boundaries. If we use the definition above, the left boundary term is

BT0 =

1

a0

n

[a0u(0, t) − ux(0, t)]2 − [ux(0, t)]2

o

, a0 = a[u(0, t)] (4.7)

We now substitute the left boundary condition from eq. (4.1) for ux(0, t) in eq.

(4.7) to obtain

BT0 =

1

a0

n

[a0u(0, t) + g0(t) − a ̃0u(0, t)]2 − [ ̃a0u(0, t) − g0(t)]2

o

. (4.8)

12
The condition a0 > 0 in eq. (4.1) produces  ̃a0 = a, and we find that

BT0 =

g0(t)

2

a0

−

1

a

[a0u(0, t) − g0(t)]2 ≤

g0(t)

2

a0

(4.9)

If u(0, t) = a[u(0, t)] = 0, then BT0 = 0, and the energy of the continuous equation

is bounded from above by the boundary data. If a < 0, then  ̃a0 = 0, which yields

the condition

BT0 =

1

a0

[a0u(0, t) + g0(t)]2 −

g0(t)

2

a0

=

g0(t)

2

|a0|

−

1

|a0|

[a0u(0, t) + g0(t)]2(4.10)

Again, the energy is bounded from above by the boundary data. The right boundary

term (again scaled by 2) is

BT1 = −

2

3

u(L, t)

3 + 2u(L, t)ux(1, t) = −a1u(L, t)2 + 2 ̃a1 − 2u(L, t)g1(t) (4.11)

where a1 = a[u(L, t)]. The same procedure is followed to show that the energy on

the right boundary is bounded. We find that

BT1 =







0 for a1 = 0

g1(t)

2

a1

−

1

a1

[a1u(L, t) + g1(t)]2for a1 > 0

g1(t)

2

|a1| −

1

|a1|

[a1u(L, t) − g1(t)]2for a1 < 0

(4.12)

Since the energy of the continuous equation is bounded by the energy imposed

through the boundary terms, the equation and boundary conditions in eq. (4.1) are

strongly well-posed.

4.2 Energy Analysis of the Semidiscrete Problem: Single Domain

The canonical splitting can be used to construct a semidiscrete operator that satisfies

the semidiscrete analog to eq. (4.3). Using the finite-difference operators D and

D2 that satisfy the SBP condition on the finite domain, we discretize the skewsymmetric form (α = 2/3) of eq. (4.2) as

ut = −

1

3

[DUu + UDu] + D2u

+ σ0P

−1

e0 [ ̃a0u1 −  (Su)

1 − g0(t)]

+ σ1P

−1

e1 [ ̃a1uN −  (Su)N − g1(t)]

U = diag(u), e0 = (1, 0, . . .)

T

, e1 = (. . . , 0, 1)T

a ̃0 =

u1 + |u1|

3

a ̃1 =

uN − |uN |

3

(4.13)

The proper values of σ0 and σ1 are determined in the following energy analysis. Both

forms of D2 in eq. (D1) (see appendix D) satisfy the energy analysis, but we use

the narrow stencil here for illustrative purposes. The simultaneous approximation

term (SAT) penalty method [29] is used in eq. (4.13) to satisfy the mixed boundary

13
conditions with specified boundary data. For the finite domain, a diagonal norm

P is chosen to ensure that matrix multiplication of the norm and the nonlinear

coefficient matrix commutes. This choice reduces the overall accuracy for hyperbolic

problems to (p + 1), where 2p is the internal accuracy and p is the accuracy at the

boundary. For parabolic problems, the narrow-stencil second derivative operator

D2 = P

−1

(R + BS) constructed with the diagonal norm, provides a global (p + 2)

order of accuracy [30]. Premultiplying eq. (4.13) by u

TP and using eq. (2.1) yields

the discrete energy

1

2

d

dtkuk

2

P = −

1

3

u

T

(QU + UQ) u + u

T Ru + uTBSu

+ u

T

e0 {σ0 [ ̃a0u1 −  (Su)

1 − g0(t)]}

+ u

T

e1 {σ1 [ ̃a1uN −  (Su)N − g1(t)]}

(4.14)

Adding eq. (4.14) to its transpose yields

d

dtkuk

2

P = −

2

3

u

TBUu + 2uT Ru + 2uTBSu

+ 2u1σ0 ( ̃a0u1 −  [Su)

1 − g0(t)]

+ 2uN σ1 ( ̃a1uN −  [Su)N − g1(t)]

(4.15)

Simplifying with eq. (2.1) yields

d

dtkuk

2

P = 2u

T Ru

+ 2u1



σ0 [ ̃a0u1 − g0(t)] + 1

3

u

2

1



− 2u1 [(σ0 + 1) (Su)

1

]

+ 2uN



σ1 [ ̃a1uN − g1(t)] −

1

3

u

2

N



+ 2uN [(−σ1 + 1) (Su)N ]

(4.16)

Above, R is negative semidefinite and ensures that the energy only decays. The

viscous boundary terms cancel for σ0 = −1 and σ1 = 1, which leaves

d

dtkuk

2

P = 2u

T Ru + BTg0 + BTg1

BTg0 = u1 [a(u1)u1 − 2 ̃a0u1 + 2g0(t)]

BTg1 = uN [2 ̃a1uN − a(uN )uN − 2g1(t)]

(4.17)

For the discrete equation to be strongly well-posed, the discrete energy must be

bounded from above by the imposed boundary data. For the left boundary when

a(u1) ≥ 0 and using completing the squares,

BTg0 = −au2

1 + 2u1g0 =

g

2

0

a

−

1

a

(au1 − g0)

2

. (4.18)

The result in eq. 4.18 is equivalent to continuous result found in eq. (4.9).

14
Completing the squares at the left boundary if a(u1) < 0 yields

BTg0 = au2

1 + 2u1g0 =

g

2

0

|a|

−

1

|a|

(au1 + g0)

2

(4.19)

Again, note the similarity between eq. 4.19 and the continuous result given in eq.

(4.10). The right boundary also mimics the continuous case, so the semidiscrete

equation is strongly well-posed. Note that the same penalty is used for the inviscid

and the viscous Burgers equations and that the equation is well-posedness in both

cases.

4.3 Energy Analysis of the Semidiscrete Problem: Multidomain

Often multiple semidiscrete domains can be used advantageously to discretize a

conservation law. Interface penalties are derived such that the semidiscrete energy

is unaffected by the exchange of information between two domains. For a two domain

problem on the interval x ∈ [−1, 1] divided at x = b, we use the discretization

x

(1) =

h

x

(1)

1

, x

(1)

2

, . . . , x

(1)

N1−1

, x

(1)

N1

i

, x

(1)

i = −1 + (i − 1) b + 1

N1 − 1

, i = 1, 2, . . . , N1

x

(2) =

h

x

(2)

1

, x

(2)

2

, . . . , x

(2)

N2−1

, x

(2)

N2

i

, x

(2)

i = b + (i − 1) 1 − b

N2 − 1

, i = 1, 2, . . . , N2

where N1 and N2 are the number of uniform cells in each of the two domains. The

solution on each domain is denoted by u

(1) and u(2), respectively. Operators that are

defined on each domain are denoted similarly. The semidiscrete Burgers equation

for the multidomain problem is

u

(1)

t = −

1

3



D

(1)U(1)u(1) + U(1)D(1)u(1)

+ D

(1)

2 u

(1)

+



P

(1)−1

e

(1)

1



k00u

(1)

N1

− k01u

(2)

1



+ l01 

S

(1)u(1)

N1

−



S

(2)u(2)

1



+



P

(1)−1 

S

(1)T

e

(1)

1

l10 

u

(1)

N1

− u

(2)

1



u

(2)

t = −

1

3



D

(2)U(2)u(2) + U(2)D(2)u(2)

+ D

(2)

2 u

(2)

+



P

(2)−1

e

(2)

0



k11u

(2)

1 − k10u

(1)

N1



+ r01 

S

(2)u(2)

1

−



S

(1)u(1)

N1



+



P

(2)−1 

S

(2)T

e

(2)

0

r10 

u

(2)

1 − u

(1)

N1



U

(1) = diag(u(1)), U(2) = diag(u(2)), e

(1)

1 = (. . . , 0, 1)T

, e

(2)

0 = (1, 0, . . .)

T

(4.20)

where the boundary condition penalties have been dropped for clarity of the interface treatment. The energy is calculated by premultiplying the two equations by



u

(1)T

P

(1) and 

u

(2)T

P

(2), respectively. The equations are then added to their

15
transpose as

d

dt 








u

(1)










2

P(1)

+










u

(2)










2

P(2)

+ 2










u

(1)

x










2

P(1)

+










u

(2)

x










2

P(2)

=

2

3



u

(2)

1

3

−



u

(1)

N1

3



+ 2



S

(1)u(1)

N1

−



S

(2)u(2)

1



+ 2u

(1)

N1



k00u

(1)

N1

− k01u

(2)

1



+ l01 S

(1)u(1)

N1

−



S

(2)u(2)

1



+ 2u

(2)

1



k11u

(2)

1 − k10u

(1)

N1



+ r01 S

(2)u(2)

1

−



S

(1)u(1)

N1



+ 2 

S

(1)u(1)

N1

l10 

u

(1)

N1

− u

(2)

1



+ 2 S

(2)u(2)

1

r10 u

(2)

1 − u

(1)

N1



(4.21)

The terms that are related to the boundary data have been removed. Because we

want to construct the interface penalty such that it conserves the energy, we must

choose the free parameters such that the right-hand side of eq. (4.21) vanishes. This

results in the following constraints:

k00 =

1

3

u

(1)

N1

, k10 = −k01, k11 = −

1

3

u

(2)

1

l10 = −l01 − 1, r01 = l01 + 1, r10 = −l01

The formal order of accuracy of the penalty terms must be the same as the order of

the boundary finite-difference approximations. This requires that the value of k01

be a function of u. We choose an equal weighting of the values on either side of the

interface:

k01 =

u

(1)

N1

+ u

(2)

1

2

(4.22)

For the viscous terms, l01 = −

1

2

yields an equal weighting of all of the viscous

penalty terms; thus, this is the chosen value. The full semidiscrete set of equations

16
for the two-domain problem is

u

(1)

t = −

1

3



D

(1)U(1)u(1) + U(1)D(1)u(1)

+ D

(1)

2 u

(1)

−



P

(1)−1

e

(1)

0

h

a ̃

(1)

0

u

(1)

1 − 



S

(1)u(1)

1

− g0(t)

i

+



P

(1)−1

e

(1)

1

( u

(1)

N1

u

(1)

N1

−

u

(1)

N1

+ u

(2)

1

2

u

(2)

1

!

−



2



S

(1)u(1)

N1

−



S

(2)u(2)

1



−



2



P

(1)−1 

S

(1)T

e

(1)

1



u

(1)

N1

− u

(2)

1



,

u

(2)

t = −

1

3



D

(2)U(2)u(2) + U(2)D(2)u(2)

+ D

(2)

2 u

(2)

−



P

(2)−1

e

(2)

0

( u

(2)

1

u

(2)

1 −

u

(1)

N1

+ u

(2)

1

2

u

(1)

N1

!

−



2



S

(2)u(2)

1

−



S

(1)u(1)

N1



+



2



P

(2)−1 

S

(2)T

e

(2)

0



u

(2)

1 − u

(1)

N1



+



P

(2)−1

e

(2)

1



a ̃

(2)

1

u

(2)

N2

− 



S

(2)u(2)

N2

− g1(t)



.

(4.23)

In this form, the boundary penalty does not contribute to the energy. The energy changes, however, as a result of the one-sided stencils at the interface. The

semidiscrete equations remain strongly well-posed.

5 Numerical Tests

The new discrete operators were implemented for the Burgers equation in a onedimensional finite-difference solver with a uniform grid. Integration in time was

conducted with a five-step, fourth-order Runge-Kutta scheme [31]. The step size in

time was chosen such that the temporal error is negligible in comparison with the

spatial error.

Two problems are evaluated. For the first problem, the accuracy of the three

operator sets that are developed herein is tested by using the skew-symmetric (α =

2/3) and the conservation (α = 1) forms of the viscous Burgers equation. The

second problem tests the accuracy of the multidomain interface closures that are

developed based on the new operator sets.

17
5.1 Accuracy Validation

The viscous Burgers equation is used to test the accuracy of the new operator when

it is applied to the nonlinear equation. The test problem is defined as

ut = P

−1∆



α ̄fc + (1 − α) ̄fe +  ̄fV



+ σ0P

−1

e0 [−a ̃0u1 −  (Su)

1 − g0(t)]

+ σ1P

−1

e1 [−a ̃1uN −  (Su)N − g1(t)]

g0(t) = v(−1, t) + |v(−1, t)|

3

v(−1, t) − vx(−1, t)

g1(t) = v(1, t) − |v(1, t)|

3

v(1, t) − vx(1, t)

v(x, t) = 4x

x

2 + 2t +

1

40

, x ∈ [x1, xN ], t ∈ [0, T]

(5.1)

where  = 1 and the initial condition is u0 = v(x, 0). The definition of the viscous

flux  ̄fV is given in D. The values for the penalty coefficients σ are those that are

identified in section 4. Note that the sign of the convective term has changed for

this problem. This is to satisfy the exact solution v(x, t). The problem is simulated

on the interval [−1, 1] up to T = 0.05. Five grids are used to determine the order

of accuracy of the set of operators. The L

2 and L∞ error norms of the conservation

form and the canonical split-form are compared in tables 1 through 3. Figure 1

shows the numerical solution for N = 64 at different times. All simulations use the

narrow-stencil viscous operator.

Table 1. Error Norms and Convergence Rates for Viscous Burgers Equation Using

Canonical Split and Conservation Forms, and Conservative Narrow-stencil (2-4-2)

Viscous Operator.

Canonical Split Conservation

N L2 Error L2 Rate L∞ Error L∞ Rate L2 Error L2 Rate L∞ Error L∞ Rate

32 1.64e-03 - 2.70e-03 - 1.29e-03 - 2.44e-03 -

64 1.10e-04 3.90 1.91e-04 3.82 8.86e-05 3.86 1.75e-04 3.80

128 6.97e-06 3.98 1.21e-05 3.97 5.66e-06 3.97 1.17e-05 3.90

256 4.38e-07 3.99 7.66e-07 3.99 3.56e-07 3.99 7.38e-07 3.99

512 2.74e-08 4.00 4.79e-08 4.00 2.23e-08 4.00 4.62e-08 4.00

Table 2. Error Norms and Convergence Rates for Viscous Burgers Equation Using

Canonical Split and Conservation Forms, and Conservative Narrow-stencil (3-6-3)

Viscous Operator.

Canonical Split Conservation

N L2 Error L2 Rate L∞ Error L∞ Rate L2 Error L2 Rate L∞ Error L∞ Rate

32 2.34e-04 - 4.62e-04 - 2.83e-04 - 6.24e-04 -

64 4.56e-06 5.69 9.15e-06 5.66 5.61e-06 5.65 1.30e-05 5.58

128 7.52e-08 5.92 1.51e-07 5.92 9.31e-08 5.91 2.19e-07 5.90

256 1.19e-09 5.98 2.40e-09 5.98 1.48e-09 5.98 3.49e-09 5.97

512 2.00e-11 5.90 3.98e-11 5.92 2.33e-11 5.99 5.37e-11 6.02

The expected order of accuracy is p + 2, where p is the boundary accuracy.

This design order is achieved for both the skew-symmetric and conservation forms

18
Table 3. Error Norms and Convergence Rates for Viscous Burgers Equation Using

Canonical Split and Conservation Forms, and Conservative Narrow-stencil (4-8-4)

Viscous Operator.

Canonical Split Conservation

N L2 Error L2 Rate L∞ Error L∞ Rate L2 Error L2 Rate L∞ Error L∞ Rate

32 1.52e-04 - 3.64e-04 - 1.89e-04 - 4.42e-04 -

64 6.97e-07 7.77 2.17e-06 7.39 8.58e-07 7.78 2.46e-06 7.49

128 3.85e-09 7.50 1.78e-08 6.93 4.45e-09 7.59 1.92e-08 7.00

256 2.74e-11 7.13 1.91e-10 6.55 2.94e-11 7.24 1.98e-10 6.60

512 1.74e-12 3.98 2.65e-12 6.17 1.73e-12 4.09 2.69e-12 6.20

x

u(x)

-1 -0.5 0 0.5 1

-10

-5

0

5

10

t = 0.00

t = 0.01

t = 0.02

t = 0.03

t = 0.04

t = 0.05

Figure 1. Numerical approximation to skew-symmetric form of viscous Burgers

equation with 64 uniform cells and the (4-8-4) narrow stencil conservative operators.

for the (2-4-2) operator. However, the (3-6-3) and (4-8-4) operators exhibit greater

accuracy than design order convergence. The (2-4-2) operator does not deviate from

design order because the maximum order that can be achieved is equivalent to the

design order. This is not the case for the (3-6-3) and (4-8-4) operators. For these

operators, the accuracy is at least design order for well resolved solutions. For the

(2-4-2) operator, the conservation form is slightly more accurate than the canonical

split-form. For the (3-6-3) and (4-8-4) operators, the canonical split-form is more

accurate.

19
5.2 Multidomain Validation

The test problem that was used to test the interfaces is described by the exact

solution

v(x, t) = 1 − tanh 

x − t − x0

2



, x ∈ [−10, 10], t ∈ [0, T] (5.2)

where  = 0.25, x0 = −5.0, and the initial condition is u0 = v(x, 0). The problem is

simulated up to T = 10.0 with four grid resolutions. The skew-symmetric form of

the discrete equations is used.

The L

2 and L∞ errors are tabulated in tables 4 through 6 for each operator on

(1) a single domain, (2) two domains with equivalent spacing, and (3) two domains

with grid spacing that changes discontinuously at the interface by a factor of r = 2.

(Herein, the jump in spacing is denoted as the “compression ratio.”) The interface

between domains is located at x = 0.0. The discrete solution with 256 uniform cells

on a single domain for the (2-4-2) operator set is shown in figure 2.

x

u(x)

-10 -5 0 5 10 0

0.5

1

1.5

2

t = 0.0

t = 2.5

t = 5.0

Figure 2. Discrete solution to viscous Burgers equation using 256 equispaced nodes

with (4-8-4) narrow operators.

Tables 4 through 6 demonstrate that the interface has a negligible effect on the

accuracy for this well-resolved test case. In addition, the treatment has been verified

to be design order accurate, even when a discontinuous grid resolution is used across

the subdomains.

20
Table 4. Error Norms and Convergence Rates of (2-4-2) Operator for Single Domain,

Two Domains with Uniform Resolution, and Two Domains with compression ratio

r = 2.

Single Domain

Resolution Error

N L2 Error L2 Rate L∞ Error L∞ Rate

64 8.89e-03 - 1.13e-02 -

128 6.37e-04 3.80 8.86e-04 3.68

256 4.19e-05 3.93 6.17e-05 3.84

512 2.66e-06 3.98 4.04e-06 3.93

Two Domains, r = 1

Resolution Error

N1 N2 L2 Error L2 Rate L∞ Error L∞ Rate

32 32 8.91e-03 - 1.14e-02 -

64 64 6.38e-04 3.80 8.87e-04 3.68

128 128 4.19e-05 3.93 6.17e-05 3.85

256 256 2.66e-06 3.98 4.04e-06 3.93

Two Domains, r = 2

Resolution Error

N1 N2 L2 Error L2 Rate L∞ Error L∞ Rate

32 64 6.26e-04 - 8.65e-04 -

64 128 4.12e-05 3.93 6.08e-05 3.83

128 256 2.61e-06 3.98 3.97e-06 3.94

256 512 1.64e-07 4.00 2.50e-07 3.99

Table 5. Error Norms and Convergence Rates of (3-6-3) Operator for Single Domain,

Two Domains with Uniform Resolution, and Two Domains with compression ratio

r = 2.

Single Domain

Resolution Error

N L2 Error L2 Rate L∞ Error L∞ Rate

64 3.48e-03 - 3.96e-03 -

128 9.88e-05 5.14 1.33e-04 4.89

256 1.89e-06 5.71 3.23e-06 5.37

512 3.13e-08 5.92 5.41e-08 5.90

Two Domains, r = 1

Resolution Error

N1 N2 L2 Error L2 Rate L∞ Error L∞ Rate

32 32 3.47e-03 - 3.91e-03 -

64 64 9.87e-05 5.14 1.33e-04 4.88

128 128 1.89e-06 5.71 3.21e-06 5.37

256 256 3.11e-08 5.92 5.37e-08 5.90

Two Domains, r = 2

Resolution Error

N1 N2 L2 Error L2 Rate L∞ Error L∞ Rate

32 64 9.82e-05 - 1.34e-04 -

64 128 1.87e-06 5.72 3.15e-06 5.41

128 256 3.11e-08 5.91 5.05e-08 5.96

256 512 5.54e-10 5.81 7.55e-10 6.07

6 Conclusion

We have developed a class of split-form finite-difference operators that satisfy the

sufficient conditions of the Lax-Wendroff theorem. These operators are applicable to

the conservation or divergence form of the conservation law, but facilitate operator

21
Table 6. Error Norms and Convergence Rates of (4-8-4) Operator for Single Domain,

Two Domains with Uniform Resolution, and Two Domains with compression ratio

r = 2.

Single Domain

Resolution Error

N L2 Error L2 Rate L∞ Error L∞ Rate

64 1.97e-03 - 2.00e-03 -

128 2.88e-05 6.10 3.56e-05 5.81

256 1.81e-07 7.31 3.37e-07 6.72

512 8.53e-10 7.73 1.36e-09 7.95

Two Domains, r = 1

Resolution Error

N1 N2 L2 Error L2 Rate L∞ Error L∞ Rate

32 32 2.16e-03 - 2.70e-03 -

64 64 2.89e-05 6.23 3.79e-05 6.15

128 128 1.84e-07 7.29 2.93e-07 7.02

256 256 1.01e-09 7.51 1.67e-09 7.46

Two Domains, r = 2

Resolution Error

N1 N2 L2 Error L2 Rate L∞ Error L∞ Rate

32 64 3.01e-05 - 4.44e-05 -

64 128 2.67e-07 6.82 4.36e-07 6.67

128 256 3.43e-09 6.28 4.63e-09 6.56

256 512 1.09e-10 4.98 4.83e-11 6.58

splitting at the discrete level to improve accuracy or robustness, depending on the

conservation law. The method for constructing a conservative split operator is

illustrated for the (2-4-2) operator with a general splitting parameter. Higher order

operators were also derived and are supplied in an accompanying text file. We also

developed fully conservative operators for linear viscous terms.

The specific operators that are derived herein were tested on the conservation and

canonical split forms of Burgers equation. Energy stability of the discrete form was

proven, including boundary and interface closures. The split and conservation forms

yield very similar results for Burgers equation. The solutions converged at least at

the design order of accuracy in all cases, even when multiple domains were used

with discontinuous grid resolutions across interfaces. The interface treatment that

was tested herein had a negligible effect on the accuracy of the solution. Nonlinear

dissipation operators for split forms will be presented in a future paper.

References

1. A. Arakawa, “Computational design for long-term numerical integration of the

equations of fluid motion: Two-dimensional incompressible flow. part i.,” Journal of Computational Physics, vol. 1, pp. 119–143, 1966.

2. W. J. Feiereisen, W. C. Reynolds, and J. H. Ferziger, “Numerical simulation of

a compressible homogeneous, turbulent shear flow,” Tech. Rep. TF-13, Stanford

University, 1981.

22
3. A. Harten, “On the symmetric form of systems of conservation laws with entropy,” Journal of Computational Physics, vol. 49, pp. 151–164, 1983.

4. E. Tadmor, “Skew-selfadjoint form of systems of conservation laws,” Journal of

Mathematical Analysis and Applications, vol. 103, pp. 428–442, 1984.

5. T. Zang, “On the rotation and skew-symmetric forms for incompressible flow

simulations,” Applied Numerical Mathematics, vol. 7, pp. 27–40, 1991.

6. M. Gerritsen and P. Olsson, “Designing an efficient solution strategy for fluid

flows,” Journal of Computational Physics, vol. 129, pp. 245–262, 1996.

7. G. A. Blaisdell, E. T. Spyropoulos, and J. H. Qin, “The effect of the formulation

of nonlinear terms on aliasing errors in spectral methods,” Applied Numerical

Mathematics, vol. 21, pp. 207–219, 1996.

8. Y. Morinishi, T. S. Lund, O. V. Vasilyev, and P. Moin, “Fully conservative

higher order finite difference schemes for incompressible flow,” Journal of Computational Physics, vol. 143, pp. 90–124, 1998.

9. F. Ducros, F. Laporte, T. Soul`eres, V. Guinot, P. Moinat, and B. Caruelle,

“High-order fluxes for conservative skew-symmetric-like schemes in structured

meshes: Applications to compressible flows,” Journal of Computational Physics,

vol. 161, pp. 114–139, 2000.

10. H. C. Yee, M. Vinokur, and M. J. Djomehri, “Entropy splitting and numerical

dissipation,” Journal of Computational Physics, vol. 162, pp. 33–81, 2000.

11. H. Yee and B. Sj ́ogreen, Designing Adaptive Low-Dissipative High Order

Schemes for Long-Time Integrations, ch. One: Turbulent Flow Computation.

Kluwer Academic Publisher, 2001.

12. A. E. Honein and P. Moin, “Higher entropy conservation and numerical stability of compressible turbulence simulations,” Journal of Computational Physics,

vol. 201, pp. 531–545, 2004.

13. C. Kennedy and A. Gruber, “Reducing aliasing formulations of the convective

terms within the navier-stokes equations for a compressible fluid,” Journal of

Computational Physics, vol. 227(3), pp. 1676–1700, 2006.

14. P. K. Subbareddy and G. V. Candler, “A fully discrete, kinetic energy consisten

finite-volume scheme for compressible flows,” Journal of Computational Physics,

vol. 228(5), pp. 1347–1364, 2009.

15. S. Pirozzoli, “Generalized conservative approximations of split convective derivative operators,” Journal of Computational Physics, vol. 229, pp. 7180–7190,

2010.

16. S. Pirozzoli, “Stabilized non-dissipative approximations of Euler equations

in generalized curvilinear coordinates,” Journal of Computational Physics,

vol. 230, pp. 2997–3014, 2011.

23
17. P. Lax and B. Wendroff, “Systems of conservation laws,” Communications on

Pure and Applied Mathematics, vol. 13, pp. 217–237, 1960.

18. A. Jameson, “The construction of discretely conservative finite volume schemes

that also globally conserve energy or entropy,” Journal of Scientific Computing,

vol. 34, pp. 152–187, 2008.

19. B. Gustafsson and H.-O. Kreiss, “Boundary conditions for time dependent problems with an artificial boundary,” Journal of Computational Physics, vol. 30,

pp. 333–351, 1979.

20. M. Sv ̈ard and J. Nordstr ̈om, “On the order of accuracy for difference approximations of initial-boundary value problems,” Journal of Computational Physics,

vol. 218, pp. 333–352, 2006.

21. M. H. Carpenter, J. Nordstr ̈om, and D. Gottlieb, “A stable and conservative

interface treatment of arbitrary spatial accuracy,” Journal of Computational

Physics, vol. 148, pp. 341–365, 1999.

22. N. K. Yamaleev and M. H. Carpenter, “A systematic methodology for constructing high-order energy stable WENO schemes,” Journal of Computational

Physics, vol. 228, pp. 4248–4272, 2009.

23. M. H. Carpenter, J. Nordstr ̈om, and D. Gottlieb, “Revisiting stable and conservative interface treatments of arbitrary spatial accuracy,” Journal of Computational Physics, vol. 45(1-3), pp. 118–150, 2010.

24. T. C. Fisher, M. H. Carpenter, N. K. Yamaleev, and S. H. Frankel, “Boundary

closures for fourth-order energy stable weighted essentially non-oscillatory finite

difference schemes,” Journal of Computational Physics, vol. 230, pp. 3727–3752,

2011.

25. Wolfram Mathematica 7. http://www.wolfram.com/products/mathematica/index.html,

2011. Accessed 06/08/2011.

26. J. Nordstr ̈om, F. Ham, M. Shoeybi, E. van der Weide, M. Sv ̈ard, K. Mattsson, G. Iaccarino, and J. Gong, “A hybrid method for unsteady fluid flow,”

Computers and Fluids, vol. 38, pp. 875–882, 2009.

27. K. Mattsson and M. H. Carpenter, “Stable and accurate interpolation operators for high-order multiblock finite-difference methods,” SIAM J. on Scientific

Computing, vol. 32, p. 2298, 2010.

28. P. Olsson and J. Oliger, “Energy and maximum norm estimates for nonlinear

conservation laws,” tech. rep., Research Institute of Advanced Computer Science, 1994.

29. M. H. Carpenter, D. Gottlieb, and S. Abarbanel, “Time-stable boundary conditions for finite-difference schemes solving hyperbolic systems: Methodology and

24
application to high-order compact schemes,” Journal of Computational Physics,

vol. 111, pp. 220–236, 1994.

30. K. Mattsson and J. Nordstr ̈om, “Summation by parts operators for finite difference approximations of second derivatives,” Journal of Computational Physics,

vol. 199, pp. 503–540, 2004.

31. M. H. Carpenter and C. A. Kennedy, “Fourth-order 2N-storage Runge-Kutta

schemes,” Tech. Rep. TM 109112, NASA, 1994.

32. P. D. Lax, Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves. Philadelphia: SIAM, 1973.

25
Appendix A

Discrete Conservation

The summation-by-parts operators are discretely conservative, as demonstrated

in the following lemma.

Lemma A.1. All differentiation matrices that satisfy the SBP convention given in

eq. (2.1) are discretely conservative in the P-norm.

Proof. Consistency of D implies that all rows in D the matrix sum to zero. Because

P is nonsingular, the following matrix relations hold:

Dsbp1 = P

−1Q1 = Q1 = 0, 1 = (1, · · · , 1)T

, 0 = (0, · · · , 0)T

(A1)

which demonstrates that all rows of Q sum to zero. Forming the dotproduct of

Q + QT with the unit vector and combining the result with that obtained in eq.

(A1) yields

(Q + Q

T

)1 = Q

T 1 = B1 = (−1, 0, . . . , 0, 1)T

(A2)

Taking the transpose of eq. (A2) and combining it with the definition of an SBP

operator yields

1

T Q = 1TPDsbp = (−1, 0, . . . , 0, 1) (A3)

Thus, eq. (A3) demonstrates that the columns of Dsbp are conservative in the P

norm (i.e., all interior columns of the matrix D sum to zero, and the left and right

boundary columns sum to −1 and 1, respectively, after premultiplication by P).

26
Appendix B

Split Continuous Equations

To utilize a split operator, the flux f in eq. (2.3) must be the product of at least

two functions: f(u) = v(u)w(u). Equation (2.3) is split by using the product rule

ut + αf(u)x + (1 − α) [v(u)w(u)x + w(u)v(u)x] = 0. (B1)

This splitting does not change the weak solution to the governing equation, nor does

it alter the Rankine Hugoniot relation. To show this, we start with the conservation

law in one dimension:

d

dt ZxR

xL

udx + f(u)|

xR

xL

= 0, x ∈ [xL, xR], t ∈ [0,∞) (B2)

with a discontinuity located at xd. At this point, note that neither eq. (2.3) nor

eq. (B1) is valid across the discontinuity. Both differential forms are valid on either

side of the discontinuity, with the mild restriction that w(u)x and v(u)x exist in all

regions where f(u)x exists. Equation (B2) can be rewritten [32] as

d

dt

x

−

Z d

xL

u dx +

d

dt ZxR

x

+

d

u dx + f(u)|

xR

xL

= 0

where u is smooth over the interval of integration. If we apply the Leibniz integral

rule, then

x

−

Z d

xL

ut dx + Su(x

−

d

) + ZxR

x

+

d

ut dx − Su(x

+

d

) + f(u)|

xR

xL

= 0

where S =

dxd

dt is the propagation speed of the discontinuity. Either eq. (2.3) or eq.

(B1) can be used to substitute for ut, because the two forms are equal for smooth

solutions. The conservation law becomes

x

−

Z d

xL

−f(u)xdx + Su(x

−

d

) + ZxR

x

+

d

−f(u)xdx − Su(x

+

d

) + f(u)|

xR

xL

= 0

Upon simplification, we get the Rankine Hugoniot relation



f(x

+

d

) − f(x

−

d

)



− S



u(x

+

d

) − u(x

−

d

)



= 0 (B3)

Thus, both eq. (2.3) and eq. (B1) satisfy the sufficient condition of the Lax-Wendroff

theorem which specifies that the continuous weak solution of the discretized equations must satisfy the Rankine Hugoniot relations.

27
Appendix C

High-Order Skew-Symmetric Fluxes

C.1 The Fourth-Order Operator: (3-6-3)

The full (3-6-3) operator is available upon request. The use of electronic format necessitates the brevity of this presentation and easier implementation. For reference,

the interior flux operator is

fi =

1

360





ui−2

ui−1

ui

ui+1

ui+2

ui+3





T 



2 0 0 1 0 0

0 −16 0 −9 1 0

0 0 74 45 −9 1

1 −9 45 74 0 0

0 1 −9 0 −16 0

0 0 1 0 0 2









ui−2

ui−1

ui

ui+1

ui+2

ui+3





(C1)

C.2 The Complete Fifth-Order Operator: (4-8-4)

The full operator can be found in the accompanying text file. The flux for the

interior is

fi =

1

840





ui−3

ui−2

ui−1

ui

ui+1

ui+2

ui+3

ui+4





T 



−1 0 0 0 −

1

2

0 0 0

0

29

3

0 0 16

3 −

1

2

0 0

0 0 −

139

3

0 −28 16

3 −

1

2

0

0 0 0 533

3

112 −28 16

3 −

1

2

−

1

2

16

3 −28 112 5333

0 0 0

0 −

1

2

16

3 −28 0 −

139

3

0 0

0 0 −

1

2

16

3

0 0 29

3

0

0 0 0 −

1

2

0 0 0 −1









ui−3

ui−2

ui−1

ui

ui+1

ui+2

ui+3

ui+4





(C2)

28
Appendix D

Conservative Finite-Difference Viscous Terms

Many conservation laws of interest also have viscous terms. We show that for the

linear viscous terms uxx the summation-by-parts narrow-stencil approximation to

the second derivative can be cast into a conservative form. The SBP finite-difference

operators for the second derivative can be written as

D2 = D

2

or D2 = P

−1

(R + BS), R = RT, v

T Rv ≤ 0 (D1)

where S represents a first derivative operator. In eq. ??, D2 = D2is called the widestencil approximation to the second derivative, which requires that S = D; D2 =

P

−1

(R + BS) is the narrow-stencil approximation to the second derivative and does

not impose the same restriction on S [30]. The narrow-stencil approximation has

superior accuracy and stability properties and should be used where possible.

The conservation form of the narrow stencil is found by constructing an approximation to the viscous flux, ux, at the flux points such that

∆ ̄fV = −∆R ̄u = − (R + BS) u (D2)

This is solved by using a linear system for an appropriate R and BS that satisfies

the SBP condition. The viscous flux

 ̄fV = −R ̄u (D3)

yields a fully conservative discrete approximation to −D2u = ∆ ̄fV . The operator

R ̄ for the 2p = 4 scheme is presented herein, while the R ̄ operators for 2p = 6 and

2p = 8 are appended verbatim in E. For the fourth-order case, we begin with the

existing fourth-order operator of Mattss ̈on et al. [30], which can be expressed with

the relation D2u = ∆R ̄u = (R + BS) u, where

Rδx =





−

9

8

59

48 −

1

12 −

1

48 0 0 0 0 0 0 0

59

48 −

59

24

59

48 0 0 0 0 0 0 0 0

−

1

12

59

48 −

55

24

59

48 −

1

12 0 0 0 0 0 0

−

1

48 0

59

48 −

59

24

4

3 −

1

12 0 0 0 0 0

0 0 −

1

12

4

3 −

5

2

4

3 −

1

12 0 0 0 0

0 0 0

.

.

.

.

.

.

.

.

.

.

.

.

.

.

. 0 0 0

0 0 0 0 −

1

12

4

3 −

5

2

4

3 −

1

12 0 0

0 0 0 0 0 −

1

12

4

3 −

59

24

59

48 0 −

1

48

0 0 0 0 0 0 −

1

12

59

48 −

55

24

59

48 −

1

12

0 0 0 0 0 0 0 0 59

48 −

59

24

59

48

0 0 0 0 0 0 0 −

1

48 −

1

12

59

48 −

9

8





(D4)

29
and

BSδx =





11

6 −3

3

2 −

1

3

0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 −

1

3

3

2 −3

11

6





(D5)

Solving for R ̄ in the relation ∆R ̄u = (R + BS) u yields

R ̄δx =





−

11

6

3 −

3

2

1

3

0 0 0 0 0

−

9

8

59

48 −

1

12 −

1

48 0 0 0 0 0

5

48 −

59

48

55

48 −

1

48 0 0 0 0 0

1

48 0 −

55

48

29

24 −

1

12 0 0 0 0

0

1

12 −

5

4

5

4 −

1

12 0 0 0 0

0 0

.

.

.

.

.

.

.

.

.

.

.

. 0 0 0

0 0 0 1

12 −

5

4

5

4 −

1

12 0 0

0 0 0 0 1

12 −

29

24

55

48 0 −

1

48

0 0 0 0 0 1

48 −

55

48

59

48 −

5

48

0 0 0 0 0 1

48

1

12 −

59

48

9

8

0 0 0 0 0 −

1

3

3

2 −3

11

6





(D6)

Note that this operational form is only design order accurate for a constant diffusivity .

A similar manipulation can be performed by using the sixth- and eighth-order

diffusion operators. These operators are appended verbatim in E.

30
Appendix E

Stencil Coefficients

The following section includes the coefficients for implementation of the split

form conservative finite difference simulations. The given coefficients can be inserted

directly into source code using “cut-and-paste”.

E.1 Definitions

• First derivative of vector u, D1u = P

−1Qu, [see eq. (2.9)] .

• Second derivative of vector u, D2u = P

−1

(R + BS)u, [see eq. (3.1)].

• Split first derivative of w(u)v(u),

αD1W v + (1 − α)(W D1v + V D1w) = P

−1∆  ̄f

, [see eq. (2.24)].

• Conservative Second Derivative, P

−1∆f = −P−1∆R ̄u, [see eq. (3.2)].

For each scheme [ (2-4-2), (3-6-3), and (4-8-4) ], the following are given:

• The flux points for the local stencil, written in terms of a general flux f = v∗w.

• The diagonal of the P-norm P scaled by the grid spacing

• The diagonal of the inverse of the P-norm P scaled by the inverse of grid

spacing

• The nonzero entries of the Q matrix for the SBP first derivative operator (for

reference)

• The nonzero entries of the matrix R ̄, used to calculate the conservative flux

terms

• The nonzero entries of the matrix R, which approximates the second derivative

• The nonzero entries of the matrix B ∗ S, which yields the boundary differentiation operator

Note that only the first s+1 entries are shown. The (s+1)th entry is the interior

operator, and just shifts indices for subsequent entries in the matrix. The respective

symmetry or skew symmetry properties of each matrix is as follows:

• The flux points are per-symmetric. Replace vi and wj with uN−(i−1) and

wN−(j−1), respectively, at right boundary i=1,2,...,N. Note that for Burgers

equation wi = ui and vj = uj/2.

31
• The P-norm P is per-symmetric.

Pi,i = PN−(i−1),N−(i−1), i = 1, 2, ..., N

• The inverse P-norm P

−1

is per-symmetric.

P

−1

i,i = P

−1

N−(i−1),N−(i−1), i = 1, 2, ..., N

• The boundary block portion of Q is per-skew-symmetric.

QN−(i−1),N−(j−1) = −Qi,j , i = 1, 2, ..., s, j = 1, 2, ..., s

• The boundary block portion of R ̄ is per-skew-symmetric.

R ̄N−(i),N−(j−1) = −R ̄

i,j , i = 0, 2, ..., s − 1, j = 1, 2, ..., s

• The boundary block portion of R is per-symmetric.

RN−(i−1),N−(j−1) = Ri,j , i = 1, 2, ..., s, j = 1, 2, ..., s

• The boundary contribution B ∗ S is per-symmetric.

BSN,N−(j−1) = B ∗ S1,j , j = 1, 2, ..., N

E.2 Coefficients

E.2.1 (2-4-2) Scheme Skew-Symmetric

Flux point j:

vj−1∗wj−1 = −0.08333333333333333333333333333∗α+0.0000000000000000000000000000

vj−1∗wj+0 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−1∗wj+1 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

vj−1∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj+0 = 0.5833333333333333333333333333∗α+0.0000000000000000000000000000

vj+0∗wj+1 = −0.6666666666666666666666666667∗α+0.6666666666666666666666666667

vj+0∗wj+2 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

vj+1∗wj−1 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

vj+1∗wj+0 = −0.6666666666666666666666666667∗α+0.6666666666666666666666666667

vj+1∗wj+1 = 0.5833333333333333333333333333∗α+0.0000000000000000000000000000

vj+1∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj+0 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

vj+2∗wj+1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj+2 = −0.08333333333333333333333333333∗α+0.0000000000000000000000000000

32
Flux point j=3

v1∗w1 = −0.03125000000000000000000000000∗α+0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w4 = 0.03125000000000000000000000000∗α−0.03125000000000000000000000000

v1 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w3 = 0.5312500000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w4 = −0.6145833333333333333333333333∗α+0.6145833333333333333333333333

v3∗w5 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

v4∗w1 = 0.03125000000000000000000000000∗α−0.03125000000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w3 = −0.6145833333333333333333333333∗α+0.6145833333333333333333333333

v4 ∗ w4 = 0.5833333333333333333333333333 ∗α + 0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w3 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w5 = −0.08333333333333333333333333333∗α+0.0000000000000000000000000000

Flux point j=2

v1∗w1 = −0.1145833333333333333333333333∗α+0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w3 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

v1∗w4 = 0.03125000000000000000000000000∗α−0.03125000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w2 = 0.6145833333333333333333333333 ∗α + 0.0000000000000000000000000000

v2∗w3 = −0.6145833333333333333333333333∗α+0.6145833333333333333333333333

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w1 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

v3∗w2 = −0.6145833333333333333333333333∗α+0.6145833333333333333333333333

v3 ∗ w3 = 0.5312500000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w1 = 0.03125000000000000000000000000∗α−0.03125000000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w4 = −0.03125000000000000000000000000∗α+0.0000000000000000000000000000

33
Flux point j=1

v1 ∗ w1 = 0.5000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w2 = −0.6145833333333333333333333333∗α+0.6145833333333333333333333333

v1∗w3 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

v1∗w4 = 0.03125000000000000000000000000∗α−0.03125000000000000000000000000

v2∗w1 = −0.6145833333333333333333333333∗α+0.6145833333333333333333333333

v2 ∗ w2 = 0.6145833333333333333333333333 ∗α + 0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w1 = 0.08333333333333333333333333333∗α−0.08333333333333333333333333333

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w3 = −0.08333333333333333333333333333∗α+0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w1 = 0.03125000000000000000000000000∗α−0.03125000000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w4 = −0.03125000000000000000000000000∗α+0.0000000000000000000000000000

Flux point j=0

v1 ∗ w1 = 0.0000000000000000000000000000 ∗ α + 1.000000000000000000000000000

The P-norm P is :

p1,1/dx = 0.35416666666666666666666666666666666667

p2,2/dx = 1.2291666666666666666666666666666666667

p3,3/dx = 0.89583333333333333333333333333333333333

p4,4/dx = 1.0208333333333333333333333333333333333

p5,5/dx = 1.0000000000000000000000000000000000000

The inverse P-norm P

−1

is:

P

−1

1,1

∗ dx = 2.8235294117647058823529411764705882353

P

−1

2,2

∗ dx = 0.81355932203389830508474576271186440678

P

−1

3,3

∗ dx = 1.1162790697674418604651162790697674419

P

−1

4,4

∗ dx = 0.97959183673469387755102040816326530612

P

−1

5,5

∗ dx = 1.0000000000000000000000000000000000000

The Qi,j ∗ dx nonzero terms are:

Q1,1 = −0.50000000000000000000000000000000000000

Q1,2 = 0.61458333333333333333333333333333333333

Q1,3 = −0.083333333333333333333333333333333333333

Q1,4 = −0.031250000000000000000000000000000000000

Q2,1 = −0.61458333333333333333333333333333333333

Q2,3 = 0.61458333333333333333333333333333333333

Q3,1 = 0.083333333333333333333333333333333333333

Q3,2 = −0.61458333333333333333333333333333333333

34
Q3,4 = 0.61458333333333333333333333333333333333

Q3,5 = −0.083333333333333333333333333333333333333

Q4,1 = 0.031250000000000000000000000000000000000

Q4,3 = −0.61458333333333333333333333333333333333

Q4,5 = 0.66666666666666666666666666666666666667

Q4,6 = −0.083333333333333333333333333333333333333

Q5,3 = 0.083333333333333333333333333333333333333

Q5,4 = −0.66666666666666666666666666666666666667

Q5,6 = 0.66666666666666666666666666666666666667

Q5,7 = −0.083333333333333333333333333333333333333

The (2-4-2) viscous terms.

The R ̄

i,j ∗ dx nonzero terms:

R ̄

0,1 = −1.8333333333333333333333333333333333333

R ̄

0,2 = 3.0000000000000000000000000000000000000

R ̄

0,3 = −1.5000000000000000000000000000000000000

R ̄

0,4 = 0.33333333333333333333333333333333333333

R ̄

1,1 = −1.1250000000000000000000000000000000000

R ̄

1,2 = 1.2291666666666666666666666666666666667

R ̄

1,3 = −0.083333333333333333333333333333333333333

R ̄

1,4 = −0.020833333333333333333333333333333333333

R ̄

2,1 = 0.10416666666666666666666666666666666667

R ̄

2,2 = −1.2291666666666666666666666666666666667

R ̄

2,3 = 1.1458333333333333333333333333333333333

R ̄

2,4 = −0.020833333333333333333333333333333333333

R ̄

3,1 = 0.020833333333333333333333333333333333333

R ̄

3,3 = −1.1458333333333333333333333333333333333

R ̄

3,4 = 1.2083333333333333333333333333333333333

R ̄

3,5 = −0.083333333333333333333333333333333333333

R ̄

4,3 = 0.083333333333333333333333333333333333333

R ̄

4,4 = −1.2500000000000000000000000000000000000

R ̄

4,5 = 1.2500000000000000000000000000000000000

R ̄

4,6 = −0.083333333333333333333333333333333333333

The Ri,j ∗ dx nonzero terms are:

R1,1 = −1.1250000000000000000000000000000000000

R1,2 = 1.2291666666666666666666666666666666667

R1,3 = −0.083333333333333333333333333333333333333

R1,4 = −0.020833333333333333333333333333333333333

R2,1 = 1.2291666666666666666666666666666666667

R2,2 = −2.4583333333333333333333333333333333333

R2,3 = 1.2291666666666666666666666666666666667

R3,1 = −0.083333333333333333333333333333333333333

R3,2 = 1.2291666666666666666666666666666666667

R3,3 = −2.2916666666666666666666666666666666667

35
R3,4 = 1.2291666666666666666666666666666666667

R3,5 = −0.083333333333333333333333333333333333333

R4,1 = −0.020833333333333333333333333333333333333

R4,3 = 1.2291666666666666666666666666666666667

R4,4 = −2.4583333333333333333333333333333333333

R4,5 = 1.3333333333333333333333333333333333333

R4,6 = −0.083333333333333333333333333333333333333

R5,3 = −0.083333333333333333333333333333333333333

R5,4 = 1.3333333333333333333333333333333333333

R5,5 = −2.5000000000000000000000000000000000000

R5,6 = 1.3333333333333333333333333333333333333

R5,7 = −0.083333333333333333333333333333333333333

The BSi,j nonzero terms are:

BS1,1 = 1.8333333333333333333333333333333333333

BS1,2 = −3.0000000000000000000000000000000000000

BS1,3 = 1.5000000000000000000000000000000000000

BS1,4 = −0.33333333333333333333333333333333333333

E.2.2 (3-6-3) Scheme Skew-Symmetric

The flux point coefficients at point j are:

vj−2∗wj−2 = 0.01666666666666666666666666667∗α+0.0000000000000000000000000000

vj−2∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−2∗wj+0 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−2∗wj+1 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

vj−2∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−2∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−1∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−1∗wj−1 = −0.1333333333333333333333333333∗α+0.0000000000000000000000000000

vj−1∗wj+0 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−1∗wj+1 = 0.1500000000000000000000000000∗α−0.1500000000000000000000000000

vj−1∗wj+2 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

vj−1∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj+0 = 0.6166666666666666666666666667∗α+0.0000000000000000000000000000

vj+0∗wj+1 = −0.7500000000000000000000000000∗α+0.7500000000000000000000000000

vj+0∗wj+2 = 0.1500000000000000000000000000∗α−0.1500000000000000000000000000

vj+0∗wj+3 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

vj+1∗wj−2 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

vj+1∗wj−1 = 0.1500000000000000000000000000∗α−0.1500000000000000000000000000

vj+1∗wj+0 = −0.7500000000000000000000000000∗α+0.7500000000000000000000000000

vj+1∗wj+1 = 0.6166666666666666666666666667∗α+0.0000000000000000000000000000

36
vj+1∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+1∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj−1 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

vj+2∗wj+0 = 0.1500000000000000000000000000∗α−0.1500000000000000000000000000

vj+2∗wj+1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj+2 = −0.1333333333333333333333333333∗α+0.0000000000000000000000000000

vj+2∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj+0 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

vj+3∗wj+1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj+3 = 0.01666666666666666666666666667∗α+0.0000000000000000000000000000

Flux point j=6

v3∗w3 = −0.01396219135802469135802469136∗α+0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w7 = 0.01396219135802469135802469136∗α−0.01396219135802469135802469136

v3 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗w4 = 0.06997685185185185185185185185∗α+0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w7 = −0.06997685185185185185185185185∗α+0.06997685185185185185185185185

v4 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w5 = −0.2056828703703703703703703704∗α+0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = 0.2223495370370370370370370370 ∗α −0.2223495370370370370370370370

v5∗w8 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

v5 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w6 = 0.6496682098765432098765432099 ∗α + 0.0000000000000000000000000000

v6∗w7 = −0.7830015432098765432098765432∗α+0.7830015432098765432098765432

v6 ∗ w8 = 0.1500000000000000000000000000 ∗α −0.1500000000000000000000000000

v6∗w9 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

v7∗w3 = 0.01396219135802469135802469136∗α−0.01396219135802469135802469136

37
v7∗w4 = −0.06997685185185185185185185185∗α+0.06997685185185185185185185185

v7 ∗ w5 = 0.2223495370370370370370370370 ∗α −0.2223495370370370370370370370

v7∗w6 = −0.7830015432098765432098765432∗α+0.7830015432098765432098765432

v7 ∗ w7 = 0.6166666666666666666666666667 ∗α + 0.0000000000000000000000000000

v7 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w5 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

v8 ∗ w6 = 0.1500000000000000000000000000 ∗α −0.1500000000000000000000000000

v8 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w8 = −0.1333333333333333333333333333∗α+0.0000000000000000000000000000

v8 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9∗w6 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

v9 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗w9 = 0.01666666666666666666666666667∗α+0.0000000000000000000000000000

Flux point j=5

v2∗w2 = −0.02272762345679012345679012346∗α+0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2∗w6 = 0.02272762345679012345679012346∗α−0.02272762345679012345679012346

v2 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗w3 = 0.08048225308641975308641975309∗α+0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w6 = −0.09444444444444444444444444444∗α+0.09444444444444444444444444444

v3∗w7 = 0.01396219135802469135802469136∗α−0.01396219135802469135802469136

v3 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w4 = −0.1737924382716049382716049383∗α+0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.2437692901234567901234567901 ∗α −0.2437692901234567901234567901

v4∗w7 = −0.06997685185185185185185185185∗α+0.06997685185185185185185185185

v4 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

38
v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w5 = 0.6160378086419753086419753086 ∗α + 0.0000000000000000000000000000

v5∗w6 = −0.8217206790123456790123456790∗α+0.8217206790123456790123456790

v5 ∗ w7 = 0.2223495370370370370370370370 ∗α −0.2223495370370370370370370370

v5∗w8 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

v6∗w2 = 0.02272762345679012345679012346∗α−0.02272762345679012345679012346

v6∗w3 = −0.09444444444444444444444444444∗α+0.09444444444444444444444444444

v6 ∗ w4 = 0.2437692901234567901234567901 ∗α −0.2437692901234567901234567901

v6∗w5 = −0.8217206790123456790123456790∗α+0.8217206790123456790123456790

v6 ∗ w6 = 0.6496682098765432098765432099 ∗α + 0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7∗w3 = 0.01396219135802469135802469136∗α−0.01396219135802469135802469136

v7∗w4 = −0.06997685185185185185185185185∗α+0.06997685185185185185185185185

v7 ∗ w5 = 0.2223495370370370370370370370 ∗α −0.2223495370370370370370370370

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7∗w7 = −0.1663348765432098765432098765∗α+0.0000000000000000000000000000

v7 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w5 = −0.01666666666666666666666666667∗α+0.01666666666666666666666666667

v8 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗w8 = 0.01666666666666666666666666667∗α+0.0000000000000000000000000000

Flux point j=4

v1 ∗w1 = 0.06527777777777777777777777778∗α+0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w5 = −0.06527777777777777777777777778∗α+0.06527777777777777777777777778

v1 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2∗w2 = −0.1170100308641975308641975309∗α+0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2∗w5 = 0.09428240740740740740740740741∗α−0.09428240740740740740740740741

v2∗w6 = 0.02272762345679012345679012346∗α−0.02272762345679012345679012346

v2 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w3 = −0.1723919753086419753086419753∗α+0.0000000000000000000000000000

39
v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.2528742283950617283950617284 ∗α −0.2528742283950617283950617284

v3∗w6 = −0.09444444444444444444444444444∗α+0.09444444444444444444444444444

v3∗w7 = 0.01396219135802469135802469136∗α−0.01396219135802469135802469136

v4 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w4 = 0.7241242283950617283950617284 ∗α + 0.0000000000000000000000000000

v4∗w5 = −0.8979166666666666666666666667∗α+0.8979166666666666666666666667

v4 ∗ w6 = 0.2437692901234567901234567901 ∗α −0.2437692901234567901234567901

v4∗w7 = −0.06997685185185185185185185185∗α+0.06997685185185185185185185185

v5∗w1 = −0.06527777777777777777777777778∗α+0.06527777777777777777777777778

v5∗w2 = 0.09428240740740740740740740741∗α−0.09428240740740740740740740741

v5 ∗ w3 = 0.2528742283950617283950617284 ∗α −0.2528742283950617283950617284

v5∗w4 = −0.8979166666666666666666666667∗α+0.8979166666666666666666666667

v5 ∗ w5 = 0.6160378086419753086419753086 ∗α + 0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w2 = 0.02272762345679012345679012346∗α−0.02272762345679012345679012346

v6∗w3 = −0.09444444444444444444444444444∗α+0.09444444444444444444444444444

v6 ∗ w4 = 0.2437692901234567901234567901 ∗α −0.2437692901234567901234567901

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w6 = −0.1720524691358024691358024691∗α+0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7∗w3 = 0.01396219135802469135802469136∗α−0.01396219135802469135802469136

v7∗w4 = −0.06997685185185185185185185185∗α+0.06997685185185185185185185185

v7 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗w7 = 0.05601466049382716049382716049∗α+0.0000000000000000000000000000

Flux point j=3

v1∗w1 = −0.09305555555555555555555555556∗α+0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w4 = 0.1583333333333333333333333333 ∗α −0.1583333333333333333333333333

v1∗w5 = −0.06527777777777777777777777778∗α+0.06527777777777777777777777778

v1 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w2 = 0.2480478395061728395061728395 ∗α + 0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2∗w4 = −0.3650578703703703703703703704∗α+0.3650578703703703703703703704

40
v2∗w5 = 0.09428240740740740740740740741∗α−0.09428240740740740740740740741

v2∗w6 = 0.02272762345679012345679012346∗α−0.02272762345679012345679012346

v2 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w3 = 0.3450077160493827160493827161 ∗α + 0.0000000000000000000000000000

v3∗w4 = −0.5173996913580246913580246914∗α+0.5173996913580246913580246914

v3 ∗ w5 = 0.2528742283950617283950617284 ∗α −0.2528742283950617283950617284

v3∗w6 = −0.09444444444444444444444444444∗α+0.09444444444444444444444444444

v3∗w7 = 0.01396219135802469135802469136∗α−0.01396219135802469135802469136

v4 ∗ w1 = 0.1583333333333333333333333333 ∗α −0.1583333333333333333333333333

v4∗w2 = −0.3650578703703703703703703704∗α+0.3650578703703703703703703704

v4∗w3 = −0.5173996913580246913580246914∗α+0.5173996913580246913580246914

v4 ∗ w4 = 0.7241242283950617283950617284 ∗α + 0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w1 = −0.06527777777777777777777777778∗α+0.06527777777777777777777777778

v5∗w2 = 0.09428240740740740740740740741∗α−0.09428240740740740740740740741

v5 ∗ w3 = 0.2528742283950617283950617284 ∗α −0.2528742283950617283950617284

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w5 = −0.2818788580246913580246913580∗α+0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w2 = 0.02272762345679012345679012346∗α−0.02272762345679012345679012346

v6∗w3 = −0.09444444444444444444444444444∗α+0.09444444444444444444444444444

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗w6 = 0.07171682098765432098765432099∗α+0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7∗w3 = 0.01396219135802469135802469136∗α−0.01396219135802469135802469136

v7 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7∗w7 = −0.01396219135802469135802469136∗α+0.0000000000000000000000000000

Flux point j=2

v1∗w1 = −0.1638888888888888888888888889∗α+0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w3 = 0.07083333333333333333333333333∗α−0.07083333333333333333333333333

v1 ∗ w4 = 0.1583333333333333333333333333 ∗α −0.1583333333333333333333333333

v1∗w5 = −0.06527777777777777777777777778∗α+0.06527777777777777777777777778

41
v1 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w2 = 0.6638888888888888888888888889 ∗α + 0.0000000000000000000000000000

v2∗w3 = −0.4158410493827160493827160494∗α+0.4158410493827160493827160494

v2∗w4 = −0.3650578703703703703703703704∗α+0.3650578703703703703703703704

v2∗w5 = 0.09428240740740740740740740741∗α−0.09428240740740740740740740741

v2∗w6 = 0.02272762345679012345679012346∗α−0.02272762345679012345679012346

v3∗w1 = 0.07083333333333333333333333333∗α−0.07083333333333333333333333333

v3∗w2 = −0.4158410493827160493827160494∗α+0.4158410493827160493827160494

v3 ∗ w3 = 0.3450077160493827160493827161 ∗α + 0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w1 = 0.1583333333333333333333333333 ∗α −0.1583333333333333333333333333

v4∗w2 = −0.3650578703703703703703703704∗α+0.3650578703703703703703703704

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w4 = 0.2067245370370370370370370370 ∗α + 0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w1 = −0.06527777777777777777777777778∗α+0.06527777777777777777777777778

v5∗w2 = 0.09428240740740740740740740741∗α−0.09428240740740740740740740741

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w5 = −0.02900462962962962962962962963∗α+0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w2 = 0.02272762345679012345679012346∗α−0.02272762345679012345679012346

v6 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w6 = −0.02272762345679012345679012346∗α+0.0000000000000000000000000000

Flux point j=1

v1 ∗ w1 = 0.5000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w2 = −0.6638888888888888888888888889∗α+0.6638888888888888888888888889

v1∗w3 = 0.07083333333333333333333333333∗α−0.07083333333333333333333333333

v1 ∗ w4 = 0.1583333333333333333333333333 ∗α −0.1583333333333333333333333333

v1∗w5 = −0.06527777777777777777777777778∗α+0.06527777777777777777777777778

v2∗w1 = −0.6638888888888888888888888889∗α+0.6638888888888888888888888889

v2 ∗ w2 = 0.6638888888888888888888888889 ∗α + 0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w1 = 0.07083333333333333333333333333∗α−0.07083333333333333333333333333

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

42
v3∗w3 = −0.07083333333333333333333333333∗α+0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w1 = 0.1583333333333333333333333333 ∗α −0.1583333333333333333333333333

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w4 = −0.1583333333333333333333333333∗α+0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w1 = −0.06527777777777777777777777778∗α+0.06527777777777777777777777778

v5 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗w5 = 0.06527777777777777777777777778∗α+0.0000000000000000000000000000

Flux point j=0

v1 ∗ w1 = 0.0000000000000000000000000000 ∗ α + 1.000000000000000000000000000

The P-norm P is:

P1,1/dx = 0.30833333333333333333333333333333333333

P2,2/dx = 1.4360879629629629629629629629629629630

P3,3/dx = 0.51331018518518518518518518518518518519

P4,4/dx = 1.3928240740740740740740740740740740741

P5,5/dx = 0.79745370370370370370370370370370370370

P6,6/dx = 1.0596064814814814814814814814814814815

P7,7/dx = 0.99238425925925925925925925925925925926

P8,8/dx = 1.0000000000000000000000000000000000000

The P-norm inverse P

−1

is:

P

−1

1,1

∗ dx = 3.2432432432432432432432432432432432432

P

−1

2,2

∗ dx = 0.69633617563145763148986927577813955737

P

−1

3,3

∗ dx = 1.9481397970687711386696730552423900789

P

−1

4,4

∗ dx = 0.71796576366960279208908093734419145754

P

−1

5,5

∗ dx = 1.2539912917271407837445573294629898403

P

−1

6,6

∗ dx = 0.94374658656471873293282359366466411797

P

−1

7,7

∗ dx = 1.0076741853467378880828532107951762264

P

−1

8,8

∗ dx = 1.0000000000000000000000000000000000000

The Qi,j ∗ dx nonzero terms are:

Q1,1 = −0.50000000000000000000000000000000000000

Q1,2 = 0.66388888888888888888888888888888888889

Q1,3 = −0.070833333333333333333333333333333333333

Q1,4 = −0.15833333333333333333333333333333333333

Q1,5 = 0.065277777777777777777777777777777777778

Q2,1 = −0.66388888888888888888888888888888888889

43
Q2,3 = 0.41584104938271604938271604938271604938

Q2,4 = 0.36505787037037037037037037037037037037

Q2,5 = −0.094282407407407407407407407407407407407

Q2,6 = −0.022727623456790123456790123456790123457

Q3,1 = 0.070833333333333333333333333333333333333

Q3,2 = −0.41584104938271604938271604938271604938

Q3,4 = 0.51739969135802469135802469135802469136

Q3,5 = −0.25287422839506172839506172839506172840

Q3,6 = 0.094444444444444444444444444444444444444

Q3,7 = −0.013962191358024691358024691358024691358

Q4,1 = 0.15833333333333333333333333333333333333

Q4,2 = −0.36505787037037037037037037037037037037

Q4,3 = −0.51739969135802469135802469135802469136

Q4,5 = 0.89791666666666666666666666666666666667

Q4,6 = −0.24376929012345679012345679012345679012

Q4,7 = 0.069976851851851851851851851851851851852

Q5,1 = −0.065277777777777777777777777777777777778

Q5,2 = 0.094282407407407407407407407407407407407

Q5,3 = 0.25287422839506172839506172839506172840

Q5,4 = −0.89791666666666666666666666666666666667

Q5,6 = 0.82172067901234567901234567901234567901

Q5,7 = −0.22234953703703703703703703703703703704

Q5,8 = 0.016666666666666666666666666666666666667

Q6,2 = 0.022727623456790123456790123456790123457

Q6,3 = −0.094444444444444444444444444444444444444

Q6,4 = 0.24376929012345679012345679012345679012

Q6,5 = −0.82172067901234567901234567901234567901

Q6,7 = 0.78300154320987654320987654320987654321

Q6,8 = −0.15000000000000000000000000000000000000

Q6,9 = 0.016666666666666666666666666666666666667

Q7,3 = 0.013962191358024691358024691358024691358

Q7,4 = −0.069976851851851851851851851851851851852

Q7,5 = 0.22234953703703703703703703703703703704

Q7,6 = −0.78300154320987654320987654320987654321

Q7,8 = 0.75000000000000000000000000000000000000

Q7,9 = −0.15000000000000000000000000000000000000

Q7,10 = 0.016666666666666666666666666666666666667

Q8,5 = −0.016666666666666666666666666666666666667

Q8,6 = 0.15000000000000000000000000000000000000

Q8,7 = −0.75000000000000000000000000000000000000

Q8,9 = 0.75000000000000000000000000000000000000

Q8,10 = −0.15000000000000000000000000000000000000

Q8,11 = 0.016666666666666666666666666666666666667

The (3-6-3) viscous terms.

44
The R ̄

i,j ∗ dx nonzero terms:

R ̄

0,1 = −2.2833333333333333333333333333333333333

R ̄

0,2 = 5.0000000000000000000000000000000000000

R ̄

0,3 = −5.0000000000000000000000000000000000000

R ̄

0,4 = 3.3333333333333333333333333333333333333

R ̄

0,5 = −1.2500000000000000000000000000000000000

R ̄

0,6 = 0.20000000000000000000000000000000000000

R ̄

1,1 = −1.1840277777777777777777777777777777778

R ̄

1,2 = 1.3277777777777777777777777777777777778

R ̄

1,3 = −0.070833333333333333333333333333333333333

R ̄

1,4 = −0.10555555555555555555555555555555555556

R ̄

1,5 = 0.032638888888888888888888888888888888889

R ̄

2,1 = 0.14375000000000000000000000000000000000

R ̄

2,2 = −1.1225212191358024691358024691358024691

R ̄

2,3 = 0.76084876543209876543209876543209876543

R ̄

2,4 = 0.25950231481481481481481481481481481481

R ̄

2,5 = −0.030216049382716049382716049382716049383

R ̄

2,6 = −0.011363811728395061728395061728395061728

R ̄

3,1 = 0.072916666666666666666666666666666666667

R ̄

3,2 = −0.29083912037037037037037037037037037037

R ̄

3,3 = −0.83790702160493827160493827160493827160

R ̄

3,4 = 1.2943016975308641975308641975308641975

R ̄

3,5 = −0.28309027777777777777777777777777777778

R ̄

3,6 = 0.051599151234567901234567901234567901235

R ̄

3,7 = −0.0069810956790123456790123456790123456790

R ̄

4,1 = −0.032638888888888888888888888888888888889

R ̄

4,2 = 0.074218750000000000000000000000000000000

R ̄

4,3 = 0.19689236111111111111111111111111111111

R ̄

4,4 = −1.5987152777777777777777777777777777778

R ̄

4,5 = 1.5127430555555555555555555555555555556

R ̄

4,6 = −0.19217013888888888888888888888888888889

R ̄

4,7 = 0.039670138888888888888888888888888888889

R ̄

5,2 = 0.011363811728395061728395061728395061728

R ̄

5,3 = −0.055981867283950617283950617283950617284

R ̄

5,4 = 0.19711805555555555555555555555555555556

R ̄

5,5 = −1.4322029320987654320987654320987654321

R ̄

5,6 = 1.4512712191358024691358024691358024691

R ̄

5,7 = −0.18267939814814814814814814814814814815

R ̄

5,8 = 0.011111111111111111111111111111111111111

R ̄

6,3 = 0.0069810956790123456790123456790123456790

R ̄

6,4 = −0.046651234567901234567901234567901234568

R ̄

6,5 = 0.21123842592592592592592592592592592593

R ̄

6,6 = −1.4271141975308641975308641975308641975

R ̄

6,7 = 1.3833236882716049382716049382716049383

R ̄

6,8 = −0.13888888888888888888888888888888888889

45
R ̄

6,9 = 0.011111111111111111111111111111111111111

R ̄

7,5 = −0.011111111111111111111111111111111111111

R ̄

7,6 = 0.13888888888888888888888888888888888889

R ̄

7,7 = −1.3611111111111111111111111111111111111

R ̄

7,8 = 1.3611111111111111111111111111111111111

R ̄

7,9 = −0.13888888888888888888888888888888888889

R ̄

7,10 = 0.011111111111111111111111111111111111111

The Ri,j ∗ dx nonzero terms:

R1,1 = −1.1840277777777777777777777777777777778

R1,2 = 1.3277777777777777777777777777777777778

R1,3 = −0.070833333333333333333333333333333333333

R1,4 = −0.10555555555555555555555555555555555556

R1,5 = 0.032638888888888888888888888888888888889

R2,1 = 1.3277777777777777777777777777777777778

R2,2 = −2.4502989969135802469135802469135802469

R2,3 = 0.83168209876543209876543209876543209877

R2,4 = 0.36505787037037037037037037037037037037

R2,5 = −0.062854938271604938271604938271604938272

R2,6 = −0.011363811728395061728395061728395061728

R3,1 = −0.070833333333333333333333333333333333333

R3,2 = 0.83168209876543209876543209876543209877

R3,3 = −1.5987557870370370370370370370370370370

R3,4 = 1.0347993827160493827160493827160493827

R3,5 = −0.25287422839506172839506172839506172840

R3,6 = 0.062962962962962962962962962962962962963

R3,7 = −0.0069810956790123456790123456790123456790

R4,1 = −0.10555555555555555555555555555555555556

R4,2 = 0.36505787037037037037037037037037037037

R4,3 = 1.0347993827160493827160493827160493827

R4,4 = −2.8930169753086419753086419753086419753

R4,5 = 1.7958333333333333333333333333333333333

R4,6 = −0.24376929012345679012345679012345679012

R4,7 = 0.046651234567901234567901234567901234568

R5,1 = 0.032638888888888888888888888888888888889

R5,2 = −0.062854938271604938271604938271604938272

R5,3 = −0.25287422839506172839506172839506172840

R5,4 = 1.7958333333333333333333333333333333333

R5,5 = −2.9449459876543209876543209876543209877

R5,6 = 1.6434413580246913580246913580246913580

R5,7 = −0.22234953703703703703703703703703703704

R5,8 = 0.011111111111111111111111111111111111111

R6,2 = −0.011363811728395061728395061728395061728

R6,3 = 0.062962962962962962962962962962962962963

R6,4 = −0.24376929012345679012345679012345679012

46
R6,5 = 1.6434413580246913580246913580246913580

R6,6 = −2.8783854166666666666666666666666666667

R6,7 = 1.5660030864197530864197530864197530864

R6,8 = −0.15000000000000000000000000000000000000

R6,9 = 0.011111111111111111111111111111111111111

R7,3 = −0.0069810956790123456790123456790123456790

R7,4 = 0.046651234567901234567901234567901234568

R7,5 = −0.22234953703703703703703703703703703704

R7,6 = 1.5660030864197530864197530864197530864

R7,7 = −2.7444347993827160493827160493827160494

R7,8 = 1.5000000000000000000000000000000000000

R7,9 = −0.15000000000000000000000000000000000000

R7,10 = 0.011111111111111111111111111111111111111

R8,5 = 0.011111111111111111111111111111111111111

R8,6 = −0.15000000000000000000000000000000000000

R8,7 = 1.5000000000000000000000000000000000000

R8,8 = −2.7222222222222222222222222222222222222

R8,9 = 1.5000000000000000000000000000000000000

R8,10 = −0.15000000000000000000000000000000000000

R8,11 = 0.011111111111111111111111111111111111111

The BSi,j nonzero terms

BS1,1 = 2.2833333333333333333333333333333333333

BS1,2 = −5.0000000000000000000000000000000000000

BS1,3 = 5.0000000000000000000000000000000000000

BS1,4 = −3.3333333333333333333333333333333333333

BS1,5 = 1.2500000000000000000000000000000000000

BS1,6 = −0.20000000000000000000000000000000000000

E.2.3 (4-8-4) Scheme Skew-Symmetric

The coefficients for flux point j are :

vj−3∗wj−3 = −0.003571428571428571428571428571∗α+0.0000000000000000000000000000

vj−3∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−3∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−3∗wj+0 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−3∗wj+1 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

vj−3∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−3∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−3∗wj+4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−2∗wj−3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−2∗wj−2 = 0.03452380952380952380952380952∗α+0.0000000000000000000000000000

vj−2∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−2∗wj+0 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

47
vj−2∗wj+1 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

vj−2∗wj+2 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

vj−2∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−2∗wj+4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−1∗wj−3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−1∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−1∗wj−1 = −0.1654761904761904761904761905∗α+0.0000000000000000000000000000

vj−1∗wj+0 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj−1∗wj+1 = 0.2000000000000000000000000000∗α−0.2000000000000000000000000000

vj−1∗wj+2 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

vj−1∗wj+3 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

vj−1∗wj+4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj−3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+0∗wj+0 = 0.6345238095238095238095238095∗α+0.0000000000000000000000000000

vj+0∗wj+1 = −0.8000000000000000000000000000∗α+0.8000000000000000000000000000

vj+0∗wj+2 = 0.2000000000000000000000000000∗α−0.2000000000000000000000000000

vj+0∗wj+3 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

vj+0∗wj+4 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

vj+1∗wj−3 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

vj+1∗wj−2 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

vj+1∗wj−1 = 0.2000000000000000000000000000∗α−0.2000000000000000000000000000

vj+1∗wj+0 = −0.8000000000000000000000000000∗α+0.8000000000000000000000000000

vj+1∗wj+1 = 0.6345238095238095238095238095∗α+0.0000000000000000000000000000

vj+1∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+1∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+1∗wj+4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj−3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj−2 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

vj+2∗wj−1 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

vj+2∗wj+0 = 0.2000000000000000000000000000∗α−0.2000000000000000000000000000

vj+2∗wj+1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj+2 = −0.1654761904761904761904761905∗α+0.0000000000000000000000000000

vj+2∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+2∗wj+4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj−3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj−1 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

vj+3∗wj+0 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

vj+3∗wj+1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+3∗wj+3 = 0.03452380952380952380952380952∗α+0.0000000000000000000000000000

vj+3∗wj+4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+4∗wj−3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

48
vj+4∗wj−2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+4∗wj−1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+4∗wj+0 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

vj+4∗wj+1 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+4∗wj+2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+4∗wj+3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

vj+4∗wj+4 = −0.003571428571428571428571428571∗α+0.0000000000000000000000000000

Flux point j=8

v4∗w4 = −0.009602757962753002435542118082∗α+0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w9 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v4 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v4 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v4 ∗w12 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗w5 = 0.04329003028680608045687410767∗α+0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w9 = −0.04329003028680608045687410767∗α+0.04329003028680608045687410767

v5 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v5 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v5 ∗w12 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w6 = −0.05535800500755857898715041572∗α+0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w9 = 0.05178657643613000755857898715∗α−0.05178657643613000755857898715

v6∗w10 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v6 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v6 ∗w12 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v7 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7∗w7 = −0.08327658231502477534223565970∗α+0.0000000000000000000000000000

v7 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w9 = 0.1178003918388342991517594692 ∗α −0.1178003918388342991517594692

v7∗w10 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

v7∗w11 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v7 ∗w12 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

49
v8 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w8 = 0.6049473149985302763080540858 ∗α + 0.0000000000000000000000000000

v8∗w9 = −0.7704235054747207524985302763∗α+0.7704235054747207524985302763

v8 ∗w10 = 0.2000000000000000000000000000∗α−0.2000000000000000000000000000

v8∗w11 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

v8∗w12 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v9∗w4 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v9∗w5 = −0.04329003028680608045687410767∗α+0.04329003028680608045687410767

v9∗w6 = 0.05178657643613000755857898715∗α−0.05178657643613000755857898715

v9 ∗ w7 = 0.1178003918388342991517594692 ∗α −0.1178003918388342991517594692

v9∗w8 = −0.7704235054747207524985302763∗α+0.7704235054747207524985302763

v9 ∗ w9 = 0.6345238095238095238095238095 ∗α + 0.0000000000000000000000000000

v9 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v9 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v9 ∗w12 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w5 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10∗w6 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v10∗w7 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

v10 ∗w8 = 0.2000000000000000000000000000∗α−0.2000000000000000000000000000

v10 ∗w9 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10∗w10 = −0.1654761904761904761904761905∗α+0.0000000000000000000000000000

v10∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10∗w12 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11 ∗w4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11 ∗w5 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11 ∗w6 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11∗w7 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v11∗w8 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

v11 ∗w9 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11∗w11 = 0.03452380952380952380952380952∗α+0.0000000000000000000000000000

v11∗w12 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v12 ∗w4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v12 ∗w5 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v12 ∗w6 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v12 ∗w7 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v12∗w8 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v12 ∗w9 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v12∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v12∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v12∗w12 = −0.003571428571428571428571428571∗α+0.0000000000000000000000000000

50
Flux point j=7

v3∗w3 = −0.09050595238095238095238095238∗α+0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w8 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v3 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v3 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w4 = 0.4166642632747543461829176115 ∗α + 0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w8 = −0.4262670212375073486184597296∗α+0.4262670212375073486184597296

v4∗w9 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v4 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v4 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w5 = −0.7125127557581674645166708659∗α+0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w8 = 0.7558027860449735449735449735 ∗α −0.7558027860449735449735449735

v5∗w9 = −0.04329003028680608045687410767∗α+0.04329003028680608045687410767

v5 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v5 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v6 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w6 = 0.4365587627026119089611153103 ∗α + 0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w8 = −0.4919167677101704879482657260∗α+0.4919167677101704879482657260

v6∗w9 = 0.05178657643613000755857898715∗α−0.05178657643613000755857898715

v6∗w10 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v6 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v7 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w7 = 0.4497956821617535903250188964 ∗α + 0.0000000000000000000000000000

v7∗w8 = −0.5330722644767783656672545561∗α+0.5330722644767783656672545561

v7 ∗ w9 = 0.1178003918388342991517594692 ∗α −0.1178003918388342991517594692

51
v7∗w10 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

v7∗w11 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v8∗w3 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v8∗w4 = −0.4262670212375073486184597296∗α+0.4262670212375073486184597296

v8 ∗ w5 = 0.7558027860449735449735449735 ∗α −0.7558027860449735449735449735

v8∗w6 = −0.4919167677101704879482657260∗α+0.4919167677101704879482657260

v8∗w7 = −0.5330722644767783656672545561∗α+0.5330722644767783656672545561

v8 ∗ w8 = 0.6049473149985302763080540858 ∗α + 0.0000000000000000000000000000

v8 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v8 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v9 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9∗w4 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v9∗w5 = −0.04329003028680608045687410767∗α+0.04329003028680608045687410767

v9∗w6 = 0.05178657643613000755857898715∗α−0.05178657643613000755857898715

v9 ∗ w7 = 0.1178003918388342991517594692 ∗α −0.1178003918388342991517594692

v9 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9∗w9 = −0.1358996959509112286890064668∗α+0.0000000000000000000000000000

v9 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v9 ∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w5 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10∗w6 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v10∗w7 = −0.03809523809523809523809523810∗α+0.03809523809523809523809523810

v10 ∗w8 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w9 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10∗w10 = 0.03452380952380952380952380952∗α+0.0000000000000000000000000000

v10∗w11 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11 ∗w3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11 ∗w4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11 ∗w5 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11 ∗w6 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11∗w7 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v11 ∗w8 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11 ∗w9 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v11∗w11 = −0.003571428571428571428571428571∗α+0.0000000000000000000000000000

Flux point j=6

v2∗w2 = −0.2107209027515327118501721676∗α+0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

52
v2 ∗ w7 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v2 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w3 = 1.090508276381540270429159318 ∗ α + 0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w7 = −1.181014228762492651381540270 ∗ α + 1.181014228762492651381540270

v3∗w8 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v3 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗w4 = −2.220792581963340891912320484∗α+ 0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w7 = 2.637456845238095238095238095 ∗ α − 2.637456845238095238095238095

v4∗w8 = −0.4262670212375073486184597296∗α+0.4262670212375073486184597296

v4∗w9 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v4 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v5 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w5 = 2.067331343301209372637944067 ∗ α + 0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = −2.779844099059376837154614932 ∗ α + 2.779844099059376837154614932

v5 ∗ w8 = 0.7558027860449735449735449735 ∗α −0.7558027860449735449735449735

v5∗w9 = −0.04329003028680608045687410767∗α+0.04329003028680608045687410767

v5 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v6 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w6 = −0.2263261349678760393046107332∗α+0.0000000000000000000000000000

v6 ∗ w7 = 0.6628848976704879482657260435 ∗α −0.6628848976704879482657260435

v6∗w8 = −0.4919167677101704879482657260∗α+0.4919167677101704879482657260

v6∗w9 = 0.05178657643613000755857898715∗α−0.05178657643613000755857898715

v6∗w10 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v7 ∗ w2 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v7 ∗ w3 = −1.181014228762492651381540270 ∗ α + 1.181014228762492651381540270

v7 ∗ w4 = 2.637456845238095238095238095 ∗ α − 2.637456845238095238095238095

v7 ∗ w5 = −2.779844099059376837154614932 ∗ α + 2.779844099059376837154614932

v7 ∗ w6 = 0.6628848976704879482657260435 ∗α −0.6628848976704879482657260435

53
v7 ∗ w7 = 0.4497956821617535903250188964 ∗α + 0.0000000000000000000000000000

v7 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v8 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w3 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v8∗w4 = −0.4262670212375073486184597296∗α+0.4262670212375073486184597296

v8 ∗ w5 = 0.7558027860449735449735449735 ∗α −0.7558027860449735449735449735

v8∗w6 = −0.4919167677101704879482657260∗α+0.4919167677101704879482657260

v8 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗w8 = 0.07187505052175191064079952969∗α+0.0000000000000000000000000000

v8 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v9 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9∗w4 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v9∗w5 = −0.04329003028680608045687410767∗α+0.04329003028680608045687410767

v9∗w6 = 0.05178657643613000755857898715∗α−0.05178657643613000755857898715

v9 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9∗w9 = −0.01809930411207692953724699756∗α+0.0000000000000000000000000000

v9 ∗w10 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w2 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w3 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w4 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w5 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10∗w6 = 0.003571428571428571428571428571∗α−0.003571428571428571428571428571

v10 ∗w7 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w8 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10 ∗w9 = 0.0000000000000000000000000000∗α+0.0000000000000000000000000000

v10∗w10 = −0.003571428571428571428571428571∗α+0.0000000000000000000000000000

Flux point j=5

v1∗w1 = −0.1239583333333333333333333333∗α+0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w6 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v1 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w2 = 1.048880436009280255312001344 ∗ α + 0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

54
v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w6 = −1.259601338760812967162173511 ∗ α + 1.259601338760812967162173511

v2 ∗ w7 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v2 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗w3 = −3.083680753049676660787771899∗α+ 0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w6 = 4.174189029431216931216931217 ∗ α − 4.174189029431216931216931217

v3 ∗ w7 = −1.181014228762492651381540270 ∗ α + 1.181014228762492651381540270

v3∗w8 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v3 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w4 = 3.940000223739187032837826489 ∗ α + 0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = −6.160792805702527924750146972 ∗ α + 6.160792805702527924750146972

v4 ∗ w7 = 2.637456845238095238095238095 ∗ α − 2.637456845238095238095238095

v4∗w8 = −0.4262670212375073486184597296∗α+0.4262670212375073486184597296

v4∗w9 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v5 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗w5 = −1.281241573365457294028722600∗α+ 0.0000000000000000000000000000

v5 ∗ w6 = 3.348572916666666666666666667 ∗ α − 3.348572916666666666666666667

v5 ∗ w7 = −2.779844099059376837154614932 ∗ α + 2.779844099059376837154614932

v5 ∗ w8 = 0.7558027860449735449735449735 ∗α −0.7558027860449735449735449735

v5∗w9 = −0.04329003028680608045687410767∗α+0.04329003028680608045687410767

v6 ∗ w1 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v6 ∗ w2 = −1.259601338760812967162173511 ∗ α + 1.259601338760812967162173511

v6 ∗ w3 = 4.174189029431216931216931217 ∗ α − 4.174189029431216931216931217

v6 ∗ w4 = −6.160792805702527924750146972 ∗ α + 6.160792805702527924750146972

v6 ∗ w5 = 3.348572916666666666666666667 ∗ α − 3.348572916666666666666666667

v6∗w6 = −0.2263261349678760393046107332∗α+0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w2 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v7 ∗ w3 = −1.181014228762492651381540270 ∗ α + 1.181014228762492651381540270

55
v7 ∗ w4 = 2.637456845238095238095238095 ∗ α − 2.637456845238095238095238095

v7 ∗ w5 = −2.779844099059376837154614932 ∗ α + 2.779844099059376837154614932

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w7 = 1.112680579832241538590744940 ∗ α + 0.0000000000000000000000000000

v7 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w3 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v8∗w4 = −0.4262670212375073486184597296∗α+0.4262670212375073486184597296

v8 ∗ w5 = 0.7558027860449735449735449735 ∗α −0.7558027860449735449735449735

v8 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w8 = −0.4200417171884185773074661964∗α+0.0000000000000000000000000000

v8 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9∗w4 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v9∗w5 = −0.04329003028680608045687410767∗α+0.04329003028680608045687410767

v9 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗w9 = 0.03368727232405307802133198959∗α+0.0000000000000000000000000000

Flux point j=4

v1 ∗ w1 = 0.4209583333333333333333333333 ∗α + 0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w5 = −0.5449166666666666666666666667∗α+0.5449166666666666666666666667

v1 ∗ w6 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v1 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗w2 = −1.964242213745065927605610145∗α+ 0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w5 = 3.013122649754346182917611489 ∗ α − 3.013122649754346182917611489

v2 ∗ w6 = −1.259601338760812967162173511 ∗ α + 1.259601338760812967162173511

v2 ∗ w7 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v2 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

56
v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w3 = 3.207549800668724279835390946 ∗ α + 0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = −6.291230553718400940623162845 ∗ α + 6.291230553718400940623162845

v3 ∗ w6 = 4.174189029431216931216931217 ∗ α − 4.174189029431216931216931217

v3 ∗ w7 = −1.181014228762492651381540270 ∗ α + 1.181014228762492651381540270

v3∗w8 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v3 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗w4 = −1.164265920256991685563114135∗α+ 0.0000000000000000000000000000

v4 ∗ w5 = 5.104266143996178718400940623 ∗ α − 5.104266143996178718400940623

v4 ∗ w6 = −6.160792805702527924750146972 ∗ α + 6.160792805702527924750146972

v4 ∗ w7 = 2.637456845238095238095238095 ∗ α − 2.637456845238095238095238095

v4∗w8 = −0.4262670212375073486184597296∗α+0.4262670212375073486184597296

v4∗w9 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v5∗w1 = −0.5449166666666666666666666667∗α+0.5449166666666666666666666667

v5 ∗ w2 = 3.013122649754346182917611489 ∗ α − 3.013122649754346182917611489

v5 ∗ w3 = −6.291230553718400940623162845 ∗ α + 6.291230553718400940623162845

v5 ∗ w4 = 5.104266143996178718400940623 ∗ α − 5.104266143996178718400940623

v5 ∗w5 = −1.281241573365457294028722600∗α+ 0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w1 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v6 ∗ w2 = −1.259601338760812967162173511 ∗ α + 1.259601338760812967162173511

v6 ∗ w3 = 4.174189029431216931216931217 ∗ α − 4.174189029431216931216931217

v6 ∗ w4 = −6.160792805702527924750146972 ∗ α + 6.160792805702527924750146972

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w6 = 3.122246781698790627362055934 ∗ α + 0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w2 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v7 ∗ w3 = −1.181014228762492651381540270 ∗ α + 1.181014228762492651381540270

v7 ∗ w4 = 2.637456845238095238095238095 ∗ α − 2.637456845238095238095238095

v7 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗w7 = −1.667163519227135298563869992∗α+ 0.0000000000000000000000000000

v7 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

57
v8 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w3 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v8∗w4 = −0.4262670212375073486184597296∗α+0.4262670212375073486184597296

v8 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w8 = 0.3357610688565549676660787772 ∗α + 0.0000000000000000000000000000

v8 ∗ w9 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9∗w4 = 0.009602757962753002435542118082∗α−0.009602757962753002435542118082

v9 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v9∗w9 = −0.009602757962753002435542118082∗α+0.0000000000000000000000000000

Flux point j=3

v1∗w1 = −0.4192916666666666666666666667∗α+0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w4 = 0.8402500000000000000000000000 ∗α −0.8402500000000000000000000000

v1∗w5 = −0.5449166666666666666666666667∗α+0.5449166666666666666666666667

v1 ∗ w6 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v1 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w2 = 1.458512047141387419165196943 ∗ α + 0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = −3.422754260886453346770807088 ∗ α + 3.422754260886453346770807088

v2 ∗ w5 = 3.013122649754346182917611489 ∗ α − 3.013122649754346182917611489

v2 ∗ w6 = −1.259601338760812967162173511 ∗ α + 1.259601338760812967162173511

v2 ∗ w7 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v2 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w3 = −0.5392203804747207524985302763∗α+0.0000000000000000000000000000

v3 ∗ w4 = 3.746770181143445032333921223 ∗ α − 3.746770181143445032333921223

v3 ∗ w5 = −6.291230553718400940623162845 ∗ α + 6.291230553718400940623162845

v3 ∗ w6 = 4.174189029431216931216931217 ∗ α − 4.174189029431216931216931217

v3 ∗ w7 = −1.181014228762492651381540270 ∗ α + 1.181014228762492651381540270

v3∗w8 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v4 ∗ w1 = 0.8402500000000000000000000000 ∗α −0.8402500000000000000000000000

v4 ∗ w2 = −3.422754260886453346770807088 ∗ α + 3.422754260886453346770807088

58
v4 ∗ w3 = 3.746770181143445032333921223 ∗ α − 3.746770181143445032333921223

v4 ∗w4 = −1.164265920256991685563114135∗α+ 0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w1 = −0.5449166666666666666666666667∗α+0.5449166666666666666666666667

v5 ∗ w2 = 3.013122649754346182917611489 ∗ α − 3.013122649754346182917611489

v5 ∗ w3 = −6.291230553718400940623162845 ∗ α + 6.291230553718400940623162845

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w5 = 3.823024570630721424372218023 ∗ α + 0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w1 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v6 ∗ w2 = −1.259601338760812967162173511 ∗ α + 1.259601338760812967162173511

v6 ∗ w3 = 4.174189029431216931216931217 ∗ α − 4.174189029431216931216931217

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗w6 = −3.038546024003737297388091039∗α+ 0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w2 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v7 ∗ w3 = −1.181014228762492651381540270 ∗ α + 1.181014228762492651381540270

v7 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w7 = 0.9702933260109599395313681028 ∗α + 0.0000000000000000000000000000

v7 ∗ w8 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w3 = 0.09050595238095238095238095238∗α−0.09050595238095238095238095238

v8 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v8∗w8 = −0.09050595238095238095238095238∗α+0.0000000000000000000000000000

Flux point j=2

v1∗w1 = −0.07820833333333333333333333333∗α+0.0000000000000000000000000000

v1 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w3 = −0.3410833333333333333333333333∗α+0.3410833333333333333333333333

v1 ∗ w4 = 0.8402500000000000000000000000 ∗α −0.8402500000000000000000000000

v1∗w5 = −0.5449166666666666666666666667∗α+0.5449166666666666666666666667

59
v1 ∗ w6 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v1 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w2 = 0.5782083333333333333333333333 ∗α + 0.0000000000000000000000000000

v2 ∗ w3 = 0.8803037138080540858318636096 ∗α −0.8803037138080540858318636096

v2 ∗ w4 = −3.422754260886453346770807088 ∗ α + 3.422754260886453346770807088

v2 ∗ w5 = 3.013122649754346182917611489 ∗ α − 3.013122649754346182917611489

v2 ∗ w6 = −1.259601338760812967162173511 ∗ α + 1.259601338760812967162173511

v2 ∗ w7 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v3∗w1 = −0.3410833333333333333333333333∗α+0.3410833333333333333333333333

v3 ∗ w2 = 0.8803037138080540858318636096 ∗α −0.8803037138080540858318636096

v3∗w3 = −0.5392203804747207524985302763∗α+0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w1 = 0.8402500000000000000000000000 ∗α −0.8402500000000000000000000000

v4 ∗ w2 = −3.422754260886453346770807088 ∗ α + 3.422754260886453346770807088

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w4 = 2.582504260886453346770807088 ∗ α + 0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w1 = −0.5449166666666666666666666667∗α+0.5449166666666666666666666667

v5 ∗ w2 = 3.013122649754346182917611489 ∗ α − 3.013122649754346182917611489

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗w5 = −2.468205983087679516250944822∗α+ 0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w1 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v6 ∗ w2 = −1.259601338760812967162173511 ∗ α + 1.259601338760812967162173511

v6 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w6 = 1.135643005427479633828840178 ∗ α + 0.0000000000000000000000000000

v6 ∗ w7 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w1 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w2 = 0.2107209027515327118501721676 ∗α −0.2107209027515327118501721676

v7 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v7∗w7 = −0.2107209027515327118501721676∗α+0.0000000000000000000000000000

60
Flux point j=1

v1 ∗ w1 = 0.5000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v1∗w2 = −0.5782083333333333333333333333∗α+0.5782083333333333333333333333

v1∗w3 = −0.3410833333333333333333333333∗α+0.3410833333333333333333333333

v1 ∗ w4 = 0.8402500000000000000000000000 ∗α −0.8402500000000000000000000000

v1∗w5 = −0.5449166666666666666666666667∗α+0.5449166666666666666666666667

v1 ∗ w6 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v2∗w1 = −0.5782083333333333333333333333∗α+0.5782083333333333333333333333

v2 ∗ w2 = 0.5782083333333333333333333333 ∗α + 0.0000000000000000000000000000

v2 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v2 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3∗w1 = −0.3410833333333333333333333333∗α+0.3410833333333333333333333333

v3 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w3 = 0.3410833333333333333333333333 ∗α + 0.0000000000000000000000000000

v3 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v3 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w1 = 0.8402500000000000000000000000 ∗α −0.8402500000000000000000000000

v4 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4∗w4 = −0.8402500000000000000000000000∗α+0.0000000000000000000000000000

v4 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v4 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5∗w1 = −0.5449166666666666666666666667∗α+0.5449166666666666666666666667

v5 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v5 ∗ w5 = 0.5449166666666666666666666667 ∗α + 0.0000000000000000000000000000

v5 ∗ w6 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w1 = 0.1239583333333333333333333333 ∗α −0.1239583333333333333333333333

v6 ∗ w2 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w3 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w4 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6 ∗ w5 = 0.0000000000000000000000000000 ∗α + 0.0000000000000000000000000000

v6∗w6 = −0.1239583333333333333333333333∗α+0.0000000000000000000000000000

Flux point j=0

v1 ∗ w1 = 0.0000000000000000000000000000 ∗ α + 1.000000000000000000000000000

The P-norm P is:

P1,1/dx = 0.29950000000000000000000000000000000000

P2,2/dx = 1.4888460333207357016880826404635928445

P3,3/dx = 0.38651394400352733686067019400352733686

61
P4,4/dx = 1.5399915674603174603174603174603174603

P5,5/dx = 0.73536072530864197530864197530864197531

P6,6/dx = 1.0203624889770723104056437389770723104

P7,7/dx = 1.0523566468253968253968253968253968254

P8,8/dx = 0.97245927028218694885361552028218694885

P9,9/dx = 1.0046093238221214411690602166792642983

P10,10/dx = 1.0000000000000000000000000000000000000

The P-norm inverse P

−1

is:

P

−1

1,1

∗ dx = 3.3388981636060100166944908180300500835

P

−1

2,2

∗ dx = 0.67166112386355419609800361200830541810

P

−1

3,3

∗ dx = 2.5872287805246011319981479868054754454

P

−1

4,4

∗ dx = 0.64935420500331280158653723817674480927

P

−1

5,5

∗ dx = 1.3598768136281481908969048637106484335

P

−1

6,6

∗ dx = 0.98004386754996649194151043751580269422

P

−1

7,7

∗ dx = 0.95024819106399044849734934488965855641

P

−1

8,8

∗ dx = 1.0283207025316559720897982867102527893

P

−1

9,9

∗ dx = 0.99541182456421479599641740148968966041

P

−1

10,10 ∗ dx = 1.0000000000000000000000000000000000000

The Qi,j ∗ dx nonzero terms are:

Q1,1 = −0.50000000000000000000000000000000000000

Q1,2 = 0.57820833333333333333333333333333333333

Q1,3 = 0.34108333333333333333333333333333333333

Q1,4 = −0.84025000000000000000000000000000000000

Q1,5 = 0.54491666666666666666666666666666666667

Q1,6 = −0.12395833333333333333333333333333333333

Q2,1 = −0.57820833333333333333333333333333333333

Q2,3 = −0.88030371380805408583186360964138741917

Q2,4 = 3.4227542608864533467708070882674057277

Q2,5 = −3.0131226497543461829176114890400604686

Q2,6 = 1.2596013387608129671621735113798605862

Q2,7 = −0.21072090275153271185017216763248509280

Q3,1 = −0.34108333333333333333333333333333333333

Q3,2 = 0.88030371380805408583186360964138741917

Q3,4 = −3.7467701811434450323339212228101116990

Q3,5 = 6.2912305537184009406231628453850676073

Q3,6 = −4.1741890294312169312169312169312169312

Q3,7 = 1.1810142287624926513815402704291593181

Q3,8 = −0.090505952380952380952380952380952380952

Q4,1 = 0.84025000000000000000000000000000000000

Q4,2 = −3.4227542608864533467708070882674057277

Q4,3 = 3.7467701811434450323339212228101116990

Q4,5 = −5.1042661439961787184009406231628453851

62
Q4,6 = 6.1607928057025279247501469723691945914

Q4,7 = −2.6374568452380952380952380952380952381

Q4,8 = 0.42626702123750734861845972957084068195

Q4,9 = −0.0096027579627530024355421180818006214832

Q5,1 = −0.54491666666666666666666666666666666667

Q5,2 = 3.0131226497543461829176114890400604686

Q5,3 = −6.2912305537184009406231628453850676073

Q5,4 = 5.1042661439961787184009406231628453851

Q5,6 = −3.3485729166666666666666666666666666667

Q5,7 = 2.7798440990593768371546149323927101705

Q5,8 = −0.75580278604497354497354497354497354498

Q5,9 = 0.043290030286806080456874107667758461409

Q6,1 = 0.12395833333333333333333333333333333333

Q6,2 = −1.2596013387608129671621735113798605862

Q6,3 = 4.1741890294312169312169312169312169312

Q6,4 = −6.1607928057025279247501469723691945914

Q6,5 = 3.3485729166666666666666666666666666667

Q6,7 = −0.66288489767048794826572604350382128160

Q6,8 = 0.49191676771017048794826572604350382128

Q6,9 = −0.051786576436130007558578987150415721844

Q6,10 = −0.0035714285714285714285714285714285714286

Q7,2 = 0.21072090275153271185017216763248509280

Q7,3 = −1.1810142287624926513815402704291593181

Q7,4 = 2.6374568452380952380952380952380952381

Q7,5 = −2.7798440990593768371546149323927101705

Q7,6 = 0.66288489767048794826572604350382128160

Q7,8 = 0.53307226447677836566725455614344503233

Q7,9 = −0.11780039183883429915175946921978668010

Q7,10 = 0.038095238095238095238095238095238095238

Q7,11 = −0.0035714285714285714285714285714285714286

Q8,3 = 0.090505952380952380952380952380952380952

Q8,4 = −0.42626702123750734861845972957084068195

Q8,5 = 0.75580278604497354497354497354497354498

Q8,6 = −0.49191676771017048794826572604350382128

Q8,7 = −0.53307226447677836566725455614344503233

Q8,9 = 0.77042350547472075249853027630805408583

Q8,10 = −0.20000000000000000000000000000000000000

Q8,11 = 0.038095238095238095238095238095238095238

Q8,12 = −0.0035714285714285714285714285714285714286

Q9,4 = 0.0096027579627530024355421180818006214832

Q9,5 = −0.043290030286806080456874107667758461409

Q9,6 = 0.051786576436130007558578987150415721844

Q9,7 = 0.11780039183883429915175946921978668010

Q9,8 = −0.77042350547472075249853027630805408583

Q9,10 = 0.80000000000000000000000000000000000000

63
Q9,11 = −0.20000000000000000000000000000000000000

Q9,12 = 0.038095238095238095238095238095238095238

Q9,13 = −0.0035714285714285714285714285714285714286

Q10,6 = 0.0035714285714285714285714285714285714286

Q10,7 = −0.038095238095238095238095238095238095238

Q10,8 = 0.20000000000000000000000000000000000000

Q10,9 = −0.80000000000000000000000000000000000000

Q10,11 = 0.80000000000000000000000000000000000000

Q10,12 = −0.20000000000000000000000000000000000000

Q10,13 = 0.038095238095238095238095238095238095238

Q10,14 = −0.0035714285714285714285714285714285714286

The (4-8-4) viscous terms.

The R ̄

i,j ∗ dx nonzero terms are:

R ̄

0,1 = −2.4500000000000000000000000000000000000

R ̄

0,2 = 6.0000000000000000000000000000000000000

R ̄

0,3 = −7.5000000000000000000000000000000000000

R ̄

0,4 = 6.6666666666666666666666666666666666667

R ̄

0,5 = −3.7500000000000000000000000000000000000

R ̄

0,6 = 1.2000000000000000000000000000000000000

R ̄

0,7 = −0.16666666666666666666666666666666666667

R ̄

1,1 = −1.1602083333333333333333333333333333333

R ̄

1,2 = 1.1564166666666666666666666666666666667

R ̄

1,3 = 0.34108333333333333333333333333333333333

R ̄

1,4 = −0.56016666666666666666666666666666666667

R ̄

1,5 = 0.27245833333333333333333333333333333333

R ̄

1,6 = −0.049583333333333333333333333333333333333

R ̄

2,1 = 0.093973038486880787332060994325567220421

R ̄

2,2 = −0.78549893930185917599538939754826348881

R ̄

2,3 = 0.046946483020436971650521028934066801314

R ̄

2,4 = 0.90729349114883760012958720175606131931

R ̄

2,5 = −0.26981952253301897863082607780986533944

R ̄

2,6 = −0.0063708948742115737446125435968063627529

R ̄

2,7 = 0.013476344052934369258658793939239849966

R ̄

3,1 = 0.078746435191005348453800791189557639584

R ̄

3,2 = −0.35162113760879514268346232087309277418

R ̄

3,3 = −0.78713436791133549855006393774752877572

R ̄

3,4 = 1.3893360044549165552850116426510172239

R ̄

3,5 = −0.45575020893114105530342782230270280369

R ̄

3,6 = 0.19808018162235668097172402106009418736

R ̄

3,7 = −0.092080135398545507498723070166851812875

R ̄

3,8 = 0.020423228581538619325140696189507115568

R ̄

4,1 = 0.032926148937526826310667271834304379849

R ̄

4,2 = −0.25191802597411192851286613885119855940

R ̄

4,3 = 0.86404038300564659933759665503107585561

64
R ̄

4,4 = −2.2705734372043358134175184621100160529

R ̄

4,5 = 1.9056036402801850141544340055650364152

R ̄

4,6 = −0.43544242414791128227318211540346767698

R ̄

4,7 = 0.20103197771146276694189556489629198797

R ̄

4,8 = −0.049509365793563383515243628194746597938

R ̄

4,9 = 0.0038411031851012009742168472327202485933

R ̄

5,1 = −0.049583333333333333333333333333333333333

R ̄

5,2 = 0.24538547964509061280525769040437984250

R ̄

5,3 = −0.46041931266967499182355759462391059263

R ̄

5,4 = 0.46614273043637761869004215693100311426

R ̄

5,5 = −1.3962385850218969553701412457089944128

R ̄

5,6 = 1.3149492005817200178282234837505867257

R ̄

5,7 = −0.13710317505277669221310366487172969582

R ̄

5,8 = 0.025083621435175386943099478252132619196

R ̄

5,9 = −0.0082166260206816635264869708001342671278

R ̄

6,2 = 0.013476344052934369258658793939239849966

R ̄

6,3 = −0.095466650296571124068030350198485565808

R ̄

6,4 = 0.28303601526589854838736615694592248492

R ̄

6,5 = −0.37189053652589385913924939420370302618

R ̄

6,6 = −0.87133674643234682661783233764660965659

R ̄

6,7 = 1.1363503084863723276987369694549351287

R ̄

6,8 = −0.10962852877844960704052530455432144803

R ̄

6,9 = 0.017245508513770457235161180548736518789

R ̄

6,10 = −0.0017857142857142857142857142857142857143

R ̄

7,3 = 0.020423228581538619325140696189507115568

R ̄

7,4 = −0.11461749320859847500014290897545707151

R ̄

7,5 = 0.24553622332394432725008432375299558042

R ̄

7,6 = −0.14002362171685620274536032506670231808

R ̄

7,7 = −1.2086425337884685391310360256213220479

R ̄

7,8 = 1.3257988165251715762302257817106939321

R ̄

7,9 = −0.15030001654212813132573693881511201603

R ̄

7,10 = 0.023611111111111111111111111111111111111

R ̄

7,11 = −0.0017857142857142857142857142857142857143

R ̄

8,4 = 0.0038411031851012009742168472327202485933

R ̄

8,5 = −0.021645015143403040228437053833879230705

R ̄

8,6 = 0.036310098576467624086671705719324766944

R ̄

8,7 = 0.094189280727723188040648358108675568993

R ̄

8,8 = −1.3644581220605526161081716637272192828

R ̄

8,9 = 1.4063261467781557067271352985638699924

R ̄

8,10 = −0.17638888888888888888888888888888888889

R ̄

8,11 = 0.023611111111111111111111111111111111111

R ̄

8,12 = −0.0017857142857142857142857142857142857143

R ̄

9,6 = 0.0017857142857142857142857142857142857143

R ̄

9,7 = −0.023611111111111111111111111111111111111

R ̄

9,8 = 0.17638888888888888888888888888888888889

65
R ̄

9,9 = −1.4236111111111111111111111111111111111

R ̄

9,10 = 1.4236111111111111111111111111111111111

R ̄

9,11 = −0.17638888888888888888888888888888888889

R ̄

9,12 = 0.023611111111111111111111111111111111111

R ̄

9,13 = −0.0017857142857142857142857142857142857143

The Ri,j ∗ dx nonzero terms are:

R1,1 = −1.1602083333333333333333333333333333333

R1,2 = 1.1564166666666666666666666666666666667

R1,3 = 0.34108333333333333333333333333333333333

R1,4 = −0.56016666666666666666666666666666666667

R1,5 = 0.27245833333333333333333333333333333333

R1,6 = −0.049583333333333333333333333333333333333

R2,1 = 1.2541813718202141206653943276589005538

R2,2 = −1.9419156059685258426620560642149301555

R2,3 = −0.29413685031289636168281230439926653202

R2,4 = 1.4674601578155042667962538684227279860

R2,5 = −0.54227785586635231196415941114319867278

R2,6 = 0.043212438459121759588720789736526970580

R2,7 = 0.013476344052934369258658793939239849966

R3,1 = −0.015226603295875438878260203136009580837

R3,2 = 0.43387780169306403331192707667517071464

R3,3 = −0.83408085093177247020058496668159557703

R3,4 = 0.48204251330607895515542444089495590463

R3,5 = −0.18593068639812207667260174449283746424

R3,6 = 0.20445107649656825471633656465690055011

R3,7 = −0.10555647945147987675738186410609166284

R3,8 = 0.020423228581538619325140696189507115568

R4,1 = −0.045820286253478522143133519355253259735

R4,2 = 0.099703111634683214170596182021894214773

R4,3 = 1.6511747509169820978876605927786046313

R4,4 = −3.6599094416592523687025301047610332768

R4,5 = 2.3613538492113260694578618278677392189

R4,6 = −0.63352260577026796324490613646356186434

R4,7 = 0.29311211311000827444061863506314380084

R4,8 = −0.069932594375102002840384324384253713506

R4,9 = 0.0038411031851012009742168472327202485933

R5,1 = −0.082509482270860159644000605167637713182

R5,2 = 0.49730350561920254131812382925557840190

R5,3 = −1.3244596956753215911611542496549864482

R5,4 = 2.7367161676407134321075606190410191672

R5,5 = −3.3018422253020819695245752512740308280

R5,6 = 1.7503916247296313001014055991540544027

R5,7 = −0.33813515276423945915499922976802168379

R5,8 = 0.074592987228738770458343106446879217134

66
R5,9 = −0.012057729205782864500703818032854515721

R6,1 = 0.049583333333333333333333333333333333333

R6,2 = −0.23190913559215624354659889646513999253

R6,3 = 0.36495266237310386775552724442542502682

R6,4 = −0.18310671517047907030267599998508062934

R6,5 = 1.0243480484960030962308918515052913866

R6,6 = −2.1862859470140668444460558213971963823

R6,7 = 1.2734534835391490199118406343266648245

R6,8 = −0.13471215021362499398362478280645406723

R6,9 = 0.025462134534452120761648151348870785916

R6,10 = −0.0017857142857142857142857142857142857143

R7,2 = −0.013476344052934369258658793939239849966

R7,3 = 0.11588987887810974339317104638799268138

R7,4 = −0.39765350847449702338750906592137955643

R7,5 = 0.61742675984983818638933371795669860660

R7,6 = 0.73131312471549062387247201257990733851

R7,7 = −2.3449928422748408668297729950762571765

R7,8 = 1.4354273453036211832707510862650153801

R7,9 = −0.16754552505589858856089811936384853482

R7,10 = 0.025396825396825396825396825396825396825

R7,11 = −0.0017857142857142857142857142857142857143

R8,3 = −0.020423228581538619325140696189507115568

R8,4 = 0.11845859639369967597435975620817732010

R8,5 = −0.26718123846734736747852137758687481113

R8,6 = 0.17633372029332382683203203078602708503

R8,7 = 1.3028318145161917271716843837299976169

R8,8 = −2.6902569385857241923383974454379132149

R8,9 = 1.5566261633202838380528722373789820085

R8,10 = −0.20000000000000000000000000000000000000

R8,11 = 0.025396825396825396825396825396825396825

R8,12 = −0.0017857142857142857142857142857142857143

R9,4 = −0.0038411031851012009742168472327202485933

R9,5 = 0.021645015143403040228437053833879230705

R9,6 = −0.034524384290753338372385991433610481230

R9,7 = −0.11780039183883429915175946921978668010

R9,8 = 1.5408470109494415049970605526161081717

R9,9 = −2.8299372578892668178382464096749811035

R9,10 = 1.6000000000000000000000000000000000000

R9,11 = −0.20000000000000000000000000000000000000

R9,12 = 0.025396825396825396825396825396825396825

R9,13 = −0.0017857142857142857142857142857142857143

R10,6 = −0.0017857142857142857142857142857142857143

R10,7 = 0.025396825396825396825396825396825396825

R10,8 = −0.20000000000000000000000000000000000000

R10,9 = 1.6000000000000000000000000000000000000

67
R10,10 = −2.8472222222222222222222222222222222222

R10,11 = 1.6000000000000000000000000000000000000

R10,12 = −0.20000000000000000000000000000000000000

R10,13 = 0.025396825396825396825396825396825396825

R10,14 = −0.0017857142857142857142857142857142857143

The BSi,j nonzero terms are:

BS1,1 = 2.4500000000000000000000000000000000000

BS1,2 = −6.0000000000000000000000000000000000000

BS1,3 = 7.5000000000000000000000000000000000000

BS1,4 = −6.6666666666666666666666666666666666667

BS1,5 = 3.7500000000000000000000000000000000000

BS1,6 = −1.2000000000000000000000000000000000000

BS1,7 = 0.16666666666666666666666666666666666667

68