∎

11institutetext: Peter Frolkovič 22institutetext: Svetlana Krišková
33institutetext: Michaela Rohová 44institutetext: Michal Žeravý
55institutetext: Dept. of Mathematics and Descriptive Geometry, Faculty of
Engineering STU
Radlinského 11, 810 05 Bratislava, Slovakia
55email<EMAIL_ADDRESS>
# Semi-implicit methods for advection equations with explicit forms of
numerical solution. ††thanks: The work was supported by the grant VEGA
1/0709/19 and APVV-19-0460.

Peter Frolkovič  Svetlana Krišková  Michaela Rohová  Michal Žeravý

(Received: / Accepted: The correct dates will be entered by the editor)

###### Abstract

We present a parametric family of semi-implicit second order accurate
numerical methods for non-conservative and conservative advection equation for
which the numerical solutions can be obtained in a fixed number of forward and
backward alternating substitutions. The methods use a novel combination of
implicit and explicit time discretizations for one-dimensional case and the
Strang splitting method in several dimensional case. The methods are described
for advection equations with a continuous variable velocity that can change
its sign inside of computational domain. The methods are unconditionally
stable in the non-conservative case for variable velocity and for variable
numerical parameter. Several numerical experiments confirm the advantages of
presented methods including an involvement of differential programming to find
optimized values of the variable numerical parameter.

###### Keywords:

advection equation semi-implicit method unconditional stability conservation
laws

††journal: JJIAM

##  1 Introduction

Implicit (or semi-implicit) numerical schemes are useful numerical methods to
solve advection dominated problems in several circumstances arbogast2020third
; carciopolo2019conservative ; frolkovivc2016numerical ; frolkovivc2018semi ;
fuhrmann2001stability ; hadjimichael2021positivity ; hahn2019iterative ;
ibolya2020numerical ; knodel2020global ; mikula2014inflow ; partl2016numerical
; polivka2014compositional ; puppo2021quinpi ; qin2018implicit ;
radu2010analysis ; zhang2021numerical . They can avoid or reduce significantly
the main disadvantage of fully explicit schemes that are implemented on a
fixed mesh with a finite stencil in numerical discretization. Such explicit
schemes based on finite difference or finite volume methods require a
stability restriction on the choice of discretization parameters
leveque2002finite , e.g. the time step. This can be disadvantageous in several
situations like highly variable velocity field or nonuniform space
discretization step frolkovivc2018semi ; may2017explicit  or long time
simulations when approaching stationary solutions qin2018implicit . Another
important type of problems when implicit schemes can be useful are stiff
differential equations involving the advection term frolkovivc2016numerical ;
izzo2017highly ; knodel2020global ; partl2016numerical .

In an ideal case, the implicit and semi-implicit methods can offer an
unconditional stability that make them convenient tool to solve numerically
the problems having previously mentioned complications. The price to pay is
that the numerical solution must be obtained in general by solving linear
algebraic system of equations. Our aim is to propose a semi-implicit second
order accurate numerical method having unconditional stability where the
obtained linear systems can be solved in a small given number of alternating
substitutions. Consequently, such scheme can compete well with fully explicit
schemes that themselves lack the advantage of unconditional stability.

In frolkovivc2018semi  a novel parametric family of semi-implicit second order
accurate numerical schemes for linear non-conservative advection is
introduced. We base our work here on this result that is extended in several
aspects. We treat here only representative advection equations, but the
derived numerical schemes can be used for more involved advection dominated
problems.

Additionally to frolkovivc2018semi , we show how one-dimensional discretized
problems with a velocity that can change its sign several times in domain can
be solved by one forward and one backward substitution as known for so called
fast sweeping methods lozano2021implicit ; zhao2005fast . The derivation of
second order scheme is based here on the discretization of an error term for
the fully implicit first order accurate upwind scheme. Consequently, it is
easy to apply some limiter procedures based on a blending between these two
types of schemes if necessary. The free parameter of this family of schemes is
introduced in a more convenient way than in frolkovivc2018semi . The parameter
can vary in space and time that we demonstrate in one numerical experiment by
letting find its optimized values using differential programming tools.

In this work we present for the first time the semi-implicit scheme also for
the linear advection equation in the conservative form. The scheme is based on
finite volume method and it offers an exact local mass balance property in a
discrete form that is confirmed also by numerical experiments. Nevertheless,
to insure such property, the local variability of the parameter in the method
is not allowed and the unconditional von Neumann stability can be shown only
for constant velocity case. Therefore, if the advection dominated problem can
be expressed in a non-conservative form, e.g. for the divergence free
velocity, we prefer the more flexible semi-implicit non-conservative method.
As confirmed by several numerical experiments, the errors in mass balance for
chosen numerical examples diminish fast with a mesh refinement.

Finally, we apply the derived one-dimensional non-conservative scheme for two-
dimensional advection equation using Strang splitting leveque2002finite ;
uccar2019operator  that preserves the second order accuracy. Such treatment
requires to solve only a given number of one-dimensional problems in
alternating directions that can be solved by our proposed numerical scheme.
One can show that such splitting scheme preserves the property of
unconditional stability uccar2019operator  that is not the case of unsplitted
version presented in frolkovivc2018semi .

We proceed as follows. In Section 2 we derive the parametric family of one-
dimensional semi-implicit non-conservative advection scheme. In Section 3 we
derive analogous conservative advection schemes. In Section 4 we introduce the
Strang splitting for problems in two dimension that preserves the accuracy and
the unconditional stability of 1D scheme. In Section 5 we illustrate all
properties on several numerical examples.

##  2 Non-conservative advection equation

The linear advection equation in the non-conservative form is given by

$\partial_{t}\phi+v\partial_{x}\phi=0\,,\quad\phi(x,0)=\phi^{0}(x)\,,$ (1)

where the velocity function $v=v(x)$ is a given continuous function. The
unknown function $\phi=\phi(x,t)$ is prescribed at $t=0$ by the given function
$\phi^{0}=\phi^{0}(x)$ and $\phi(x,t)$ should be determined for $t>0$ and
$x\in(0,L)\subset R$. The boundary values are prescribed by given functions
$\phi_{0}=\phi_{0}(t)$ and $\phi_{L}=\phi_{L}(t)$ only if an inflow situation
occurs at the boundary, namely,

$\phi(0,t)=\phi_{0}(t)\,,\,\,\hbox{if }v(0)\geq
0\,,\quad\phi(L,t)=\phi_{L}(t)\,,\,\,\hbox{if }v(L)\leq 0\,.$ (2)

The model equation (1) can be solved by finding so called characteristic
curves, the characteristics, generated by $v$ and using the fact that the
solution $\phi$ is constant along characteristics leveque2002finite .

In what follows, we use the following common notations of finite difference
methods. We denote $x_{i}=ih$, $i=0,1,\ldots,I$ for a chosen $I$, where
$h=L/I$, and $t^{n}=n\tau$, $n=0,1,\ldots,N$ for a chosen $N$ and $\tau>0$.
Our aim is to find the approximations $\phi_{i}^{n}$ of $\phi(x_{i},t^{n})$.
The initial values are given by $\phi_{i}^{0}=\phi^{0}(x_{i})$, and for $i=0$
and $i=I$ they can be replaced by the values using the boundary conditions
(2).

###### Remark 1

To simplify our presentation, we consider a particular form of the velocity
$v$ being a piecewise linear function of $x$ that can be determined by its
discrete values $v_{i}=v(x_{i})$ for each mesh. Such situation is typical if
the velocity field is obtained by a numerical solution of some flow equation
or when the velocity function is approximated from its values in mesh points.
We allow that the velocity can change its sign inside of $(0,L)$ when the
points for which $v=0$ can be then easily determined from its piecewise linear
form.

Later we distinguish the case of so called expanding characteristics for
points $\bar{x}\in(x_{i},x_{i+1})$ such that $v(\bar{x})=0$ and
$v^{\prime}(\bar{x})>0$. In this case, the regions $x\leq\bar{x}$ and
$x\geq\bar{x}$ are ”separated” from each other that may require a special
treatment in numerical schemes, see later. Note that if the velocity $v$ in
(1) depends also on $t$, we use an approximation by fixing its value on each
time subinterval $[t^{n},t^{n+1}]$ using a representative value, e.g., at
$t=t^{n}+\tau/2$. In such case, the positions of zero points $\bar{x}$ have to
be redetermined in each time subinterval.

###  2.1 First order accurate fully implicit scheme

The scheme can be derived using the backward Euler method for the time
discretization and the upwind one-sided finite difference for the space
discretization,

$\frac{\phi_{i}^{n+1}-\phi_{i}^{n}}{\tau}+v_{i}^{+}\,\frac{\phi_{i}^{n+1}-\phi_{i-1}^{n+1}}{h}+v_{i}^{-}\,\frac{\phi_{i+1}^{n+1}-\phi_{i}^{n+1}}{h}=0,$

where $v^{+}:=\max\\{0,v\\}$ and $v^{-}:=\min\\{0,v\\}$. Denoting (signed)
Courant numbers at mesh points

$C_{i}:=\frac{\tau v_{i}}{h}\,,$ (3)

we can write the scheme in the concise form

$\phi_{i}^{n+1}=\left(1+|C_{i}|\right)^{-1}\left(\phi_{i}^{n}+|C_{i}|\phi_{i\mp
1}^{n+1}\right)$ (4)

with $\mp=-\operatorname{sgn}(C_{i})$. The scheme (4) can be used for
$n=0,1,\ldots,N$ and for the inner nodes with $i=1,2,\ldots,I-1$. For boundary
nodes one applies either the Dirichlet boundary conditions (2),

$\displaystyle\phi_{0}^{n+1}=\phi_{0}(t^{n+1})\,,\,\,\hbox{if }C_{0}\geq 0\,,$
(5) $\displaystyle\phi_{I}^{n+1}=\phi_{L}(t^{n+1})\,,\,\,\hbox{if }C_{I}\leq
0\,,$ (6)

or one simplifies (4) for the outflow boundaries to

$\displaystyle\phi_{0}^{n+1}=\left(1-C_{0}\right)^{-1}\left(\phi_{0}^{n}-C_{0}\phi_{1}^{n+1}\right)\,,\,\,\hbox{if
}C_{0}<0\,,$ (7)
$\displaystyle\phi_{I}^{n+1}=\left(1+C_{I}\right)^{-1}\left(\phi_{I}^{n}+C_{I}\phi_{I-1}^{n+1}\right)\,,\,\,\hbox{if
}C_{I}>0\,.$ (8)

Finally, we modify (4) for the special case of expanding characteristics as
described in Remark 1. Namely, if there exists $i^{*}$ such that

$\displaystyle v_{i^{*}}<0\hbox{ and }v_{1+i^{*}}>0\,,$

then there is a point $\bar{x}_{i^{*}}\in(x_{i^{*}},x_{1+i^{*}})$ where the
piecewise linear velocity attains the zero value. For such point with
expanding characteristics we decouple the computations of $\phi_{i^{*}}^{n+1}$
and $\phi_{1+i^{*}}^{n+1}$ from each other. Namely, the scheme (4) is replaced
for $i=i^{*}$ and $i=1+i^{*}$ by the explicit definitions

$\displaystyle\begin{array}[]{l}\phi_{i^{*}}^{n+1}=\left(1-C_{i^{*}}\right)^{-1}\left(\phi_{i^{*}}^{n}-C_{i^{*}}\phi_{i^{**}}^{n}\right)\,,\\\\[4.30554pt]
\phi_{1+i^{*}}^{n+1}=\left(1+C_{1+i^{*}}\right)^{-1}\left(\phi_{1+i^{*}}^{n}+C_{1+i^{*}}\phi_{i^{**}}^{n}\right)\,.\end{array}$
(11)

To use it, we determine the location of $\bar{x}_{i^{*}}$ from the linear
interpolation,

$\bar{x}_{i^{*}}=x_{i^{*}}+\frac{v_{i^{*}}}{v_{i^{*}}-v_{1+i^{*}}}h\,.$

Afterwards, we define the interpolated value
$\phi_{i^{**}}^{n}\approx\phi(\bar{x}_{i^{*}},t^{n})$ by

$\phi_{i^{**}}^{n}:=\phi_{i^{*}}^{n}+\frac{\bar{x}_{i^{*}}-x_{i^{*}}}{h}(\phi_{1+i^{*}}^{n}-\phi_{i^{*}}^{n})$

and

$C_{1+i^{*}}:=\frac{\tau v_{1+i^{*}}}{x_{1+i^{*}}-\bar{x}_{i^{*}}}\,,\quad
C_{i^{*}}:=\frac{\tau v_{i^{*}}}{\bar{x}_{i^{*}}-x_{i}}\,.$

The schemes (11) is then determined using the approximation that
$\phi(x_{i^{*}},t)\equiv\phi_{i^{**}}^{n}$ for $t\in[t^{n},t^{n+1}]$ as
$v(\bar{x}_{i^{*}})=0$.

Formally, the schemes (4) - (11) represent a system of linear algebraic
equations. Nevertheless, each $i$-th equation in (4) - (11) contains at most
one neighbor value, either $\phi_{i-1}^{n+1}$ or $\phi_{i+1}^{n+1}$ (or none).
Such linear equations can be solved using one forward and backward
substitution as known for fast sweeping methods zhao2005fast .

Alternatively, one can formulate the substitutions analogously to
lozano2021implicit  when the discrete equations are split into two fractional
steps that must be applied in the defined order for $i$, namely

$\displaystyle\phi_{i}^{n+1/2}=\left(1+C_{i}^{+}\right)^{-1}\left(\phi_{i}^{n}+C_{i}^{+}\phi_{i-1}^{n+1}\right)\,,\,\,i=0,1,\ldots,I$
(12)
$\displaystyle\phi_{i}^{n+1}=\left(1-C_{i}^{-}\right)^{-1}\left(\phi_{i}^{n+1/2}-C_{i}^{-}\phi_{i+1}^{n+1}\right)\,,\,\,i=I,I-1,\ldots,0\,.$
(13)

The boundary conditions have to be appropriately included, namely taking (5)
instead of (12) for $i=0$ if $C_{0}\geq 0$, and (6) instead of (13) for $i=I$
if $C_{I}\leq 0$.

We can now summarize the advantages of the first order accurate implicit
upwind scheme. The numerical solution is defined by (12) - (13) explicitly as
a convex combinations of the neighbor values for arbitrary large Courant
numbers. Therefore, the scheme is unconditionally stable and it insures a
discrete minimum and maximum principle for any choice of $h$ and $\tau$. This
is especially critical if (11) has to be used as the corresponding Courant
numbers $C_{i^{*}}$ or $C_{1+i^{*}}$ can be arbitrary large.

The main disadvantage is the low accuracy of the method that can be
demonstrated for many test examples of practical interest. Therefore, we
derive an extension of this method in a form of the second order accurate
semi-implicit methods for which the numerical solution can be again determined
explicitly.

###  2.2 Parametric family of second order accurate semi-implicit schemes

We begin by presenting an error term of the first order accurate scheme (4).
To do so, we express the values $\phi(x_{i},t^{n})$ and $\phi(x_{i\pm
1},t^{n+1})$ using Taylor series at $(x_{i},t^{n+1})$. First, we have for
$t=t^{n+1}$

$\displaystyle\phi(x_{i\pm 1},\cdot)=\phi(x_{i},\cdot)\pm
h\partial_{x}\phi(x_{i},\cdot)+\frac{h^{2}}{2}\partial_{xx}\phi(x_{i},\cdot)+\mathcal{O}(h^{3})\,.$
(14)

Second, we have for $x=x_{i}$

$\displaystyle\phi(\cdot,t^{n})=\phi(\cdot,t^{n+1})-\tau\partial_{t}\phi(\cdot,t^{n+1})+\frac{\tau^{2}}{2}\partial_{tt}\phi(\cdot,t^{n+1})+\mathcal{O}(\tau^{3})=\quad$
(15) $\displaystyle\phi(\cdot,t^{n+1})+\tau
v_{i}\partial_{x}\phi(\cdot,t^{n+1})-\frac{\tau^{2}}{2}v_{i}\partial_{tx}\phi(\cdot,t^{n+1})+\mathcal{O}(\tau^{3})\,,$

where one exploits that $\partial_{t}\phi=-v\partial_{x}\phi$ and
$\partial_{tt}\phi=-v\partial_{tx}\phi)$. Such approach is often called Lax-
Wendroff procedure leveque2002finite . Opposite to its standard form when all
time derivatives in (15) are replaced by spatial derivatives using the
equation (1), we allow also the mixed derivatives in (15).

Using (15) and (14) we obtain that the discrete values $\phi(x_{i},t^{n+1})$
of exact solution fulfill the first order accurate scheme (4) with the leading
error term

$e_{2}(x_{i},t^{n+1}):=-\frac{\tau^{2}}{2}v_{i}\partial_{tx}\phi(x_{i},t^{n+1})+\frac{h\tau}{2}|v_{i}|\partial_{xx}\phi(x_{i},t^{n+1})\,.$
(16)

Now, to extend the scheme (4) to be second order accurate, we have to
approximate the derivatives in $e_{2}$ with at least first order accurate
approximations. Our aim is to derive a parametric family of semi-implicit
upwind schemes that has a convenient stencil in its implicit part and that are
unconditionally stable frolkovivc2018semi .

To do so we define the parametric “upwind based” approximations of
$\partial_{x}\phi(x_{i},t^{*})$ for $*=n$ or $*=n+1$,

$\begin{array}[]{l}h\,\partial_{x}^{\alpha-}\phi_{i}^{*}:=\alpha(\phi_{i}^{*}-\phi_{i-1}^{*})+(1-\alpha)(\phi_{i+1}^{*}-\phi_{i}^{*})\\\\[4.30554pt]
h\,\partial_{x}^{\alpha+}\phi_{i}^{*}:=\alpha(\phi_{i+1}^{*}-\phi_{i}^{*})+(1-\alpha)(\phi_{i}^{*}-\phi_{i-1}^{*})\,.\end{array}$
(17)

The approximation $\partial_{x}^{\alpha-}\phi_{i}^{*}$ will be used in (16) if
$v_{i}>0$ and $\partial_{x}^{\alpha+}\phi_{i}^{*}$ if $v_{i}<0$. Additionally,
applying the standard backward finite difference for the time derivative, and
the upwind finite difference for the spatial derivatives, we obtain

$\displaystyle\begin{array}[]{l}e_{2}(x_{i},t^{n+1})\approx-\dfrac{\tau}{2}(v_{i}^{+}\partial_{x}^{\alpha-}+v_{i}^{-}\partial_{x}^{\alpha+})(\phi_{i}^{n+1}-\phi_{i}^{n})\\\\[4.30554pt]
+\dfrac{\tau}{2}v_{i}^{+}\partial_{x}^{\alpha-}\left(\phi_{i}^{n+1}-\phi_{i-1}^{n+1}\right)-\dfrac{\tau}{2}v_{i}^{-}\partial_{x}^{\alpha+}\left(\phi_{i+1}^{n+1}-\phi_{i}^{n+1}\right)\,.\end{array}$
(20)

Doing it this way we see that the terms with
$\partial_{x}^{\alpha\mp}\phi_{i}^{n+1}$ cancel.

Combining now the first order and the second order accurate approximations, we
obtain the semi-implicit scheme of the form

$\displaystyle\phi_{i}^{n+1}+C_{i}^{+}\left(\phi_{i}^{n+1}-\phi_{i-1}^{n+1}-\frac{h}{2}\partial_{x}^{\alpha-}\phi_{i-1}^{n+1}\right)+$
(21) $\displaystyle
C_{i}^{-}\left(\phi_{i+1}^{n+1}-\phi_{i}^{n+1}-\frac{h}{2}\partial_{x}^{\alpha+}\phi_{i+1}^{n+1}\right)=\phi_{i}^{n}-C_{i}^{+}\frac{h}{2}\partial_{x}^{\alpha-}\phi_{i}^{n}-C_{i}^{-}\frac{h}{2}\partial_{x}^{\alpha+}\phi_{i}^{n}\,,$

where one can clearly distinguish between the contributions of two
approximations having the different order of accuracy. One can write (21) in
the concise form using $\mp=-\operatorname{sgn}(C_{i})$,

$\displaystyle\phi_{i}^{n+1}+|C_{i}|\left(\phi_{i}^{n+1}-\phi_{i\mp
1}^{n+1}\mp\frac{h}{2}\partial_{x}^{\alpha\mp}\phi_{i\mp
1}^{n+1}\right)=\phi_{i}^{n}\mp|C_{i}|\frac{h}{2}\partial_{x}^{\alpha\mp}\phi_{i}^{n}\,.$
(22)

In frolkovivc2018semi , using a different parametric form, it is shown that
the scheme is unconditionally stable for $\alpha\geq 0$ according to the von
Neumann stability analysis. The value $\alpha=0$ choose “downwind” one-sided
finite difference in (17) and $\alpha=1$ the upwind one. The case $\alpha=0.5$
results in the central finite difference in (17). In general, the parameters
$\alpha$ in (17) can be different for each $i$ and $n$ that we do not
emphasize in its notation, but that we exploit later in numerical experiments.
In particular, the variable choice

$\alpha=\frac{2+|C_{i}|}{6}\,,$ (23)

can improve numerical results for many examples as shown in frolkovivc2018semi
and as confirmed also for other numerical experiments here. The motivation
behind (23) is that in the case of constant velocity $v$ in (1) the third
order accuracy of (21) can be proved frolkovivc2018semi .

To obtain the numerical solution $\phi_{i}^{n+1}$ from (22), we can write it
in the form

$\displaystyle\begin{array}[]{l}\phi_{i}^{n+1}=\left(1+\frac{1+\alpha}{2}|C_{i}|\right)^{-1}\Big{(}\phi_{i-1}^{n+1}\frac{1+2\alpha}{2}C_{i}^{+}-\phi_{i-2}^{n+1}\frac{\alpha}{2}C_{i}^{+}\\\\[8.61108pt]
-\,\phi_{i+1}^{n+1}\frac{1+2\alpha}{2}C_{i}^{-}+\phi_{i+2}^{n+1}\frac{\alpha}{2}C_{i}^{-}+\phi_{i}^{n}-C_{i}^{+}\frac{h}{2}\partial_{x}^{\alpha-}\phi_{i}^{n}-C_{i}^{-}\frac{h}{2}\partial_{x}^{\alpha+}\phi_{i}^{n}\Big{)}.\end{array}$
(26)

Analogously to (4), we can obtain the values $\phi_{i}^{n+1}$ from (26) using
one forward and one backward substitution, or, similarly to (12) - (13),
splitting it into two fractional steps with the prescribed order of the
evaluation,

$\displaystyle\begin{array}[]{l}\phi_{i}^{n+1/2}=\dfrac{2\phi_{i}^{n}+\phi_{i-1}^{n+1}(1+2\alpha)C_{i}^{+}-\phi_{i-2}^{n+1}\alpha
C_{i}^{+}-C_{i}^{+}h\partial_{x}^{\alpha-}\phi_{i}^{n}}{2+(1+\alpha)C_{i}^{+}},\\\\[6.45831pt]
i=0,1,\ldots,I\,,\end{array}$ (29)
$\displaystyle\begin{array}[]{l}\phi_{i}^{n+1}=\dfrac{2\phi_{i}^{n+1/2}-\phi_{i+1}^{n+1}(1+2\alpha)C_{i}^{-}+\phi_{i+2}^{n+1}\alpha
C_{i}^{+}-C_{i}^{-}h\partial_{x}^{\alpha+}\phi_{i}^{n}}{2-(1+\alpha)C_{i}^{-}},\\\\[6.45831pt]
i=I,I-1,\ldots,0\,.\end{array}$ (32)

The scheme (29) - (32) represents the explicit form of numerical solution for
the proposed family of semi-implicit non-conservative scheme to solve (1). The
boundary conditions have to be appropriately included, namely taking (5)
instead of (29) for $i=0$ if $C_{0}\geq 0$, and (6) instead of (32) for $i=I$
if $C_{I}\leq 0$. If for $i=1$ or $i=I$ one has $C_{1}>0$ or $C_{I-1}<0$ then
one can set $\alpha=0$ in (29) or (32). In the case that for some $i$ one has
$C_{i}\leq 0$ and $C_{i+1}>0$, one may use (11) instead.

##  3 Conservative advection equation

The linear advection equation in the conservative form is written as

$\partial_{t}\phi+\partial_{x}\left(v\phi\right)=0\,,\quad\phi(x,0)=\phi^{0}(x)\,.$
(33)

The same assumptions on the input functions, the initial and boundary
conditions as in the non-conservative case apply also here.

To use a conservative finite difference (or a finite volume) method, we divide
the interval $(0,L)$ to subintervals (the ”control volumes”)
$(x_{i-1/2},x_{i+1/2})$, where the ”face points” are given by $x_{i+1/2}=ih$,
$i=0,1,\ldots,I$ using the discretization step $h=L/I$. The points $x_{i}$ are
now shifted compared to the notation in Section 2, namely

$x_{i}=ih-h/2\,,\,\,i=1,2,\ldots I.$

Our aim is now to find the approximations

$\Phi_{i}^{n}\approx\frac{1}{h}\int_{x_{i-1/2}}^{x_{i+1/2}}\phi(x,t^{n})\,dx\approx\phi(x_{i},t^{n})\,,$

for $i=1,2,\ldots,I$. For the initial conditions, we consider
$\Phi_{i}^{0}=\phi^{0}(x_{i})$ that is a second order accurate approximation
of the above integrals. The velocity is evaluated in points $x_{i-1/2}$, i.e.
$v_{i+1/2}:=v(x_{i+1/2})$, in particular the boundary fluxes are given by

$v_{1/2}=v(x_{1/2})=v(0)\,,\,\,\,v_{I+1/2}=v(x_{I+1/2})=v(L)\,.$

###  3.1 First order accurate fully implicit scheme

We define the scheme in a locally conservative form

$\displaystyle\frac{\Phi_{i}^{n+1}-\Phi_{i}^{n}}{\tau}+$ (34)
$\displaystyle\frac{v_{i+1/2}^{+}\Phi_{i}^{n+1}+v_{i+1/2}^{-}\Phi_{i+1}^{n+1}-v_{i-1/2}^{+}\Phi_{i-1}^{n+1}-v_{i-1/2}^{-}\Phi_{i}^{n+1}}{h}=0\,.$

Indexing now the signed Courant numbers at the faces

$C_{i+1/2}:=\frac{\tau v_{i+1/2}}{h}\,,\,\,i=0,1,\ldots,I\,,$

we can rewrite the equations (34) for $i=1,2,\ldots,I$ to the form

$\displaystyle\Phi_{i}^{n+1}=\frac{\Phi_{i}^{n}-C_{i+1/2}^{-}\Phi_{i+1}^{n+1}+C_{i-1/2}^{+}\Phi_{i-1}^{n+1}}{1+C_{i+1/2}^{+}-C_{i-1/2}^{-}}\,,$
(35)

where we define, formally,

$\displaystyle\Phi_{0}^{n+1}=\phi_{0}(t^{n+1}),\quad\Phi_{I+1}^{n+1}=\phi_{L}(t^{n+1})\,.$

The scheme (35) defines the first order accurate conservative implicit upwind
method for (33). The scheme represents a system of linear algebraic equations
that can be solved using one forward and one backward substitution as
described in the previous section. This can be viewed as the most important
advantage of the first order accurate upwind scheme together with its locally
conservative form. It is important to note that the case of zero velocity with
diverging characteristics as described in Remark 1 is captured by the
conservative scheme (35) automatically and no special treatment is required
here.

The scheme can be written using two fractional time steps with the prescribed
order,

$\displaystyle\Phi_{i}^{n+1/2}=\left(1+C_{i+1/2}^{+}\right)^{-1}\left(\Phi_{i}^{n}+C_{i-1/2}^{+}\Phi_{i-1}^{n+1}\right)\,,\,\,i=0,1,\ldots,I$
$\displaystyle\Phi_{i}^{n+1}=\left(1-C_{i-1/2}^{-}\right)^{-1}\left(\Phi_{i}^{n+1/2}-C_{i+1/2}^{-}\Phi_{i+1}^{n+1}\right)\,,\,\,i=I,I-1,\ldots,0\,.$

The main disadvantage is again the low accuracy that motivates us to extend
the scheme in a form of second order accurate semi-implicit method.

###  3.2 Parametric class of second order accurate semi-implicit schemes

To derive the error term of the first order accurate scheme (34), we express
again the values $\phi(x_{i},t^{n})$ using Taylor series at $(x_{i},t^{n+1})$.
Instead of (15), we obtain now

$\displaystyle\phi(x_{i},t^{n})=\phi(x_{i},t^{n+1})+$
$\displaystyle\tau\partial_{x}(v(x_{i})\phi(x_{i},t^{n+1}))-\frac{\tau^{2}}{2}\partial_{tx}\left(v(x_{i})\phi(x_{i},t^{n+1})\right)+\mathcal{O}(\tau^{3})\,,$

where we exploited that $\partial_{t}\phi=-\partial_{x}(v\phi)$ and
$\partial_{tt}\phi=-\partial_{tx}(v\phi)$. Using additionally (14) and

$v_{i\pm
1/2}=v_{i}\pm\frac{h}{2}v^{\prime}_{i}\,,\quad\partial_{x}(v_{i}\phi(x_{i},t^{n+1}))=v_{i}\partial_{x}\phi(x_{i},t^{n+1})+v^{\prime}_{i}\phi(x_{i},t^{n+1})\,,$

we obtain the following form of the second order error term

$\displaystyle
e_{2}(x_{i},t^{n+1}):=-\frac{\tau^{2}}{2}\partial_{tx}(v_{i}\phi(x_{i},t^{n+1}))\pm\frac{\tau
h}{2}\partial_{x}\left(v_{i}\partial_{x}\phi(x_{i},t^{n+1})\right)\,,$

where $\pm=\operatorname{sgn}(v_{i})$. Now denoting analogously to (17)

$\Phi_{i,\alpha-}^{*}=\alpha\Phi_{i}^{*}+(1-\alpha)\Phi_{i+1}^{*}\,,\quad\Phi_{i,\alpha+}^{*}=(1-\alpha)\Phi_{i}^{*}+\alpha\Phi_{i+1}^{*}\,,\quad$

we can apply the following approximations,

$\displaystyle\partial_{x}(v(x_{i})\phi(x_{i},t^{n+1}))\approx$ $\displaystyle
v_{i+1/2}^{+}\,\Phi_{i,\alpha-}^{n+1}-v_{i-1/2}^{+}\,\Phi_{i-1,\alpha-}^{n+1}+v_{i+1/2}^{-}\,\Phi_{i,\alpha+}^{n+1}-v_{i-1/2}^{-}\,\Phi_{i-1,\alpha+}^{n+1}\,,$

and

$\displaystyle hv_{i+1/2}\partial_{x}\phi(x_{i+1/2},t^{n+1})\approx$
$\displaystyle
v_{i+1/2}^{+}(\Phi_{i,\alpha-}^{n+1}-\Phi_{i-1,\alpha-}^{n+1})+v_{i+1/2}^{-}(\Phi_{i+1,\alpha+}^{n+1}-\Phi_{i,\alpha+}^{n+1})\,.$

Together with the above approximations we use the backward finite difference
in time and the central difference for the first space derivative in the
second term of $e_{2}$. After some algebraic manipulations when several terms
cancel, we obtain

$\displaystyle e(x_{i},t^{n+1})\approx$
$\displaystyle-\frac{\tau}{2h}\left(v_{i+1/2}^{+}\left(\Phi_{i-1,\alpha-}^{n+1}-\Phi_{i,\alpha-}^{n}\right)+v_{i+1/2}^{-}\left(\Phi_{i+1,\alpha+}^{n+1}-\Phi_{i,\alpha+}^{n}\right)\right.$
$\displaystyle\left.-v_{i-1/2}^{+}\left(\Phi_{i-2,\alpha-}^{n+1}-\Phi_{i-1,\alpha-}^{n}\right)-v_{i-1/2}^{-}\left(\Phi_{i,\alpha+}^{n+1}-\Phi_{i-1,\alpha+}^{n}\right)\right)\,.$

Putting together the first order and the second order approximations, the
second order accurate semi-implicit conservative scheme can be written in the
form

$\Phi_{i}^{n+1}+\frac{\tau}{h}\left(F_{i+1/2}^{n+1/2}-F_{i-1/2}^{n+1/2}\right)=\Phi_{i}^{n}\,,$
(36)

where the numerical fluxes for $i=0,1,\ldots,I$ are defined by

$\displaystyle
F_{i+1/2}^{n+1/2}:=v_{i+1/2}^{+}\left(\Phi_{i}^{n+1}-\frac{1}{2}\Phi^{n+1}_{i-1,\alpha-}+\frac{1}{2}\Phi^{n}_{i,\alpha-}\right)+$
(37) $\displaystyle
v_{i+1/2}^{-}\left(\Phi_{i+1}^{n+1}-\frac{1}{2}\Phi^{n+1}_{i+1,\alpha+}+\frac{1}{2}\Phi^{n}_{i,\alpha+}\right)\,.$

Concerning the values $\Phi_{0}^{n+1}$ and $\Phi_{I+1}^{n+1}$ that occur in
(37) for $i=1$ if $v_{3/2}>0$, and for $i=I$ if $v_{I-1/2}<0$, we use the
linear extrapolation,

$\Phi_{0}^{n+1}=2\phi_{0}(t^{n+1})-\Phi_{1}^{n+1}\,,\quad\Phi_{I+1}^{n+1}=2\phi_{L}(t^{n+1})-\Phi_{I}^{n+1}\,.$
(38)

Analogously, if $v_{1/2}<0$ and $v_{I+1/2}>0$, we again linearly extrapolate
the missing values by

$\Phi_{0}^{n}=2\phi_{0}(t^{n})-\Phi_{1}^{n}\,,\quad\Phi_{I+1}^{n}=2\phi_{L}(t^{n})-\Phi_{I}^{n}\,.$

Concerning the inflow fluxes $F_{1/2}^{n+1/2}$ and $F_{I+1/2}^{n+1/2}$ at the
boundary, we apply a second order accurate approximation, e.g.,

$\displaystyle F_{1/2}^{n+1/2}=v_{1/2}\phi_{0}(t^{n+1/2}),\hbox{ if
}v_{1/2}\geq 0\,,\quad$ $\displaystyle
F_{I+1/2}^{n+1/2}=v_{I+1/2}\phi_{L}(t^{n+1/2}),\hbox{ if }v_{I+1/2}\leq 0\,.$

We formulate now the scheme (36) in the form suitable for the fast sweeping
method. First we substitute (37) to (36),

$\displaystyle\Phi_{i}^{n+1}+\frac{1}{2}C_{i+1/2}^{+}\left((1+\alpha)\Phi_{i}^{n+1}-\alpha\Phi^{n+1}_{i-1}+\alpha\Phi_{i}^{n}+(1-\alpha)\Phi_{i+1}^{n}\right)$
$\displaystyle+\left.\frac{1}{2}C_{i+1/2}^{-}\left((1+\alpha)\Phi_{i+1}^{n+1}-\alpha\Phi^{n+1}_{i+2}+(1-\alpha)\Phi_{i}^{n}+\alpha\Phi_{i+1}^{n}\right)\right.-$
$\displaystyle\frac{1}{2}C_{i-1/2}^{+}\left((1+\alpha)\Phi_{i-1}^{n+1}-\alpha\Phi^{n+1}_{i-2}+\alpha\Phi_{i-1}^{n}+(1-\alpha)\Phi_{i}^{n}\right)-$
$\displaystyle\frac{1}{2}C_{i-1/2}^{-}\left((1+\alpha)\Phi_{i}^{n+1}-\alpha\Phi^{n+1}_{i+1}+(1-\alpha)\Phi_{i-1}^{n}+\alpha\Phi_{i}^{n}\right)=\Phi_{i}^{n}\,.$

Collecting all terms, we can write

$\displaystyle\left(2+(1+\alpha)C_{i+1/2}^{+}-(1+\alpha)C_{i-1/2}^{-}\right)\Phi_{i}^{n+1}=2\Phi_{i}^{n}$
(39)
$\displaystyle+\left.C_{i+1/2}^{+}\left(\alpha\Phi_{i-1}^{n+1}-\alpha\Phi_{i}^{n}-(1-\alpha)\Phi_{i+1}^{n}\right)\right.$
$\displaystyle-\left.C_{i+1/2}^{-}\left((1+\alpha)\Phi_{i+1}^{n+1}-\alpha\Phi_{i+2}^{n+1}+\alpha\Phi_{i+1}^{n}+(1-\alpha)\Phi_{i}^{n}\right)\right.$
$\displaystyle+\left.C_{i-1/2}^{+}\left((1+\alpha)\Phi_{i-1}^{n+1}-\alpha\Phi_{i-2}^{n+1}+\alpha\Phi_{i-1}^{n}+(1-\alpha)\Phi_{i}^{n}\right)\right.$
$\displaystyle-
C_{i-1/2}^{-}\left(\alpha\Phi_{i+1}^{n+1}-\alpha\Phi_{i}^{n}-(1-\alpha)\Phi_{i-1}^{n}\right)\,.$

Dividing (39) by the term before $\Phi_{i}^{n+1}$ we obtain the formula to be
used with the fast sweeping method. In the case that $C_{i-1/2}<0$ and
$C_{i+1/2}>0$ one can use the first order scheme (35) instead of (39).

For $i=1$ and $i=I$ we modify the scheme according to the boundary conditions
as mentioned before, namely for $i=1$

$\displaystyle\left(2+(1+2\alpha)C_{3/2}^{+}-(1+\alpha)C_{1/2}^{-}\right)\Phi_{1}^{n+1}=2\Phi_{1}^{n}$
(40)
$\displaystyle+\left.C_{3/2}^{+}\left(2\alpha\phi_{0}(t^{n+1})-\alpha\Phi_{1}^{n}-(1-\alpha)\Phi_{2}^{n}\right)\right.$
$\displaystyle-\left.C_{3/2}^{-}\left((1+\alpha)\Phi_{2}^{n+1}-\alpha\Phi_{3}^{n+1}+\alpha\Phi_{2}^{n}+(1-\alpha)\Phi_{1}^{n}\right)\right.$
$\displaystyle+\left.C_{1/2}^{+}\phi_{0}(t^{n+1/2})-C_{1/2}^{-}\left(\alpha\Phi_{2}^{n+1}+(1-2\alpha)\Phi_{1}^{n}-2(1-\alpha)\phi_{0}(t^{n})\right)\right.\,,$

and analogously for $i=I$,

$\displaystyle\left(2+(1+\alpha)C_{I+1/2}^{+}-(1+2\alpha)C_{I-1/2}^{-}\right)\Phi_{I}^{n+1}=2\Phi_{I}^{n}\quad$
(41)
$\displaystyle+\left.C_{I+1/2}^{+}\left(\alpha\Phi_{I-1}^{n+1}+(1-2\alpha)\Phi_{I}^{n}-2(1-\alpha)\phi_{L}(t^{n})\right)-C_{I+1/2}^{-}\phi_{L}(t^{n+1/2})\right.$
$\displaystyle+\left.C_{I-1/2}^{+}\left((1+\alpha)\Phi_{I-1}^{n+1}-\alpha\Phi_{I-2}^{n+1}+\alpha\Phi_{I-1}^{n}+(1-\alpha)\Phi_{I}^{n}\right)\right.$
$\displaystyle-
C_{I-1/2}^{-}\left(2\alpha\phi_{L}(t^{n+1})-\alpha\Phi_{I}^{n}-(1-\alpha)\Phi_{I-1}^{n}\right)\,.$

Note that for $i=2$ with $C_{3/2}>0$ and $i=I-1$ with $C_{I-1/2}<0$ one has to
use also (38).

The derived scheme is exactly mass conservative at the discrete level and it
is second order accurate as it is suggested also by numerical experiments. In
the case of constant velocity $v$ in (33) the scheme is equivalent to the
parametric family of nonconservative semi-implicit scheme (26) when also the
von Neumann unconditional stability is valid.

We note that although the parameter $\alpha$ can be chosen freely in each time
interval, e.g. $\alpha\in[0,1]$, it shall not vary with respect to $i$ if the
discrete form of local mass balance property shall be fulfilled.

##  4 Two dimensional case

We now apply the so called Strang splitting to solve the non-conservative
advection equation in two-dimensional case, but the idea can be applied in
more dimensional cases and for the conservative form, too.

The advection equation is now given in the form

$\partial_{t}\phi+\vec{v}\cdot\nabla\phi=0\,,\quad\phi(x,y,0)=\phi^{0}(x,y)\,,(x,y)\in\Omega\,,$
(42)

where $\Omega\subset R^{2}$ takes here a simple form of a square and
$\vec{v}=(v_{1}(x,y),v_{2}(x,y))$ is a given velocity vector field. The
boundary conditions are defined depending on the flow regime at the boundary
$\partial\Omega$ with the values prescribed only at the inflow part,

$\phi(x,y,t)=\phi_{in}(x,y,t)\hbox{ if }\vec{n}(x,y)\cdot\vec{v}(x,y)\leq
0\,,\,\,(x,y)\in\partial\Omega\,,$ (43)

where $\vec{n}$ is the outward normal vector.

We use analogous notation to derive numerical approximations as in the
previous section with the addition that $y_{j}=jh$ for $j=0,1,\ldots,I$.
Furthermore, we denote $t^{n+1/2}=t^{n}+\tau/2$.

The idea of the time splitting method is to approximate and split the problem
(42) into two subproblems that are coupled only by the choice of initial
conditions, and that are solved separately in a specified sequence. Let us
explain the simplest variant in details.

The first subproblem takes the form of one dimensional advection equations for
the parameter $y\in(0,L)$,

$\displaystyle\partial_{t}\phi(x,y,t)+v_{1}(x,y)\partial_{x}\phi(x,y,t)=0\,,\,\,x\in(0,L)\,,$
(44)

and the second subproblem takes the analogous form for the parameter
$x\in(0,L)$,

$\displaystyle\partial_{t}\phi(x,y,t)+v_{2}(x,y)\partial_{y}\phi(x,y,t)=0\,,\,\,y\in(0,L)\,.$
(45)

The splitting in time is realized as follows. Let the solution
$\phi(x,y,t^{n})$ (or its approximation) of (42) be available at some time
$t=t^{n}$, $0\leq n<N$. To obtain an approximation of $\phi(x,y,t^{n+1})$ we
do three steps. First, the subproblem (44) is solved for $t\
\in(t^{n},t^{n+1/2})$ with the initial condition defined by $\phi(x,y,t^{n})$.
Afterwards, the subproblem (45) is solved for $t\in(t^{n},t^{n+1})$ with the
initial condition defined by the solution of (44) at $t=t^{n+1/2}$. Finally,
the subproblem (44) is solved now for $t\in(t^{n+1/2},t^{n+1})$ with the
initial condition defined by the solution of (45) at $t=t^{n+1}$. The result
is the desired approximation of $\phi(x,y,t^{n+1})$.

To discretize the advection equation (42) also in space, we consider the first
subproblem (44) only for $y=y_{j}$, $j=0,1,\ldots,J$, and the second one (45)
only for $x=x_{i}$, $i=0,1,\ldots,I$. For the resulting one-dimensional
advection problems we can apply the numerical method from the previous
section, when each of the resulting discrete algebraic systems can be solved
in one forward and one backward substitution.

##  5 Numerical experiments

In following numerical experiments we want to illustrate the properties of the
derived semi-implicit schemes using some standard benchmarks. For the non-
conservative advection, the methods are implemented in $C$-language for two-
dimensional problems. The numerical solutions of one-dimensional conservative
advection is implemented in Python. The experiments with optimized choice of
parameters using automatic differentiation is realized with Python and its
library PyTorch.

For examples having available exact solution on a whole time interval we
compute the global discrete $L_{1}$ errors $E$ that takes in one-dimensional
case the form

$E=E(h,\tau)=h\tau\sum_{i,n}|\phi_{i}^{n}-\phi(x_{i},t^{n})|$ (46)

and analogously for two-dimensional case. If an exact solution is available
only at the final time $t^{N}=T$, we compute the error in the form

$E^{N}=E^{N}(h,\tau)=h\sum_{i}|\phi_{i}^{N}-\phi(x_{i},T)|$ (47)

and analogously in two-dimensional case.

###  5.1 One dimensional nonconservative advection

Figure 1: The exact solution for $t=0$, $t=T/2$, and $t=T$ (left) and the
comparison of exact and numerical solution at $t=T$ (right) for the example
5.1.

The following example contains all important features of the non-conservative
linear advection equation (1) with variable velocity and general boundary
conditions (2). The velocity changes its sign twice inside of computational
interval with converging and expanding characteristics and it prescribes
variable inflow and outflow boundary conditions.

The example if formally treated as two-dimensional. The computational domain
is $\Omega=(-\pi/2,3\pi/2)^{2}$, the initial function $\phi(x,y,0)=\sin(x)$
and $T=1.2$. The velocity function $\vec{v}=(v_{1}(x,y),v_{2}(x,y))$ is
depending on the spatial coordinate $x$ only with $v_{1}(x,y)=\sin(x)$ and
$v_{2}(x,y)=0$. The exact solution is given by

$\phi=\sin(2\arctan(e^{-t}\tan(x/2)))\,.$

We solve the example with discretization steps resulting in the maximal
Courant number being approximately $3.81$. In Table 1 we present the global
discrete errors $E$ in (46) for two interesting choices of the parameter
$\alpha$, namely $\alpha=0.5$ and $\alpha$ defined by (23). One can see that
the results are of a good accuracy even for Courant number larger than $1$.
The EOC is approaching $2$ from above in both cases with the results for the
variable $\alpha$ slightly better than for the fixed value.

$I$ | $N$ | $E,\alpha=0.5$ | EOC |  $E,\alpha$ in (23) | EOC
---|---|---|---|---|---
40 | 1 | 0.810861 | - | 0.556925 | -
80 | 2 | 0.167179 | 2.278 | 0.099711 | 2.481
160 | 4 | 0.035211 | 2.247 | 0.018519 | 2.428
320 | 8 | 0.007858 | 2.163 | 0.003831 | 2.273
Table 1: Numerical errors and EOC for two choice of $\alpha$ in example 5.1.

Finally, we compute the example on the finest mesh using only one time step
resulting in the maximal Courant number being approximately $30.5$. The
obtained numerical solution is compared with the one obtained on the coarsest
mesh in Figure 2. One can see that the large Courant numbers do not result in
any instabilities.

Figure 2: Comparison of exact and numerical solutions obtained only with one
time step for example 5.1. The maximal Courant number is $3.81$ (left) and
$30.5$ (right).

###  5.2 Optimization of parameters by automatic differentiation

To illustrate the possibilities of variable parameter $\alpha=\alpha_{i}^{n}$
in (22), we compute an example with variable velocity, where we let the code
optimize the values of $\alpha_{i}^{n}$ at each time step and at each grid
point. To do so we use the library PyTorch NEURIPS2019_9015  to implement a
straightforward gradient descent method to minimize the loss function

$J=J(I,N)=h\tau\sum_{i,n}(\min\\{0,\phi_{i}^{n}\\})^{2}$

with respect to parameters $\alpha_{i}^{n}$.

In particular, we set initially $\alpha_{i}^{n}=0.5$ for $i=1,2,\ldots,I-1$
and $n=1,2,\ldots,N-1$ for which we compute the values $\phi_{i}^{n}$ of
numerical solution. Afterwards using the automatic differentiation available
in PyTorch we obtain the gradient of $J$ with respect to all values of
$\alpha_{i}^{n}$. Next we subtract the gradient multiplied by a parameter (the
”learning rate”) $\eta$ from the values of $\alpha_{i}^{n}$ that results in a
smaller value of $J$. In theory, one can continue with this procedure up to a
point when the decrease of $J$ is not substantial. In our case, we use only
one step of such optimization.

To show clearly this idea we choose the following example. The domain is
$\Omega=(-2,12)$, $T=2\pi/\sqrt{3}$, the initial function
$\phi(x,0)=\exp(-2x^{2})$ and the velocity function $v(x)=2+\sin(x)$. One can
show that $\phi(x,T)=\phi^{0}(x-2\pi)$.

In Table 2 we summarize the results. For three consecutively refined meshes we
present the value of $J$ before and after the optimization step, and
analogously the errors $E^{N}$ in (47). One can see that the unphysical
oscillations can be decreased significantly with a slight improvement in the
precision.

$I$ | $N$ | $\eta\cdot 10^{6}$ | $\,\,J_{b}\cdot 10^{-3}$ | $\,\,J_{a}\cdot 10^{-3}$ | $\,\,E_{b}^{N}$ | $\,\,E_{a}^{N}$
---|---|---|---|---|---|---
70 | 50 | 0.2 | 3.68 | 0.0768 | 0.521 | 0.511
140 | 100 | 8.0 | 1.12 | 0.0156 | 0.197 | 0.190
280 | 200 | 160.0 | 0.0664 | 0.00354 | 0.0533 | 0.0448
Table 2: The learning rate $\eta$, the values of the loss function $J$, and
the errors $E^{N}$ in (47) using three meshes for example 5.2. The indices $b$
and $a$ means ”before” and ”after” optimization.

To show the influence of optimized values $\alpha$ visually, we present the
figures for $I=70$ and $N=50$ representing a rather coarse time and space
discretization. In Figure 3 we present the numerical solutions at
$t^{N/2}=T/2$ and $t^{N}=T$ for the fixed choice of all $\alpha_{i}^{n}\equiv
0.5$ when one can observe clearly some unphysical negative values.
Furthermore, the numerical solutions at the same times obtained after one
optimization step are plotted together with the values of $\alpha_{i}^{N/2}$
and $\alpha_{i}^{N}$. One can clearly observe that the largest unphysical
oscillation are suppressed. Analogously, the same results are presented also
for the refined discretization steps with $I=140$ and $N=100$.

Figure 3: Numerical solution and variable parameters for $t=T/2$ (left) and
for $t=T$ (right) for example 5.2. The top row is obtained with $I=70$ and the
bottom one with $I=140$.

###  5.3 One dimensional conservative advection

The following example illustrates the applicability of semi-implicit scheme
for the advection equation in the conservative form (33). To check the mass
conservation property, we choose an example with zero velocity at boundary
points. Inside of the computational interval
$\Omega=(-\frac{\pi}{2},\frac{5\pi}{2})$ the variable velocity $v(x)=\cos(x)$
changes twice its sign. The initial condition has the form $\phi(x,0)=\cos(x)$
and $T=1$. The exact solution is given by

$\phi=\cos^{2}\left(2\tan^{-1}\left(\tanh\left(\frac{1}{2}\left(t-2\tanh^{-1}\left(\tan\left(\frac{x}{2}\right)\right)\right)\right)\right)\right)\,,$

see the left picture in Figure 4.

Figure 4: The exact solution at $t=0$, $t=T/2$, and $t=T$ (left) and the
comparison of exact and numerical solution at $t=T$ (right) for example 5.3.

We solve the example with discretization steps such that the maximal Courant
number is approximately $4.24$. The global discrete errors $E$ in (46) are
presented in Table 3 for $\alpha=0.5$ and $\alpha=1$. The numerical results
are stable even for Courant number larger than $1$, see the right picture in
Figure 4 for a visual comparison. The EOC is approximately $2$ in both cases
with the choice $\alpha=1$, as expected, being slightly better for the example
with larger Courant numbers.

$I$ | $N$ | $E,{\alpha=0.5}$ | EOC | $E,{\alpha=1}$ | EOC
---|---|---|---|---|---
40 | 1 | 0.9610 | - | 0.7013 | -
80 | 2 | 0.2750 | 1.81 | 0.1941 | 1.85
160 | 4 | 0.0651 | 2.08 | 0.0442 | 2.13
320 | 8 | 0.0150 | 2.12 | 0.0098 | 2.17
Table 3: Numerical errors for example 5.3 with the maximal Courant number $4$.

In Table 4 we present the results for the same example with the four times
smaller time step resulting in the maximal Courant number being approximately
$1$. Of course, the precision of results increases, and, moreover, the choice
$\alpha=0.5$ gives now better results than the choice $\alpha=1$.

$I$ | $N$ | $E,{\alpha=0.5}$ | EOC | $E,{\alpha=1}$ | EOC
---|---|---|---|---|---
40 | 4 | 0.1181 | - | 0.1683 | -
80 | 8 | 0.0256 | 2.20 | 0.0461 | 1.87
160 | 16 | 0.0054 | 2.26 | 0.011 | 2.02
320 | 32 | 0.0012 | 2.16 | 0.0028 | 2.02
Table 4: Numerical errors for example 5.3 with the maximal Courant number $1$.

Note that we always obtain a perfect mass conservation at the discrete level
by checking

$\sum_{i}\Phi_{i}^{n}=\sum\Phi_{i}^{0}\,,\quad n=1,2,\ldots,N$

with a difference given only by rounding errors around $10^{-15}$.

Finally, to illustrate the stability of our scheme, we compute the example on
the finest mesh with only one time step resulting in the maximal Courant
number approximately $34$, see the right picture in Figure 5 with no
instabilities occurring.

Figure 5: Numerical solutions obtained with only one time step for example
5.3. The maximal Courant number equals $4.24$ (left) and $33.95$ (right).

###  5.4 Two-dimensional example with diagonally variable velocity

With the next example we solve an analogous problem to Example 5.1 with the
velocity $\vec{v}=(\sin(2\pi(x+y)/2),\sin(2\pi(x+y)/2))$ that is variable only
in a diagonal direction of the domain $\Omega=(-1,2)\times(-1,2)$. We use the
initial condition $\phi(x,y,0)=\sin(2\pi(x+y)/2)$ and $T=0.24$.

The example is chosen to test the accuracy of the operator splitting method as
described in Section 4. To compute the global error $E$ analogous to (46) we
use always a half time step twice in the subproblem (44) when compared to the
time step of subproblem (45) as described in Section 4. Consequently, the
maximal Courant number in the $x$ direction is a half of the maximal Courant
number in $y$ direction. Note that for explicit methods a typical stability
restriction is given by the sum of directional Courant numbers that is
required to be smaller than one leveque2002finite . We compute the example
with the maximal Courant number in $y$ direction being $1.6$.

The exact solution is given by

$\phi=\sin(2\arctan(e^{-2\pi t}\tan(\pi\frac{x+y}{2}))).$

The numerical results are presented in Figure 6 and 7 and in Table 5. We see
that the method for two typical choices of parameter $\alpha$ are second order
accurate for this example.

Figure 6: The initial condition for Example 5.4 viewed from a top and a side.

Figure 7: The numerical solution for Example 5.4 in time $t=T$ for $I=160$ viewed from a top and a side. $I$ | $N$ | $E,\alpha=0$ | EOC |  $E,\alpha$ in (23) | EOC
---|---|---|---|---|---
20 | 1 | 0.0874 | - | 0.0838 | -
40 | 2 | 0.0179 | 2.29 | 0.0173 | 2.27
80 | 4 | 0.00319 | 2.49 | 0.00302 | 2.52
160 | 8 | 0.000624 | 2.36 | 0.000569 | 2.41
Table 5: Numerical errors for two choices of $\alpha$ for example 5.4.

###  5.5 Example with deformation velocity

To test our method for a nontrivial case, we choose two-dimensional example
with a deformation velocity in which the initial profile of solution is
deformed significantly in time. For this example we check quantitatively not
only the numerical errors, but also two other numerical artifacts - negative
unphysical oscillations and a violation of mass conservation in a discrete
form. As we show, the both of them are visible for a coarse mesh, but these
numerical errors decrease rapidly with the mesh refinement.

The domain is the unit square $\Omega=(0,1)\times(0,1)$ and $T=1$. The
velocity vector $\vec{v}$, see Figure 8, is defined by

$\displaystyle v_{1}(x,y,t)=-4\cos$ $\displaystyle(\pi t)$
$\displaystyle\sin^{2}(2\pi x)\sin(2\pi y)\cos(2\pi y)\,,$ $\displaystyle
v_{2}(x,y,t)=4\cos$ $\displaystyle(\pi t)$ $\displaystyle\sin^{2}(2\pi
y)\sin(2\pi x)\cos(2\pi x)\,.$

Note that the time dependency of $\vec{v}$ in numerical simulations is
resolved by evaluating $v$ for each $n$ at $t=t^{n}+\tau/2$. The velocity has
the zero divergence and it is equal zero at the boundary, so the integral of
the initial function (the mass) shall be conserved in time.

Figure 8: The velocity field $\vec{v}$ at $t=0$ for example 5.5.

We consider two different initial conditions in this example, see Figure 9.
First, the Gaussian is chosen

$\phi(x,y,0)=e^{-100\left((x-0.5)^{2}+(y-0.5)^{2}\right)}\,,$ (48)

and, second, the signed distance function is considered,

$\phi(x,y,0)=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}\,.$ (49)

Figure 9: Initial profiles (48) (left) and (49) (right) for example 5.5.

In both cases, the initial profile of the solution $\phi$ is deformed up to
time $t=0.5$ when the direction of velocity change its sign, so the
deformation is reversed afterwards and the initial profile shall be recovered
at time $t=1$.

We compute the example with maximal Courant numbers being $0.76$ in
$y$-direction. The numerical solutions for the finest mesh at the time of
maximal deformation are plotted in Figure 10 and, moreover, at the final time
when the initial profile should be reproduced, see Figure 11. In Tables 6 and
7 we compute the error $E^{N}$ analogous to (47) obtained for each numerical
solution as the difference between its values at the initial and the final
time.

Figure 10: The maximal deformation at $t=0.5$ for (48) (left) and (49) (right)
for example 5.5.

Figure 11: Solution at $t=1$ and $I=320$ for (48) (left) and (49) (right) for example 5.5 $I$ | $\Delta t$ | $E,\alpha=0.5$ | EOC |  $E,\alpha$ in (23) | EOC
---|---|---|---|---|---
40 | 100 | 0.01088 | - | 0.00928 | -
80 | 200 | 0.00507 | 1.10 | 0.00415 | 1.16
160 | 400 | 0.00177 | 1.52 | 0.00138 | 1.59
320 | 800 | 0.00042 | 2.09 | 0.00030 | 2.18
Table 6: Numerical error $E^{N}$ in (47) for the initial condition (48) in example 5.5. $I$ | $\Delta t$ | $E,\alpha=0.5$ | EOC |  $E,\alpha$ in (23) | EOC
---|---|---|---|---|---
40 | 100 | 0.01692 | - | 0.01355 | -
80 | 200 | 0.00458 | 1.89 | 0.00351 | 1.95
160 | 400 | 0.00092 | 2.32 | 0.00067 | 2.38
320 | 800 | 0.00014 | 2.76 | 0.00001 | 2.80
Table 7: Numerical errors $E^{N}$ in (47) with the initial condition (49) in
example 5.5.

Next we evaluate the numerical artifact in the form of unphysical negative
oscillations in the case of initial condition (48). We plot the minimal value
for each numerical solution per each time step at all grid levels in Figure
12. We note that the extremal (rounded) values for each mesh with respect to
$n$ are $-0.0677$, $-0.0275$, $-0.0108$, and $-0.00163$ from the coarsest to
the finest mesh.

Finally, the plot of a difference between the initial (conserved) mass and the
actual one at each time step for all grid levels can be found in Figure 12.
Note that the initial mass (the integral of initial function) is approximated
by the value $0.031416$ for (48) and $n=0$ using

$M=M(I,n)=h^{2}\sum\limits_{i,j}\phi^{n}_{ij}\,.$ (50)

Figure 12: The minimum (left) and the difference to the initial mass (right)
for numerical solution of example 5.5.

##  6 Conclusions

We present the novel semi-implicit parametric family of one-dimensional
numerical schemes for conservative and non-conservative advection equation.
Using the Strang splitting for the advection in several dimensions, one can
obtain the numerical solutions of advection equations using an (in advance
known) fixed number of alternating substitutions. As the schemes are second
order accurate in time and space with unconditional von Neumann stability,
they can be considered as a good alternative to standard explicit and implicit
schemes for advection dominated problems.

## References

  * (1) Arbogast, T., Huang, C.S., Zhao, X., King, D.N.: A third order, implicit, finite volume, adaptive Runge–Kutta WENO scheme for advection–diffusion equations.  Computer Methods in Applied Mechanics and Engineering 368, 113–155 (2020)
  * (2) Carciopolo, L.D., Bonaventura, L., Scotti, A., Formaggia, L.: A conservative implicit multirate method for hyperbolic problems.  Computational Geosciences 23(4), 647–664 (2019)
  * (3) Frolkovič, P., Lampe, M., Wittum, G.: Numerical simulation of contaminant transport in groundwater using software tools of $r^{3}t$.  Computing and Visualization in Science 18(1), 17–29 (2016)
  * (4) Frolkovič, P., Mikula, K.: Semi-implicit second order schemes for numerical solution of level set advection equation on Cartesian grids.  Applied Mathematics and Computation 329, 129–142 (2018)
  * (5) Fuhrmann, J., Langmach, H.: Stability and existence of solutions of time-implicit finite volume schemes for viscous nonlinear conservation laws.  Applied Numerical Mathematics 37(1-2), 201–230 (2001)
  * (6) Hadjimichael, Y., Ketcheson, D.I., Lóczi, L.: Positivity preservation of implicit discretizations of the advection equation.  arXiv preprint arXiv:2105.07403 (2021)
  * (7) Hahn, J., Mikula, K., Frolkovič, P., Medl’a, M., Basara, B.: Iterative inflow-implicit outflow-explicit finite volume scheme for level-set equations on polyhedron meshes.  Computers & Mathematics with Applications 77(6), 1639–1654 (2019)
  * (8) Ibolya, G., Mikula, K.: Numerical solution of the 1d viscous Burgers’ and traffic flow equations by the inflow-implicit/outflow-explicit finite volume method.  In: Proceedings of ALGORITMY, pp. 191–200. STU Bratislava (2020)
  * (9) Izzo, G., Jackiewicz, Z.: Highly stable implicit–explicit Runge–Kutta methods.  Applied Numerical Mathematics 113, 71–92 (2017)
  * (10) Knodel, M.M., Kräutle, S., Knabner, P.: Global implicit solver for multiphase multicomponent flow in porous media with multiple gas phases and general reactions.  In: International Conference on Finite Volumes for Complex Applications, pp. 595–603. Springer (2020)
  * (11) LeVeque, R.J.: Finite volume methods for hyperbolic problems, vol. 31.  Cambridge university press (2002)
  * (12) Lozano, E., Aslam, T.D.: Implicit fast sweeping method for hyperbolic systems of conservation laws.  Journal of Computational Physics 430, 110039 (2021)
  * (13) May, S., Berger, M.: An explicit implicit scheme for cut cells in embedded boundary meshes.  Journal of Scientific Computing 71(3), 919–943 (2017)
  * (14) Mikula, K., Ohlberger, M., Urbán, J.: Inflow-implicit/outflow-explicit finite volume methods for solving advection equations.  Applied Numerical Mathematics 85, 16–37 (2014)
  * (15) Pártl, O., Beneš, M., Frolkovič, P., Illangasekare, T., Smits, K.: Numerical modeling of non-isothermal gas flow and NAPL vapor transport in soil.  Computer Physics Communications 202, 175–187 (2016)
  * (16) Paszke, A., et. al.: Pytorch: An imperative style, high-performance deep learning library.  In: H. Wallach, et. al. (eds.) Advances in Neural Information Processing Systems 32, pp. 8024–8035. Curran Associates, Inc. (2019)
  * (17) Polívka, O., Mikyška, J.: Compositional modeling in porous media using constant volume flash and flux computation without the need for phase identification.  Journal of Computational Physics 272, 149–169 (2014)
  * (18) Puppo, G., Semplice, M., Visconti, G.: Quinpi: integrating conservation laws with CWENO implicit methods.  arXiv preprint arXiv:2102.00741 (2021)
  * (19) Qin, T., Shu, C.W.: Implicit positivity-preserving high-order discontinuous Galerkin methods for conservation laws.  SIAM Journal on Scientific Computing 40(1), A81–A107 (2018)
  * (20) Radu, F.A., Pop, I.S., Attinger, S.: Analysis of an Euler implicit-mixed finite element scheme for reactive solute transport in porous media.  Numerical Methods for Partial Differential Equations: An International Journal 26(2), 320–344 (2010)
  * (21) Uçar, Y., Yağmurlu, N.M., Çelikkaya, İ.: Operator splitting for numerical solution of the modified Burgers’ equation using finite element method.  Numerical Methods for Partial Differential Equations 35(2), 478–492 (2019)
  * (22) Zhang, L., Ge, Y.: Numerical solution of nonlinear advection diffusion reaction equation using high-order compact difference method.  Applied Numerical Mathematics 166, 127–145 (2021)
  * (23) Zhao, H.: A fast sweeping method for eikonal equations.  Mathematics of computation 74(250), 603–627 (2005)


alarm there is something wrong in the protocol such that we can perform
analysis on the leader or even increase network waiting time.

If we view LSC from a broader perspective that takes clients into
considerations, the next section shows that the LSC parameters $\eta$ has a
great impact on the amount of censorship resistance to client’s experience and
system efficiency. In general, the less is the $\eta$ value, the better is the
client’s experience and system efficiency.

##  19 Navigator protocol

###  19.1 Properties

A replica is available for a DA-CR instance if it can insert some transactions
into its mini-block before sending it to a leader. Given $f$ malicious
replicas in the system, the Navigator protocol satisfies the following
properties:

  * •

Strict Transaction Inclusion in DA layer: Given a DA-CR instance with $t$
guaranteed honest mini-block inclusion, if an honest client sends $n-t+1$
copies of a transaction to distinct available replicas, then the instance is
either Incomplete, or Complete with the client’s transaction.

  * •

Probabilistic Transaction Inclusion in DA layer Under Honest Leader: Given a
DA-CR instance with $t$ guaranteed honest mini-blocks inclusion under an
honest leader, if an honest client sends $x$ copies of a transaction to random
distinct available replicas and the instance Completes, the probability of
inclusion monotonically increases with $x$ and $t$. If network latency of
honest replicas are i.i.d, the probability increases exponentially with $x$
and $t$.

  * •

Probabilistic Transaction Inclusion in DA layer Under Malicious Leader: Given
a DA-CR instance with $t$ guaranteed honest mini-blocks inclusion under a
malicious leader, if an honest client sends $x$ copies of one transaction to
random distinct available replicas and the instance Completes, the probability
of inclusion monotonically increases with $x$, up to $1$ when $x=n-t+1$. When
$x\leq f$ the malicious leader can censor the transaction with probability
$1$, but to censor $u$ distinct transactions each of $x$ copies altogether,
the probability of censoring decreases exponentially as $u$ increases.

  * •

2 $\Delta$ Probabilistic Confirmation in DA layer: Given a DA-CR guarantees at
least $t$ mini-blocks inclusion from honest replicas, if an honest client
sends $x$ copies of transaction for inclusion, the client can compute the
probability of inclusion right after 2 $\Delta$ network latency based on
number of responses from $x$ replicas and belief on the leader’s type.

Those properties are proved in Appendix 19.8.

###  19.2 Algorithm

The Algorithm is presented in Algorithm 5. Appendix 19.3 provides a concrete
implementation for the procedure call at line 4, which allows client to get a
list of available replicas.

Algorithm 5 Navigator

1:At the client:

2:procedure $submit$( $x$, $tx$)

3: $v\leftarrow currentBlockNum$

4: $S\leftarrow getAvailableReplicas(x,v)$, $cn\leftarrow 0$

5: if $|S|<x$ then

6: return $FAIL,|S|$

7: for $j\text{ in }S$ do

8: $r\leftarrow\text{send }\langle tx,v\rangle\text{ to replica }j$

9: $cn\leftarrow cn+r$

10: return $cn$

11:At the replica  $p$ :

12:upon $\textit{receiving }\langle tx,v\rangle$ do

13: if $(v>\textit{myView})\lor(v=\textit{myView }\wedge\text{waiting for
Signal})$ then

14: $\text{insert }\langle tx,v\rangle\text{ to mempool}$ and respond 1

15: else

16: respond 0

###  19.3 Client and Replicas time synchronization

We present a time synchronization protocol in Algorithm 6. Essentially, the
algorithm requires the client to periodically pings each replicas about their
status about when they would submit the next mini-block. Assuming the worst
cast network delay $\Delta$ after GST, the client can assure they can make the
deadline. If the network has not reached GST, the BFT protocol might not reach
agreement. The accuracy of the status synchronization depends on the ping
frequency, which is a tunable parameter for each client.

Algorithm 6 Navigator Synchronization

1:$S_{v}\leftarrow\emptyset$, $v\leftarrow\textit{0}$ ${\cdots INF}$, $timer$

2:At the client:

3:procedure $cronJob$

4: for $i\leftarrow 1\cdots n$ do

5: $s\leftarrow send\langle
PING,v,\textit{timer.now}\rangle\quad\quad\text{\\# timeout 2}\Delta$

6: $S_{v}\leftarrow S_{v}\cup s$

7: $\textit{\\# }\Delta\textit{ is latency in network model}$

8: $S_{v}\leftarrow\\{s:s\in S_{v},\textit{timer.now}<s.deadline\\}$

9:At the replica  $p$ :

10:upon $\textit{receiving }\langle PING,v,tn\rangle$ do

11: if $(v>\textit{myView})\lor(v=\textit{myView }\wedge\textit{currAt =
SIGNAL})$ then

12: $\textit{deadline}\leftarrow tn+2\Delta+\textit{timeToNextSignal}$

13: $send\langle PONG,deadline\rangle$

14: else

15: $send\langle PONG,0\rangle$

16:procedure getAvailableReplicas( $x$, $tv$)

17: return $\textit{shuffle}(S_{tv})[x:]$ #  return only x replicas

Figure 6: The probability of inclusion given the leader is malicious with
varying $f$. The x-axis is a ratio of the number of copies of transactions and
$2f$

###  19.4 Empirical Evaluation on Inclusion Probability when Leader is honest

Figure 7 plots the probability of censorship for a transaction using a DA-CR
with $t=f+1$ and $n=3f+1$ when the leader is honest, such that $t\leq q\leq
n-f$, where $f=10$ and $n=3f+1=31$.

Figure 7: Under honest leader, the probability of a transaction being censored
as a function of $(x,q)$, where $1\leq x\leq n$, $f<q\leq 2f+1$

###  19.5 Analysis on probability of inclusion given malicious leader with
$f+1\geq x$

Figure 6 plots the probability of inclusion with varying number of replicas
$n=3f+1$. The inclusion can be satisfied with high probability even the client
sends less than $2f+1$ copies, approximately $1.8f$. The probability of
inclusion decreases as the clients send less number of copies of transactions.

###  19.6 Upper Bound on Total Censored transactions

In Appendix 19.7, we prove with Hoeffding’s inequality and union bound that
the probability of removing $\alpha$ fraction out of total $T$ transactions is
upper bounded by $2{n\choose n-f}\exp{\\{-2(p+\alpha-1)^{2}T\\}}$, where
$p=1-\prod_{i=0}^{x-1}\frac{n-f}{n-i}$. As $T$ increases, the probability
decreases exponentially.

###  19.7 Analysis on the upper bound on the total number of censored
transaction

Suppose there is a total $T$ unique transactions, denoted as $t_{1},..t_{T}$
by all clients, where each transaction has $x$ copies which are randomly sent
to replicas without replacement. For a small constant ratio $\alpha$, it is
exponentially unlikely that the leader would be able to censor more than
$\alpha T$ transaction.

Define a set of mini-blocks $S$ of cardinality $n-f$, we want to compute the
probability that a malicious leader can find a particular $S$ such that the
union of all transactions from mini-blocks in set $S$ has a size less than
$(1-\alpha)T$. For the same reason we discussed in Section 5.2.2, the
probability of inclusion is entirely depend on how likely a client sends its
transaction to honest replicas; since adversary can always remove those
transaction by selectively accepting other honest mini-block. The probability
of such set $S$ exists can be written as Equation 7, where every transaction
corresponds to an independent random indicator variable,
$\mathbbm{1}_{t_{i}\in S}$. The variable is independent because every
transaction runs its own instance of Algorithm 5. It is a Bernoulli random
variables that evaluates to 1 if condition $t_{i}\in S$ holds true; otherwise
0. The probability of the condition being true is
$p=1-\prod_{i=0}^{x-1}\frac{n-f}{n-i}$.

$\displaystyle\Pr(\textit{remove }\alpha T)$
$\displaystyle=\Pr_{{|S|=n-f}}\Bigl{\\{}\exists
S:\sum_{i=1}^{T}\mathbbm{1}_{t_{i}\in S}<(1-\alpha)T\Bigr{\\}}$ (7)
$\displaystyle\leq{n\choose
n-f}\Pr\Bigl{\\{}S:\sum_{i=1}^{T}\mathbbm{1}_{t_{i}\in
S}<(1-\alpha)T\Bigr{\\}}$ (8)

The inequality from Equation 7 to 8 comes from Boole’s inequality. To get the
probability on the right high side, we use a random variable
$\frac{\sum_{i=1}^{T}\mathbbm{1}_{t_{i}\in S}}{T}=X$, whose expected value is
$X$ is $p$, by rearranging terms of Equation 8, we get

$\displaystyle\Pr\Bigl{\\{}S:\sum_{i=1}^{T}$
$\displaystyle\mathbbm{1}_{t_{i}\in
S}<(1-\alpha)T\Bigr{\\}}=\Pr\Bigl{\\{}X<1-\alpha\Bigr{\\}}$ (9)
$\displaystyle=\Pr\Bigl{\\{}E[X]-X>E[X]+\alpha-1\Bigr{\\}}$ (10)
$\displaystyle\leq\Pr\Bigl{\\{}\left|E[X]-X\right|>E[X]+\alpha-1\Bigr{\\}}$
(11) $\displaystyle\leq 2\exp{\\{-2(p+\alpha-1)^{2}T\\}}$ (12)

In Equation 11, we let $\alpha>1-E[X]=\prod_{i=0}^{x-1}\frac{n-f}{n-i}$,
inequality is due to
$\Pr\Bigl{\\{}\left|E[X]-X\right|>E[X]+\alpha-1\Bigr{\\}}=\Pr\Bigl{\\{}E[X]-X>E[X]+\alpha-1\Bigr{\\}}+\Pr\Bigl{\\{}X-E[X]<-E[X]-\alpha+1\Bigr{\\}}$.
Inequality 12 comes from Hoeffding’s Inequality.

The probability of removing $\alpha T$ transactions is therefore upper bounded
by ${2{n\choose n-f}\exp{\\{-2(p+\alpha-1)^{2}T\\}}}$.

As it can be shown that given $300$ replicas and $200$ transactions per mini-
block, a malicious leader can censor more than $5\%$ of the transactions.

###  19.8 Navigator Properties Proof

In this section we prove the four properties stated in Section 5.

Strict Transaction Inclusion in DA layer: Suppose there is a DA-CR instance
with $t$ guaranteed honest mini-block inclusion, and the client sends $n-t+1$
copies of the same transaction to distinct replicas which are available for
the instance round, then suppose the leader is malicious and only includes the
minimal $t$ mini-blocks from honest replicas, which is required by the DA-CR.
By pigeonhole principle, there must be at least one transaction received by
one of the $t$ honest replicas whose mini-block would get included. Similar
argument can be made for honest leader.

Figure 8: The probability of that , $X=i$, $i$ number of honest replicas
receiving the transaction given $x$ copies was sent. We show one conditional
probability when $q=16$ assuming network latency is i.i.d to illustrate the
point that both $q,x$ needs to work together to increase the probability of
$\Pr{(IN)}$.

Probabilistic Transaction Inclusion in DA layer Under Honest leader: Suppose a
DA-CR instance with $t$ guaranteed honest mini-blocks inclusion, if an honest
client sends $x$ copies of a transaction to random distinct available
replicas. We first show that the probability of inclusion increases with $x$
and $t$ in the feasible region where a leader cannot exclude all $x$ copies
without making the instance Incomplete. From Section 5 and Table 3, we can
define the feasible region as $(0<x<n,t\leq q\leq n-f)$ when the leader is
honest; By our first properties, we already knows when $x\geq n-t+1$, the
probability of inclusion is $1$. We now focus on the rest of the case.

When the leader is honest and $(0<x<n-t+1,t\leq q<n-f)$, Section 5.2.1 shows
that

$\displaystyle\Pr(\textit{IN; x,q})$
$\displaystyle=1-\sum_{i=1}^{x}\Pr(\textit{None }|X=i;q)\Pr{(X=i;x)}$
$\displaystyle=1-\sum_{i=1}^{x}\Pr(\textit{None }|X=i;q)\frac{{n-f\choose
i}{f\choose x-i}}{{n\choose x}}$ (13)

where $X$ is a random variable indicating the number of honest replicas that
has client’s transaction. $P(q,i)=1-\Pr(\textit{None }|X=i;q)$ is the
probability none of $X$ honest mini-block containing the client’s transaction.
We note that $P(q,i)$ is monotonically increasing in $q$, it is because as
there are more guaranteed honest mini-block inclusion, the chance that one of
the $i$ mini-blocks containing the transaction gets included is non-
decreasing. Since $q\geq t$, $P(q,i)$ increases as $t$ increases. Because
$\Pr{(X=i;x)}$ does not depend on $t$, the overall probability $\Pr{(IN)}$ as
a function of $q$ increases as $P(q,i)$ increases. Therefore $\Pr{(IN)}$ is
monotonically increasing as $q$ increases.

We now argue $\Pr{(IN)}$ is monotonically increasing as $x$ increases. We
first show that $P(q,i)$ is monotonically increasing in $i$ while fixing $q$.
Since the honest leader is collecting $q$ honest mini-blocks out of $n-f$
honest replicas. If the number of honest replicas whose mini-block containing
the transaction increases, then the chance of $P(q,i)$ is at least non-
decreasing. Now in the Equation 13, there is a conditional sum over all
possible value which $X$ can take up to $x$; in the worst case when
$\Pr(\textit{None }|X=i;q)$ is a constant ratio, then the final probability
$\Pr{(IN)}$ is also be a constant. But as the conditional probability starts
to decrease, the overall $\Pr{(IN)}$ increases. In addition to that, as $x$
increases, the distribution of $P(X=i;x)$ gets shifted rightward, such that it
is more common (higher probability) that more honest replicas receiving the
transaction. The intuition behind the phenomena is as $x$ increases, the
numerator terms of $\frac{{n-f\choose i}{f\choose x-i}}{{n\choose x}}$ finds
its peak also at higher $i$. Figure 8 plots distribution for different $x$
when $f=10$. Hence if the conditional probability drop fast enough, a client
only needs to send a few copies $x$, such that the permutation term (binomial
like) raises up just after the conditional term drops, the probability can be
kept very low. In Figure 8, the conditional probability is using i.i.d
assumption with number of included mini-block $q=16$, out of $21$ possible
honest mini-blocks with $n=31$.

Now, we show the second statement that if the network latency is i.i.d, the
inclusion probability increases exponentially with $x$ and $t$. Since network
latency is i.i.d, the conditional part can be computed as

$\displaystyle\centering\Pr(\textit{None }|X=i;q)\@add@centering$
$\displaystyle=\prod_{j=0}^{q-1}\frac{n-f-i-j}{n-f-j}$
$\displaystyle=\prod_{j=0}^{q-1}\frac{a-i}{a}$
$\displaystyle<(1-\frac{i}{n-f})^{q}$ $\displaystyle\approx
e^{-q\frac{i}{n-f}}$

Hence we know the term $\Pr(\textit{None }|X=i;q)$ is exponentially small when
$q$ increases. However, we need also show that the summation term in Equation
13 is small.

First note that, the permutation term, $\Pr{(X=i;x)}$, is always less than
$1$, so we can just focus on the case when $i$ is small, because when $i$ is
large, the conditional probability part is already very small.

Now we use the binomial upper bound, ${n\choose k}<(\frac{ne}{k})^{k}$[17],
for numerator, and binomial lower bound, $(\frac{n}{k})^{k}<{n\choose k}$, for
denominator to show that in worse case, when $i$ is small the probability is
exponential decreasing when $x$ increases

$\displaystyle\centering\frac{{n-f\choose i}{f\choose x-i}}{{n\choose
x}}\@add@centering$
$\displaystyle<\frac{(\frac{(n-f)e}{i})^{i}(\frac{fe}{x-i})^{x-i}}{(\frac{n}{x})^{x}}$
$\displaystyle=\frac{(n-f)^{i}f^{x-i}x^{x}e^{x}}{i^{i}(x-i)^{x-i}n^{x}}$
$\displaystyle=\frac{(n-f)^{i}f^{x}x^{x}e^{x}}{i^{i}x^{x}n^{x}}$
$\displaystyle\approx\frac{(n-f)^{i}}{i^{i}}(\frac{e}{3})^{x}$

Clearly for small $i$, as $x$ increases the value is upper bound by a
exponentially decreasing value. Hence it completes the proof.

Probabilistic Transaction Inclusion in DA layer Under Malicious Leader: The
first part of the statement that the probability increases as $x$ increases
comes directly from the fact from Section 5.2.2 where $f+1\leq x<n-t+1$ that

$\Pr(IN)=\sum_{i=f+1}^{x}\frac{{n-f\choose i}{f\choose x-i}}{{n\choose x}}$

Clearly, $\Pr(IN)$ increases as $x$ increases. If $x\leq f$, the probability
is $0$, because the malicious leader can choose the rest $f+1$ honest replicas
which does not have the transaction.

To show that a malicious leader cannot censor all $u$ unique transactions
altogether, we use the result from Section 19.6, where
$p=1-\prod_{i=0}^{x-1}\frac{n-f}{n-i}$

$\Pr(\textit{Censor all u txs})=2{n\choose
n-f}\exp{\\{-2(p+\alpha-1)^{2}T\\}}$

As we can see, given a fixed pair of $n,x$, the probability that the censoring
attack succeeds is decreasing exponentially in either $T$, $\alpha$ or $x$. By
letting $u=\alpha T$, we get our property.

2 $\Delta$ Probabilistic Confirmation in DA layer: By Algorithm 5, the clients
is able to get the response if some replicas are able to submit transactions.
The client can calculate its probability of inclusion based on the response.
Although some responses might be sent from malicious replicas, the calculation
stated in the previous properties can apply directly to derive the probability
result, since they already assume the malicious replicas are not including the
transaction. Since each categories have the explicit formula to compute the
probability, the client can use the received the responses, the number of
duplication and current belief on the leader to compute the probability of
inclusion.

##  20 Spamming Attack and Transaction De-Duplication

BigDipper expects clients to duplicate its transactions among replicas to
increase censorship resistance. However, without proper mechanism clients are
inclined to spam the BFT replicas for increasing the chance of inclusion. A
direct way of controlling the spamming problem is to increase the cost of
transaction. We could use a price mechanism, like EIP1559 [37] deployed in
Ethereum, which uses posted price when the network is not congested, and
otherwise use first price auction. But the price mechanism alone does not the
problem, since the price is raised for everyone. The final cost for a client
should also correspond to the number of copies a transaction has been sent.
This can be enforced by requiring the execution layer (virtual machine) to
also charge fees by transactions counts as a meta information when calculating
the spent resources.

Another method that does not rely on in-protocol mechanism like 1559 is to use
market mechanism entirely, such that every BFT replicas is free to post its
own prices. The $n$ replicas creates a competitive market. If we consider
there are 2 types of transaction hold by various clients, and let the BFT
replicas being the first mover in the information asymmetry game, a natural
solution is to use screening mechanism. We can find a separating subgame
perfect Nash equilibrium if the number of high value transaction type is
relatively small compared to the low value type transaction. The screening
mechanism has been been well studied in the field of contract theory [36].
However in the screening model, the BFT replicas would lose some profits from
high value transaction, because the high value type would mimic the type of
low value transactions. Then the replicas can use the value to sort the
transaction and only include the ones that fit into the mini-blocks.

For spamming concern at network level, we can rely on standard approaches
including content distribution network, proxy networks which have been already
used in the real world scenario.

###  20.1 De-duplication

As client sends multiple copies of the same transaction, it might be worrisome
that use of bandwidth is inefficient. However, as how to handle the spamming
problem, the usage of block space is essentially economical. If a client
desires more properties that take more block space, the client should pay for
more. It is beneficial if the system can remove the inefficiency entirely, but
given that BigDipper is hyperscalable, having some copies of transaction is
considered acceptable. For most time-insensitive transaction, censorship
resistance is not a strict requirement.

Figure 9: Transaction reordering from 5 mini-blocks, each contains varying
number of transactions. The commitment is denoted by first two letter of their
color ($c_{re}=0$). $c_{ord}$ is the commitment of ordered transactions.

##  21 Zero knowledge permutation argument

###  21.1 Permutation Argument

We define the property that an ordering function has to satisfy in order to
ensure censorship resistance of the protocol. We focus on an ordering function
that can be not deterministically derived from the data, and requires
exogenous input and leader’s coordination. Those include partial block
auction, whole block auction, VRF.

  * •

Correct Permutation: If an honest replicas has decided on some ordered
transactions for some view in the BFT protocol, then the decided order must be
a correct permutation of the dispersed data without corrupting integrity of
any transactions.

Suppose $B_{ord}$ is the resultant ordering after the exogenous ordering
function. In the simple case, if the entire ordered block $B_{ord}$ is
broadcast by the leader and each replica has retrieved all data by using the
Retrieve invocation, each replica can enforce the correct permutation property
by checking if the ordered dispersed block is a subset of the ordered block.
We now discuss the case when no replica retrieves the entire data. It is
achieved with permutation argument. A visualization of a permutation argument
is shown in Figure 9. The leader needs to sends three items for replicas to
construct the permutation circuit used for verifying validity, which include a
permutation annotation(function), annotations of the transactions size and
hash inside each mini-block and a polynomial commitment after the ordering.
The annotations of transactions $A_{i}$ for mini-block $i$ can be provided and
encoded by replicas when creating the mini-block. To prevent the leader from
corrupting the annotation of transactions, the leader has to provide witnesses
against the commitment list for all $A_{i}$ to every replicas.

Then the leader publishes the commitment of the ordered transactions by
sending the actual ordered block $B_{ord}$, or its commitment $c_{ord}$. A
permutation circuit with Halo2 [43] Plonkish arithmetizaton can be built to
ensure a correct permutation with inputs from $A_{i}$ and $A_{ord}$, such that
all transaction boundary are preserved.

##  22 Scooping Protocol

###  22.1 Scooping Algorithm

The scooping algorithm is shown in Algorithm 7.

Algorithm 7 Scooping

1:At the replicas $p$:

2:procedure verifyCard( $commit,v$)

3: return ms.verify(aggSigs[ $commit,v$])

Figure 10: A flow diagram depicting an adaptation of BigDipper with Hotstuff.
The procedures from Hotstuff-2 protocol is marked in black; Card-$7$ protocol
is marked with red, and Scooping protocol is marked in blue. The signature
validation procedure has to be verified after receiving the settlement
message.

##  23 Integration with Hotstuff-2

Figure 10 depicts a flow diagram that integrates Card-$7$ with Hotstuff-2. We
encourage readers to read about Hotstuff-2[30] protocol. $C_{v}(\cdot)$
represents a certificate at view $v$, $C_{v}(C_{v}{\cdot})$ represents a
double certificates, which indicates the block at view $v$ can be safely
committed by anyone. $k$ denotes block height, since no all view can
successfully commit a block, hence $k$ is different from $v$.

On the diagram, the blue bar represents the Hotstuff-2 leader replica, the
green bar at the bottom represents a non-leader replicas. Procedures needs to
be carried out after receiving the messages are written close to the message
arrow, and the sequence of execution is from top to down.

Because we are using verifiable information dispersal, all block $B$ from the
Hotstuff-2 protocol actually means the metadata overhead described in the
Card-$7$ protocol, which are commitment list and attestation set.

All replicas need to three algorithms: Hotstuff-2 (in black and purple),
Card-$7$ (in red) and Scooping (in blue). The time flows from left to right,
transitioning from view $v-1$ to $v$. Starting from Collection message from
the left, it belongs to the Card-$7$ protocol when a replica tries to get its
mini-block included by the leader; the vote2 message belongs to Hotstuff-2
protocol, the leader needs to collect at least $2f+1$ vote2 messages before
proceeding; vote2 message is sent at a stage that the replica has locked some
block, and is planning to vote the second time to commit the block. The
subscript indicates that the block $B_{k-1}$ at blockchain height $k-1$
happens at view number $v-1$.

When the leader collects at least $2f+1$ Collections messages and $2f+1$ vote2
messages. The leader performs the 2D RS encoding and sends the Dispersal
message of the Card-$7$ protocol, at the same time the leader also sends the
Propose message along with the double certificates to the replica. This is the
official start of the new block, the block is tagged with view number $v$; in
the diagram we colored all procedures relates to this block in purple. At the
same time, the previous block $B_{k-1}$ has just been able to commit after
creating a double certificate by collecting signatures from $2f+1$ vote2
messages, the double certificate are piggybacked to the replicas.

After receiving those messages at the replica side, the replica would directly
commit the block that has double certificates, in addition it would indirectly
commit all blocks that are prefix of the block just being committed. In the
diagram, we represents them as $B_{k^{\prime\prime}}$. The non-leader replica
then needs to perform the dispersal check which corresponds to the kzg
commitment checking and signature checking from Card-$7$ Protocol. After
everything passes, the replica sends the Approve message of Card-$7$ and vote
message of Hotstuff-2 to the leader.

The leader waits to receive messages from $2f+1$ replicas and creates the
combined signature for Agreement message of Card-$7$ and similarly a
certificate $C_{v}(B_{k})$ for prepare message of Hotstuff-2. Optionally the
leader can assign order of transaction in this step; however, for us we use
deterministic ordering so the ordering message should be empty.

At the replica side, after receiving the first certificate for $B_{k}$, the
first step is to check if the validate procedure from the Scooping protocol is
successful. If it is not, then the data is not available, so the replica
should not lock the certificate. Otherwise, the replica proceeds with regular
Hotstuff-2 procedures. At the same time, the leader is preparing for a new
view at $v+1$, so the replica should prepare its new mini-block $b_{k+1}$ and
sends to the leader. And we just complete a cycle of optimistic path when
nothing wrong happens.

##  24 Properties proofs to BigDipper-$7$ Hotstuff-2 integration

###  24.1 Properties proofs to Safety and Liveness

We prove that the integrated protocol satisfies safety and liveness properties
of BFT protocol, which has great similarity to the proof in [30]. All
properties listed in Section 10.1 can be proved by applying properties proven
in Card-$7$ and Navigator protocols.

Safety: We prove by contradiction, suppose two blocks $B_{l},B_{l}^{\prime}$
commit at the same height $l$. Since DA-CR has Correctness, Commitment binding
and Data-tampering Resistance properties, both block must be made available
and properly certified in different BFT instances, otherwise it cannot pass
Alg 1 line 7. Suppose both are committed as result of direct commits of
$B_{k},B_{k}^{\prime}$ at view $v,v^{\prime}$ respectively. Suppose $v\leq
v^{\prime}$ (if $v=v^{\prime}$, then suppose $k\leq k^{\prime}$), by Lemma 4
from Appendix 24.3, certificate $C_{v^{\prime}}(B_{k^{\prime}}^{\prime})$ must
extend $B_{k}$, hence $B_{l}=B_{l}^{\prime}$.

Liveness: The integrated protocol uses identical pacemaker as Hotstuff-2.
Since DA-CR has Termination property and Hotstuff-2 is live by itself, the
protocol continues to next view if a malicious leader stalls the protocol or
the network is asynchronous. When the leader is honest and the network is
synchronous, in the optimistic path where the double certificate is received,
the leader can collect at least $n-f$ mini-blocks (line 11) from replicas. If
the honest leader is not in the optimistic path, the leader must wait for
$3\Delta$ before proposing, Line 3 requires the replica to send the collection
message along with the locked certificate. Consequently the leader always have
data to propose in both cases. The Termination property of DA-CR ensures
dispersal can finish with type Complete, so BFT can make progress.

###  24.2 Properties proofs to Addditional Censorship resistance properties

Strict Transaction Inclusion: Because Card-$7$ provides Inclusion-(t) with
$t=f+1$, the strict inclusion directly come from Strict transaction in DA
layer property provided by the Navigator algorithm. Since the integration with
Hotstuff-2 is safe and live, under an honest client, a client can have
guarantee of inclusion with probability $1$ by sending $n-t+1$ copies of
transaction. If a malicious leader excludes the transaction that has been sent
to $n-t+1$ replicas, also due to Strict transaction in DA layer property
provided by the Navigator algorithm, the DA-CR results in an Incomplete state,
which no valid combined signature can be generated, hence the integrated
BigDipper-$7$ Hotstuff-2 integration would not proceed because no honest
replica can pass the Scooping Algorithm, which is inserted into the consensus
path in Alg.Scooping line 11.

Probabilistic Transaction Inclusion Under Honest Leader: The proof is almost
identical as above. Since the integration with Hotstuff-2 is safe and live, a
client can send $x$ copies of transaction to probabilistically reach
sufficient replicas, the rest follow through the Probabilistic Transaction
Inclusion in DA layer Under Honest Leader from Navigator. Since the integrated
protocol is live and safe, the probability of client’s transaction is included
increases exponentially with $x$ and $t$.

Probabilistic Transaction Inclusion Under Malicious Leader: The proof is
almost identical as above. If the leader is malicious and the instance
Completes, we use can use Probabilistic Transaction Inclusion in DA layer
Under Malicious Leader from the Navigator protocol, to derive that the
inclusion probability is monotonically increasing with $x$, and when $x\leq
f$, the probability of censoring $u$ transactions altogether decreasing
exponentially as $u$ increases.

$2\Delta$ Probabilistic Confirmation for Clients: It is a direct consequence
of Card property and Navigator property from Section 5.

Hyperscale Throughput: Since all BigDipper protocols use verifiable
information dispersal, the throughput of reaching consensus to a sequence of
bits in not bottleneck from downloading the whole data of size $O(bn)$, . The
system capacity is summation of bandwidth from all replicas, which is
$O(C_{band}n)$, where $C_{band}$is bandwidth for each replica. Then the limit
on rate of information dispersal for each replicas needs to be $O(b)\leq
C_{band}$, so at maximum the rate of reaching consensus on a sequence of bits
is only a constant of system capacity.

###  24.3 Lemma in support for safety proofs of BigDipper Hotstuff-2
integration

We show that lemma4.1 from Hotstuff-2[30] still holds in the integrated
protocol.

###### Lemma 4

If a party directly commits block $B_{k}$ in view $v$, then a certified block
that ranks no lower than $C_{v}(B_{k})$ must equal or extend $B_{k}$.

###### Proof

We use the same approach as the Hotstuff-2 proof, we are mostly concerned that
the data is unavailable, we present the whole proof and shows the data has to
be available for any committed blocks. Suppose we consider a block
$B_{k^{\prime}}^{\prime}$ that is certified in view $v^{\prime}$ to produce
certificate $C_{v^{\prime}}(B_{k^{\prime}}^{\prime})$ has a rank higher than
$C_{v}(B_{k})$ if $v^{\prime}>v$. We will use induction on view $v^{\prime}$.

At the base case $v=v^{\prime}$, $C_{v^{\prime}}(B_{k^{\prime}}^{\prime})$ is
certified in view $v$. But it is impossible, since no honest replica would
vote twice per a single view.

At the induction step, since $C_{v}(C_{v}(B_{k}))$ exists, there must be a set
of replicas $S$ whose size is no less than $2f+1$ that possess $C_{v}(B_{k})$.
Since a double certificate can only be created with at least $f+1$ honest
replicas, we know the data $B_{k}$ is available. All replicas in $S$ must have
locked on $B_{k}$ or a block that extends $B_{k}$ at the end of $v$. By
induction assumption, any certified block that ranks equally or higher than
$v$ must equal or extend $B_{k}$. By the end of $v^{\prime}$, the replicas in
$S$ are still locked at $B_{k}$ or blocks that extends $B_{k}$. The leader has
two ways to enter a new view. In the first case when the leader has a double
certificate, then at least $f+1$ honest replicas have ensured the data is
available, so it is safe to extend on top of it. In the second case, the
leader finds out the highest certified block by waiting $3\Delta$ to receive
all honest replicas. Since the leader verifies that all locked certificate
must have a valid aggregate signature corresponding to it, this step is
ensured by Alg 1 line 5. So we know it is safe for the leader to build on top
of it. Resume the Hotstuff-2 proof. Consider a proposal at view
$v^{\prime}+1$, if the leader makes a proposal $B_{k^{\prime}}^{\prime}$ that
does not extend $B_{k}$, its certificate must be ranked lower than
$C_{v}(B_{k})$. As the result, no honest replica would vote for it, so the
certificate $C_{v^{\prime}+1}(B_{k^{\prime}}^{\prime})$ cannot be formed.
Hence if $B_{k^{\prime}}^{\prime}$ was committed,
$C_{v^{\prime}+1}(C_{v^{\prime}+1}(B_{k^{\prime}}^{\prime}))$ is formed,
$B_{k^{\prime}}^{\prime}$ must extend $B_{k}$.

To clarify special cases, when a malicious leader correctly runs the Card
protocol but refuses to broadcast the aggregate signature to anyone, the
malicious leader eventually gets replaced. However, since no one has the
aggregate signature, a new block proposed by an honest leader would not build
on top of that block after waiting $3\Delta$. If a malicious leader broadcasts
the aggregate signature, but refuses to create the double certificates, any
honest replica can send the locked certificate along with the aggregate
signature to convince the leader to build on top of it.