# The energy method for high-order invariants in shallow water wave equations

Qifeng Zhang  Tong Yan  Guang-hua Gao  Department of Mathematics, Zhejiang
Sci-Tech University, Hangzhou, 310018, China  Department of Mathematics,
Nanjing University of Posts and Telecommunications, Nanjing, 210096, China

###### Abstract

Third order dispersive evolution equations are widely adopted to model one-
dimensional long waves and have extensive applications in fluid mechanics,
plasma physics and nonlinear optics. Among them are the KdV equation, the
Camassa–Holm equation and the Degasperis–Procesi equation. They share many
common features such as complete integrability, Lax pairs and bi-Hamiltonian
structure. In this paper we revisit high-order invariants for these three
types of shallow water wave equations by the energy method in combination of a
skew-adjoint operator $(1-\partial_{xx})^{-1}$. Several applications to seek
high-order invariants of the Benjamin-Bona-Mahony equation, the regularized
long wave equation and the Rosenau equation are also presented.

###### keywords:

Energy method; High-order invariant; Shallow water wave equation

##  1 Introduction

A family of third order dispersive evolution equations of the form

$\displaystyle u_{t}-\alpha^{2}u_{xxt}+\gamma
u_{xxx}+c_{0}u_{x}=(c_{1}u^{2}+c_{2}u_{x}^{2}+c_{3}uu_{xx})_{x},\quad x\in
R,\;t>0$ (1.1)

frequently appeared in the simulation of the shallow water waves, see e.g.,
[1], where $\alpha$, $\gamma$ and $c_{i}$ $(i=0,1,2,3)$ are real constants;
$u$ denotes a horizontal velocity field with the independent spatial variable
$x$ and temporal variable $t$.

A typical such equation (1.1) with $\alpha^{2}=c_{0}=c_{2}=c_{3}=0$,
$c_{1}=2$, $\gamma=-2$ is _the KdV equation_

$\displaystyle u_{t}-4uu_{x}-2u_{xxx}=0,\quad x\in R,\;t>0,$ (1.2)

which describes the unidirectional propagation of waves at the free surface of
shallow water under the influence of gravity. The first four invariants of
(1.2) are respectively as (see e.g., [2], although there is a minor typo in
the coefficient of the fourth invariant, it does not affect the reading of
this classic review)

$\displaystyle M_{1}=\int_{R}u\mathrm{d}x,\quad
M_{2}=\int_{R}u^{2}\mathrm{d}x,\quad
M_{3}=\int_{R}\Big{(}u_{x}^{2}-\frac{2}{3}u^{3}\Big{)}\mathrm{d}x,\quad
M_{4}=\int_{R}\Big{(}u_{xx}^{2}-\frac{10}{3}uu_{x}^{2}+\frac{5}{9}u^{4}\Big{)}\mathrm{d}x.$

Taking $\alpha^{2}=c_{3}=1$, $\gamma=c_{0}=0$, $c_{1}=-\frac{3}{2}$,
$c_{2}=\frac{1}{2}$, we have another example called _the Camassa–Holm
equation_ [3]

$\displaystyle u_{t}-u_{xxt}+3uu_{x}=2u_{x}u_{xx}+uu_{xxx},\quad x\in
R,\;t>0,$ (1.3)

which models the unidirectional propagation of shallow water waves over a flat
bottom. The first three invariants are listed as follows

$\displaystyle E_{1}=\int_{R}(u-u_{xx})\mathrm{d}x,\quad
E_{2}=\frac{1}{2}\int_{R}(u^{2}+u_{x}^{2})\mathrm{d}x,\quad
E_{3}=\frac{1}{2}\int_{R}u(u^{2}+u_{x}^{2})\mathrm{d}x.$

The third example by assigning $\alpha^{2}=c_{2}=c_{3}=1$, $\gamma=c_{0}=0$,
$c_{1}=-2$ is called _the Degasperis–Procesi equation_

$\displaystyle u_{t}-u_{xxt}+4uu_{x}=3u_{x}u_{xx}+uu_{xxx},\quad x\in
R,\;t>0,$ (1.4)

which can be regarded as a model for nonlinear shallow water dynamics [4]. The
frequently discussed invariants are

$\displaystyle H_{1}=\int_{R}(u-u_{xx})\mathrm{d}x,\quad
H_{2}=\int_{R}(u-u_{xx})v\mathrm{d}x,\quad H_{3}=\int_{R}u^{3}\mathrm{d}x,$

where $4v-v_{xx}=u$.

Up to now, there have been thousands of papers focusing on the theoretical and
numerical studies on these three equations. It is worth mentioning that the
invariant-preserving property is a key index of the success for numerical
methods. However, high-order invariants are usually difficult to preserve
numerically. Liu _et al._ also pointed out “it appears a rather difficult task
to preserve all three conservation laws” in [5]. In this work, higher-order
invariants of these equations will be re-derived in view of the energy method,
which may be possible to provide some thoughts for invariant-preserving
numerical methods. Actually, the energy method originated from conservation
laws in physics was first proposed in 1928 by Courant, Friedrichs and Lewy
[6]. From then on, it has been widely applied to the mathematical and
numerical analysis of nonlinear evolution equations. We trust the readers with
[7] instead of a long list of references to relevant works.

The rest of the paper is arranged as follows. In Section 2, combining the
energy method and a skew-adjoint operator, we show the high-order invariants
for the KdV equation, the Camassa–Holm equation and the Degasperis–Procesi
equation, respectively. Then we list several applications for seeking some
high-order invariants of other types of the shallow water wave equations in
Section 3.

##  2 Main results

In what follows, we directly show that $M_{i}$ $(i=1,2,3,4)$, $E_{i}$
$(i=1,2,3)$ and $H_{i}$ $(i=1,2,3)$ are invariants of (1.2), (1.3) and (1.4)
subjected to the periodic boundary conditions based on the energy method,
respectively.

###  2.1 Invariants of the KdV equation

Proof: (I) Multiplying by $1$, $u$ and $(u^{2}+u_{xx})$, respectively, with
(1.2), we have $M_{i}$ $(i=1,2,3)$. In what follows, we show the fourth
invariant $M_{4}$ of the KdV equation by the energy method.

Multiplying both sides of (1.2) by
$2u_{xxxx}+\frac{10}{3}u_{x}^{2}+\frac{20}{3}uu_{xx}+\frac{20}{9}u^{3}$ and
integrating the result, we have

$\displaystyle 0=$
$\displaystyle\int_{R}\Big{(}2u_{xxxx}+\frac{10}{3}u_{x}^{2}+\frac{20}{3}uu_{xx}+\frac{20}{9}u^{3}\Big{)}\cdot
u_{t}\mathrm{d}x$
$\displaystyle-\int_{R}\Big{(}2u_{xxxx}+\frac{10}{3}u_{x}^{2}+\frac{20}{3}uu_{xx}+\frac{20}{9}u^{3}\Big{)}\cdot(4uu_{x}+2u_{xxx})\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{R}\Big{[}2u_{xx}u_{xxt}-\frac{10}{3}(u_{t}u_{x}^{2}+2uu_{x}u_{xt})+\frac{20}{9}u^{3}u_{t}\Big{]}\mathrm{d}x$
$\displaystyle-\int_{R}\Big{(}2u_{xxxx}+\frac{10}{3}u_{x}^{2}+\frac{20}{3}uu_{xx}+\frac{20}{9}u^{3}\Big{)}\cdot(4uu_{x}+2u_{xxx})\mathrm{d}x$
$\displaystyle=$ $\displaystyle\frac{{\rm d}}{{\rm
d}t}M_{4}-8\int_{R}uu_{x}u_{xxxx}\mathrm{d}x-\frac{40}{3}\int_{R}uu_{x}^{3}\mathrm{d}x-\frac{80}{3}\int_{R}u^{2}u_{x}u_{xx}\mathrm{d}x-\frac{80}{9}\int_{R}u^{4}u_{x}\mathrm{d}x$
$\displaystyle-4\int_{R}u_{xxx}u_{xxxx}\mathrm{d}x-\frac{20}{3}\int_{R}u_{xxx}u_{x}^{2}\mathrm{d}x-\frac{40}{3}\int_{R}uu_{xx}u_{xxx}\mathrm{d}x-\frac{40}{9}\int_{R}u^{3}u_{xxx}\mathrm{d}x.$
(2.1)

It remains to check that the sum of all the integral terms in the above
equation is zero. Calculating each term in (2.1) using the integration by
parts, we have

$\displaystyle-8\int_{R}uu_{x}u_{xxxx}\mathrm{d}x=-20\int_{R}u_{x}u_{xx}^{2}\mathrm{d}x,$
(2.2)
$\displaystyle-\frac{80}{3}\int_{R}u^{2}u_{x}u_{xx}\mathrm{d}x=\frac{80}{3}\int_{R}uu_{x}^{3}\mathrm{d}x,$
(2.3) $\displaystyle-\frac{80}{9}\int_{R}u^{4}u_{x}\mathrm{d}x=0,$ (2.4)
$\displaystyle-4\int_{R}u_{xxx}u_{xxxx}\mathrm{d}x=0,$ (2.5)
$\displaystyle-\frac{20}{3}\int_{R}u_{xxx}u_{x}^{2}\mathrm{d}x=\frac{40}{3}\int_{R}u_{x}u_{xx}^{2}\mathrm{d}x,$
(2.6)
$\displaystyle-\frac{40}{3}\int_{R}uu_{xx}u_{xxx}\mathrm{d}x=\frac{20}{3}\int_{R}u_{x}u_{xx}^{2}\mathrm{d}x,$
(2.7)
$\displaystyle-\frac{40}{9}\int_{R}u^{3}u_{xxx}\mathrm{d}x=-\frac{40}{3}\int_{R}uu_{x}^{3}\mathrm{d}x.$
(2.8)

Substituting (2.2)–(2.8) into (2.1), we have $\frac{{\rm d}}{{\rm
d}t}M_{4}=0$, which completes the proof.

###### Remark 1

Suppose the general form of the KdV equation is

$\displaystyle u_{t}-auu_{x}-bu_{xxx}=0,$

and the corresponding high-order invariant

$\displaystyle M(t)=\int_{R}(u_{xx}^{2}-Auu_{x}^{2}+Bu^{4})\mathrm{d}x.$

Using the same method above, we could derive

$\displaystyle\left\\{\begin{array}[]{ll}5a=3Ab,\\\
12Bb=Aa,\end{array}\right.$

which can be rewritten as

$\displaystyle\frac{a}{b}=\frac{3A}{5}=\frac{12B}{A}.$

Therefore, it follows

$\displaystyle A^{2}=20B.$

For instance, when $a=-6$, $b=-1$, we have $A=10$, $B=5$, which deduces to the
KdV equation as

$\displaystyle u_{t}+6uu_{x}+u_{xxx}=0,$

with a fourth-order invariant

$\displaystyle M(t)=\int_{R}(u_{xx}^{2}-10uu_{x}^{2}+5u^{4})\mathrm{d}x.$

###  2.2 Invariants of the Camassa–Holm equation

Proof: Multiplying by $1$ and $u$ on both sides of (1.3), respectively, and
then integrating the results, which implies $E_{1}$ and $E_{2}$ through the
integration by parts. Below, we prove $E_{3}$ by the energy method. Firstly,
noticing that (1.3) can be written with a skew-adjoint operator
$(1-\partial_{xx})^{-1}$ as

$\displaystyle
u_{t}+uu_{x}+\partial_{x}(1-\partial_{xx})^{-1}\Big{(}u^{2}+\frac{1}{2}u_{x}^{2}\Big{)}=0.$

Let $g=(1-\partial_{xx})^{-1}\Big{(}u^{2}+\frac{1}{2}u_{x}^{2}\Big{)}$. Then
we see from the above equation that (1.3) is equivalent to

$\displaystyle u_{t}+uu_{x}+g_{x}=0,$ (2.9) $\displaystyle
g-g_{xx}=u^{2}+\frac{1}{2}u_{x}^{2}.$ (2.10)

Multiplying (2.9) by $3u^{2}+u_{x}^{2}-2(uu_{x})_{x}$ and integrating the
result on both sides, we have

$\displaystyle 0$
$\displaystyle=\int_{R}(u_{t}+uu_{x}+g_{x})\cdot(3u^{2}+u_{x}^{2}-2(uu_{x})_{x})\mathrm{d}x$
$\displaystyle=\int_{R}u_{t}\cdot(3u^{2}+u_{x}^{2}-2(uu_{x})_{x})\mathrm{d}x+\int_{R}(uu_{x}+g_{x})\cdot(3u^{2}+u_{x}^{2}-2(uu_{x})_{x})\mathrm{d}x$
$\displaystyle\triangleq A+B.$ (2.11)

Calculating each term derives that

$\displaystyle A$
$\displaystyle=\int_{R}u_{t}\cdot(3u^{2}+u_{x}^{2}-2(uu_{x})_{x})\mathrm{d}x$
$\displaystyle=\int_{R}u_{t}\cdot(3u^{2}+u_{x}^{2})\mathrm{d}x+\int_{R}2uu_{x}\cdot
u_{xt}\mathrm{d}x$ $\displaystyle=\int_{R}u_{t}\cdot
3u^{2}\mathrm{d}x+\int_{R}u_{t}\cdot
u_{x}^{2}\mathrm{d}x+\int_{R}u\cdot(u_{x}^{2})_{t}\mathrm{d}x$
$\displaystyle=\int_{R}(u^{3})_{t}\mathrm{d}x+\int_{R}(u\cdot
u_{x}^{2})_{t}\mathrm{d}x$ $\displaystyle=\frac{\rm d}{{\rm
d}t}\int_{R}(u^{3}+uu_{x}^{2})\mathrm{d}x$ (2.12)

and

$\displaystyle B$
$\displaystyle=\int_{R}(uu_{x}+g_{x})\cdot(3u^{2}+u_{x}^{2}-2(uu_{x})_{x})\mathrm{d}x$
$\displaystyle=\int_{R}u\cdot
u_{x}^{3}\mathrm{d}x+\int_{R}g_{x}\cdot(3u^{2}+u_{x}^{2})\mathrm{d}x-\int_{R}g_{x}\cdot
2(uu_{x})_{x}\mathrm{d}x$ $\displaystyle=\int_{R}u\cdot
u_{x}^{3}\mathrm{d}x+\int_{R}g_{x}\cdot(3u^{2}+u_{x}^{2})\mathrm{d}x+2\int_{R}g_{xx}\cdot
uu_{x}\mathrm{d}x$ $\displaystyle=\int_{R}u\cdot
u_{x}^{3}\mathrm{d}x+\int_{R}g_{x}\cdot(3u^{2}+u_{x}^{2})\mathrm{d}x+2\int_{R}(g-u^{2}-\frac{1}{2}u_{x}^{2})\cdot
uu_{x}\mathrm{d}x$
$\displaystyle=\int_{R}g_{x}\cdot(3u^{2}+u_{x}^{2})\mathrm{d}x+2\int_{R}g\cdot
uu_{x}\mathrm{d}x$
$\displaystyle=\int_{R}g_{x}\cdot(3u^{2}+u_{x}^{2})\mathrm{d}x-\int_{R}g_{x}\cdot
u^{2}\mathrm{d}x$
$\displaystyle=\int_{R}g_{x}\cdot(2u^{2}+u_{x}^{2})\mathrm{d}x$
$\displaystyle=2\int_{R}g_{x}\cdot(g-g_{xx})\mathrm{d}x=0.$ (2.13)

Substituting (2.12) and (2.13) into (2.11), we have

$\displaystyle\frac{\rm d}{{\rm d}t}\int_{R}(u^{3}+uu_{x}^{2})\mathrm{d}x=0,$

which implies $E_{3}$.

###  2.3 Invariants of the Degasperis–Procesi equation

Proof: Integrating on both sides of (1.4), it easily obtains $H_{1}$. Then we
show invariants $H_{2}$ and $H_{3}$ of (1.4), respectively. Firstly let
$g=(1-\partial_{xx})^{-1}\Big{(}\frac{3}{2}u^{2}\Big{)}$, then (1.4) is
equivalent to

$\displaystyle u_{t}+uu_{x}+g_{x}=0,$ (2.14) $\displaystyle
g-g_{xx}=\frac{3}{2}u^{2}.$ (2.15)

Multiplying by $2u-6v$ on both sides of (2.14) and then integrating the
result, we have

$\displaystyle 0$
$\displaystyle=\int_{R}(u_{t}+uu_{x}+g_{x})\cdot(2u-6v)\mathrm{d}x$
$\displaystyle=\int_{R}u_{t}\cdot(2u-6v)\mathrm{d}x+\int_{R}uu_{x}\cdot(2u-6v)\mathrm{d}x+\int_{R}g_{x}\cdot(2u-6v)\mathrm{d}x$
$\displaystyle\triangleq C+D.$ (2.16)

The each term in the above identity is estimated as

$\displaystyle C$
$\displaystyle=\int_{R}u_{t}\cdot(2u-6v)\mathrm{d}x=2\int_{R}u_{t}\cdot
u\mathrm{d}x-6\int_{R}u_{t}\cdot v\mathrm{d}x=2\int_{R}u_{t}\cdot
u\mathrm{d}x-6\int_{R}(4v_{t}-v_{xxt})\cdot v\mathrm{d}x$
$\displaystyle=2\int_{R}u_{t}\cdot u\mathrm{d}x-24\int_{R}v_{t}\cdot
v\mathrm{d}x-6\int_{R}v_{xt}\cdot v_{x}\mathrm{d}x=\frac{{\rm d}}{{\rm
d}t}\int_{R}(u^{2}-12v^{2}-3v_{x}^{2})\mathrm{d}x$ $\displaystyle=\frac{{\rm
d}}{{\rm d}t}\int_{R}\Big{(}u^{2}-3(4v-v_{xx})\cdot
v\Big{)}\mathrm{d}x=\frac{{\rm d}}{{\rm
d}t}\int_{R}(u^{2}-3uv)\mathrm{d}x=\frac{{\rm d}}{{\rm
d}t}\int_{R}u\cdot(u-3v)\mathrm{d}x$ $\displaystyle=\frac{{\rm d}}{{\rm
d}t}\int_{R}u\cdot(v-v_{xx})\mathrm{d}x=\frac{{\rm d}}{{\rm
d}t}\int_{R}(u-u_{xx})\cdot v\mathrm{d}x$ (2.17)

and

$\displaystyle D$
$\displaystyle=\int_{R}uu_{x}\cdot(2u-6v)\mathrm{d}x+\int_{R}g_{x}\cdot(2u-6v)\mathrm{d}x$
$\displaystyle=-6\int_{R}uu_{x}\cdot
v\mathrm{d}x+\int_{R}g_{x}\cdot(2u-6v)\mathrm{d}x$
$\displaystyle=3\int_{R}u^{2}\cdot
v_{x}\mathrm{d}x+\int_{R}g_{x}\cdot(2u-6v)\mathrm{d}x$
$\displaystyle=2\int_{R}(g-g_{xx})\cdot
v_{x}\mathrm{d}x+\int_{R}g_{x}\cdot(2u-6v)\mathrm{d}x$
$\displaystyle=2\int_{R}g\cdot v_{x}\mathrm{d}x-2\int_{R}g_{xx}\cdot
v_{x}\mathrm{d}x+\int_{R}g_{x}\cdot(2v-2v_{xx})\mathrm{d}x$
$\displaystyle=2\int_{R}g\cdot v_{x}\mathrm{d}x+2\int_{R}g_{x}\cdot
v\mathrm{d}x-2\int_{R}g_{xx}\cdot v_{x}\mathrm{d}x-2\int_{R}g_{x}\cdot
v_{xx}\mathrm{d}x$
$\displaystyle=2\int_{R}(gv)_{x}\mathrm{d}x-2\int_{R}(g_{x}\cdot
v_{x})_{x}\mathrm{d}x=0.$ (2.18)

Substituting (2.17) and (2.18) into (2.16), we have

$\displaystyle\frac{{\rm d}}{{\rm d}t}\int_{R}(u-u_{xx})\cdot v\mathrm{d}x=0,$

which implies $H_{2}$.

Finally, we show $H_{3}$. Multiplying (2.14) on both sides by $u^{2}$ and
integrating the result, it yields by noting (2.15)

$\displaystyle 0$ $\displaystyle=\int_{R}(u_{t}+uu_{x}+g_{x})\cdot
u^{2}\mathrm{d}x$ $\displaystyle=\int_{R}u_{t}\cdot
u^{2}\mathrm{d}x+\int_{R}u^{3}\cdot u_{x}\mathrm{d}x+\int_{R}g_{x}\cdot
u^{2}\mathrm{d}x$
$\displaystyle=\int_{R}\Big{(}\frac{1}{3}u^{3}\Big{)}_{t}\mathrm{d}x+\frac{2}{3}\int
g_{x}\cdot(g-g_{xx})\mathrm{d}x$ $\displaystyle=\frac{1}{3}\frac{{\rm d}}{{\rm
d}t}\int_{R}u^{3}\mathrm{d}x,$

which implies the invariant $H_{3}$.

##  3 Applications to other periodic nonlinear dispersive waves

###  3.1 Benjamin-Bona-Mahony equation

Consider the Benjamin-Bona-Mahony equation [8] of the form

$\displaystyle u_{t}-u_{xxt}+u_{x}+\varepsilon uu_{x}=0,\quad x\in R.$ (3.1)

It can be written as

$\displaystyle
u_{t}+\partial_{x}(1-\partial_{xx})^{-1}\Big{(}u+\frac{\varepsilon}{2}u^{2}\Big{)}=0,\quad
x\in R.$

Let $g=(1-\partial_{xx})^{-1}\Big{(}u+\frac{\varepsilon}{2}u^{2}\Big{)}$, then
the equation (3.1) turns out to be

$\displaystyle u_{t}+g_{x}=0,$ (3.2) $\displaystyle
g-g_{xx}=u+\frac{\varepsilon}{2}u^{2}.$ (3.3)

Multiplying both sides of (3.2) by $u^{2}$ and integrating the result, and
then using (3.3), we have

$\displaystyle 0$ $\displaystyle=\int_{R}(u_{t}+g_{x})\cdot
u^{2}\mathrm{d}x=\int_{R}u_{t}\cdot u^{2}\mathrm{d}x+\int_{R}g_{x}\cdot
u^{2}\mathrm{d}x$ $\displaystyle=\int_{R}u_{t}\cdot
u^{2}\mathrm{d}x+\frac{2}{\varepsilon}\int_{R}g_{x}\cdot(g-g_{xx}-u)\mathrm{d}x=\int_{R}u_{t}\cdot
u^{2}\mathrm{d}x-\frac{2}{\varepsilon}\int_{R}g_{x}\cdot u\mathrm{d}x$
$\displaystyle=\int_{R}u_{t}\cdot
u^{2}\mathrm{d}x+\frac{2}{\varepsilon}\int_{R}u_{t}\cdot
u\mathrm{d}x=\frac{{\rm d}}{{\rm
d}t}\int_{R}\Big{(}\frac{1}{3}u^{3}+\frac{1}{\varepsilon}u^{2}\Big{)}\mathrm{d}x,$

which indicates

$\displaystyle\int_{R}\frac{1}{3}\Big{(}u^{3}+\frac{1}{\varepsilon}u^{2}\Big{)}\mathrm{d}x$

is a three-order invariant for (3.1).

###  3.2 Regularized long wave equation

Consider the regularized long wave equation [9] of the form

$\displaystyle u_{t}-\mu u_{xxt}+u_{x}+u^{p}u_{x}=0,$ (3.4)

where $\mu>0$ is a positive constant. When $p=2$, it is called modified
regularized long wave equation; when $p\geqslant 3$, it is called generalized
regularized long wave equation. Similar to the foregoing argument, (3.4) can
be written as an equivalent form of

$\displaystyle u_{t}+g_{x}=0,$ (3.5) $\displaystyle g-\mu
g_{xx}=u+\frac{1}{p+1}u^{p+1}.$ (3.6)

Multiplying both sides of (3.5) by $u^{p+1}$, integrating the result, and then
using (3.6), we have

$\displaystyle 0$ $\displaystyle=\int_{R}(u_{t}+g_{x})\cdot
u^{p+1}\mathrm{d}x=\int_{R}u_{t}\cdot u^{p+1}\mathrm{d}x+\int_{R}g_{x}\cdot
u^{p+1}\mathrm{d}x$ $\displaystyle=\int_{R}u_{t}\cdot
u^{p+1}\mathrm{d}x+({p+1})\int_{R}g_{x}\cdot(g-\mu g_{xx}-u)\mathrm{d}x$
$\displaystyle=\int_{R}u_{t}\cdot u^{p+1}\mathrm{d}x-({p+1})\int_{R}g_{x}\cdot
u\mathrm{d}x$ $\displaystyle=\int_{R}u_{t}\cdot
u^{p+1}\mathrm{d}x+({p+1})\int_{R}u_{t}\cdot u\mathrm{d}x$
$\displaystyle=\frac{{\rm d}}{{\rm
d}t}\int_{R}\Big{(}\frac{1}{p+2}u^{p+2}+\frac{{p+1}}{2}u^{2}\Big{)}\mathrm{d}x,$

which indicates

$\displaystyle\int_{R}\Big{(}\frac{1}{p+2}u^{p+2}+\frac{{p+1}}{2}u^{2}\Big{)}\mathrm{d}x$

is a high-order invariant for (3.4). This corrects an invariant $I_{3}$ in
Example 4 appeared in [10] (pp. 492).

###  3.3 Rosenau equation

Consider the Rosenau equation [11]

$\displaystyle u_{t}+u_{xxxxt}+u_{x}+uu_{x}=0,$ (3.7)

which is equivalent to

$\displaystyle u_{t}+g_{x}=0,$ (3.8) $\displaystyle
g+g_{xxxx}=u+\frac{1}{2}u^{2}.$ (3.9)

Multiplying both sides of (3.8) by $u^{2}$ and noticing (3.9), similar to the
argument in the above, we have a third-order invariant for (3.7) of the form

$\displaystyle\int_{R}\Big{(}\frac{1}{3}u^{3}+u^{2}\Big{)}\mathrm{d}x.$

## Acknowledgement

We appreciate Prof. Zhi-zhong Sun for many useful discussions. This work is
dedicated to Prof. Zhi-zhong Sun on the occasion of his 60th birthday. The
work is supported by Natural Science Foundation of Zhejiang Province (Grant
No. LZ23A010007).

## References

## References

  * [1] J. Escher, Y. Liu, Z. Yin, Global weak solutions and blow-up structure for the Degasperis-Procesi equation. J. Funct. Anal., 241 (2006) 457–485.
  * [2] T. Tao, Low-regularity global solutions to nonlinear dispersive equations. Surveys in analysis and operator theory (Canberra, 2001), 19–48, Proc. Centre Math. Appl. Austral. Nat. Univ., 40, Austral. Nat. Univ., Canberra, (2002).
  * [3] R. Camassa, D. D. Holm, An integrable shallow water equation with peaked solitons. Phys. Rev. Lett., 71 (1993) 1661–1664.
  * [4] A. Degasperis, M. Procesi, Asymptotic integrability, in: A. Degasperis, G. Gaeta (Eds.). Symmetry and Perturbation Theory, World Scientific, Singapore, (1999) 23–37.
  * [5] H. Liu, Y. Xing, An invariant preserving discontinuous Galerkin method for the Camassa-Holm equation. SIAM J. Sci. Comput., 38 (2016) A1919–A1934.
  * [6] R. Courant, K.O. Friedrichs, H. Lewy, Über die partiellen Differenzenglei-chungen der mathematischen physik, Math. Ann., 100 (1928) 32–74.
  * [7] Z. Sun. Finite Difference Methods for Nonlinear Evolution Equations, Science Press, Beijing, (2018).
  * [8] L.A. Medeiros, G.P. Menzala, Existence and uniqueness for periodic solutions of the Benjamin-Bona-Mahony equation. SIAM J. Math. Anal., 8(5) (1977) 792–799.
  * [9] C.E. Seyler, D.L. Fenstermacher, A symmetric regularized-long-wave equation. The Physics of Fluids, 27(4) (1984) 4–7.
  * [10] A. Ghiloufi, K. Omrani, New conservative difference schemes with fourth-order accuracy for some model equation for nonlinear dispersive waves. Numer. Methods Partial Differential Equation, 34 (2018) 451–500.
  * [11] M.A. Park, On the Rosenau equation. Math. Appl. Comput., 9 (1990) 145–152.


# Self-Supervised Modality-Agnostic Pre-Training of Swin Transformers

###### Abstract

Unsupervised pre-training has emerged as a transformative paradigm, displaying
remarkable advancements in various domains. However, the susceptibility to
domain shift, where pre-training data distribution differs from fine-tuning,
poses a significant obstacle. To address this, we augment the Swin Transformer
to learn from different medical imaging modalities, enhancing downstream
performance. Our model, dubbed SwinFUSE (Swin Multi-Modal Fusion for
UnSupervised Enhancement), offers three key advantages: (i) it learns from
both Computed Tomography (CT) and Magnetic Resonance Images (MRI) during pre-
training, resulting in complementary feature representations; (ii) a domain-
invariance module (DIM) that effectively highlights salient input regions,
enhancing adaptability; (iii) exhibits remarkable generalizability, surpassing
the confines of tasks it was initially pre-trained on. Our experiments on two
publicly available 3D segmentation datasets show a modest 1-2% performance
trade-off compared to single-modality models, yet significant out-performance
of up to 27% on out-of-distribution modality. This substantial improvement
underscores our proposed approach’s practical relevance and real-world
applicability. Code is available at: https://github.com/devalab/SwinFUSE

Index Terms—  self-supervision, multi-modal, domain adaptation, 3D image
segmentation

##  1 Introduction

Supervised deep learning excels at medical image segmentation using lots of
labeled data [1, 2, 3]. The shortage of professional radiologists and their
limited time and annotation efficiency makes it difficult to get huge medical
picture collections with exact annotations. Thus, routine clinical usage of
supervised-learning-based segmentation techniques is limited. Recent research
has focused on self-supervised learning (SSL), which uses many unlabeled
images to learn the general aspects of medical images. Fully supervised model
fine-tuning uses a minimal quantity of labeled data [4]. Effective self-
supervised medical image segmentation depends on pre-training quality.

Contrastive learning is a self-supervised pre-training method that minimizes
the latent space distance of pairs of similar images (typically produced from
the same original image using different data augmentation processes) and
maximizes the distance of pairs of dissimilar ones [5]. Contrastive learning
methods address domain shift, enhancing model applicability in downstream
segmentation networks through consistent data augmentation strategies ensuring
similar input distributions. Existing self-supervised learning has two
limitations [6].

Fig. 1: Visual interpretation of SwinFUSE’s attention weights (darker shades
indicate higher relevance) for a BraTS21 MRI and the model’s segmentation
output.

  * •

Domain Shift: Upstream pre-training uses modified images, affecting downstream
segmentation network input distributions. General features from pre-trained
models may not apply to segmentation networks.

  * •

Multi-modality: Current techniques often rely on single-modal data, missing
the benefits of multiple modalities. Multi-modal images offer diverse
perspectives and augment network segmentation information.

Fig. 2: Outline of our proposed pre-training pipeline. Sub-volumes are
randomly created from input images and augmented with random inner cutouts and
rotations ($x_{i},x_{j}$). Each augmentation passes through the patch
partition layer to generate embeddings, which are fed to the DIM. The output
from the DIM is extracted as kernel densities and forwarded to the Swin
Transformer.

Vision Transformers (ViTs) transformed medical image analysis and computer
vision. Transformers thrive in pre-text tasks, large-scale training, and
layer-based global and local knowledge learning. ViTs simulate long-range
global information using self-attention blocks and encode visual
representations from patches, unlike Convolutional Neural Networks (CNNs) with
small receptive fields. A hierarchical ViT with Shifted Windows (Swin) for
local self-attention computing with non-overlapping windows was developed by
Liu et al. [7]. Linear architecture has been found to be more efficient than
ViT’s quadratic self-attention layers. Swin UNETR [8] merges feature maps at
various sizes using transformer-encoded spatial representations in
convolution-based decoders and achieves state-of-the-art (SOTA) performance in
BTCV multi-organ segmentation and Medical Segmentation Decathlon (MSD)
challenges [9]. The training paradigm uses proxy activities to learn human
anatomical patterns. We extend this intuition to allow complementary feature
learning from multiple imaging modalities.

This paper presents SwinFUSE, a modality-invariant self-supervised pre-
training approach for medical image analysis. We utilize Swin UNETR’s
contrastive learning, masked volume pinpointing, and 3D rotation prediction as
proxy tasks for pre-training. Additionally, we introduce a DIM for concurrent
feature learning from CT and MRI data. The DIM implicitly identifies relevant
input areas and directs them to the Swin Transformer encoder using attention
maps (Fig. 1). We train the network on the SynthRad dataset [10] and retain
the DIM and encoder for later fine-tuning. Our 3D image segmentation
experiments involve the BraTS21 [11] and MSD datasets. We fine-tune the entire
network to demonstrate generalization across domains, validating its efficacy
for each task, including organ segmentation in the MSD dataset.

##  2 Method

###  2.1 Datasets

#### SynthRAD

comprises registered brain and pelvis CT images with cone-beam CT and MRI
images, serving the purpose of synthetic CT generation for radiotherapy
planning [10]. We focus on a subset of 180 patients, utilizing T1-weighted
gradient-echo MRIs, with some using contrast.

#### BraTS21

consists of multi-modal MRI scans of glioma, with a total of 1254 patients
[11]. The sequences acquired include T1, T2, T1CE, and FLAIR. Segmentation
classes include peritumoral edematous/invaded tissue, tumor core, and necrotic
tumor core.

#### MSD

contains 2,633 3D images collected from various anatomical regions,
modalities, and medical image sources for segmentation purposes [9]. It covers
data on body organs or parts like the Brain, Heart, Liver, Lung, Pancreas,
Prostate, Hepatic Vessel, Hippocampus, Spleen, and Colon.

###  2.2 Pre-training

We augment the Swin UNETR architecture using a novel Domain Invariance Module,
trained to learn which features to highlight, conditioned on the input type.
The training dataset consists of CTs and MRIs; volumes are sampled randomly.
During each iteration, 3D patches measuring $x_{n}\in\mathbb{R}^{96\times
96\times 96}$ voxels undergo augmentation with random inner cutout and
rotation. These patches are then projected into a C-dimensional space ($C=48$)
using an embedding layer leading to the DIM as shown in Fig. 2. The two
embeddings from respective augmentation are each fed into a 4-layer deep
Multi-Head Attention Block (MHA) with $3,6,12,24$ heads respectively like Co-
Attention [12]. The $Q,K,V$ embedding dimensions increase by an exponent of
$2$ with the base layer having dimensions $x\in\mathbb{R}^{48}$. The query
vector from the first embedding is fed as query input to the second MHA block
and vice versa. Each block in the DIM is initialized with an embedding
dimension of $2304$. The attention weights are scaled with the original
embeddings, and the resulting average embedding is sent as input to the Swin
Transformer Encoder. The DIM is constructed as given below, where $P$ denotes
the patch partitions being fed into each MHA block ($\Psi$).

$DIM:P_{1}\cdot\Psi_{1}(Q_{2},K_{1},V_{1})+\\\
P_{2}\cdot\Psi_{2}(Q_{1},K_{2},V_{2})$ (1)

We train the model on the SynthRad dataset using the AdamW optimizer and a
warm-up cosine scheduler with 500 iterations on two RTX 3090’s. We set the
initial learning rate for the pre-training experiments to $4e^{-4}$ and a
decay of $1e^{-5}$. We implement our model using PyTorch and MONAI
111https://monai.io/.

###  2.3 Loss Function

We aim to minimize the loss of Swin UNETR’s encoder using multiple pre-
training objectives, including masked volume inpainting, 3D image rotation,
and contrastive coding. Additionally, we maximize an extra loss term that,
akin to the approach in [13], employs non-parametric density estimation
through kernel density estimation (KDE) and density matching via Jenson-
Shannon divergence (JSD). This density-matching loss is a regularizer,
ensuring that the feature distribution overlap between the source and target
datasets is minimized. The KDE, denoted as $p_{est}(X)$, is formulated as
follows:

$p_{est}(X)=\frac{1}{N}\sum_{n=1}^{N}K\left(\frac{\lVert
X-X_{n}\rVert_{2}}{\sigma}\right)$ (2)

where $X_{1},X_{2},X_{3},\cdots,X_{N}$ is the number of sampled points from
the encoded feature space, the output from the MHA block in our model, and $K$
is a Gaussian kernel. The bandwidth parameter ($\sigma$) is estimated to be
the mean of the distance between the nearest neighbors in the feature space.

Our loss term for density matching, $\mathcal{L}_{\text{JSD}}$, given density
of each MHA block outputs as $p_{1}$ and $p_{2}$, is given as follows:

$JSD_{p_{1},p_{2}}=\frac{1}{2}\\{KL[p_{1},M]+KL[p_{2},M]\\}$ (3)

where KL is the KL divergence between the two distributions and $M$ is the
average of both density estimates. The final loss:

$\mathcal{L}_{\text{total}}=\mathcal{L}_{\text{inpaint}}+\mathcal{L}_{\text{contrast}}+\mathcal{L}_{\text{rot}}-\mathcal{L}_{\text{JSD}}$
(4)

###  2.4 Fine-tuning

In the downstream task, such as 3D image segmentation, we fine-tune the
complete Swin UNETR model by removing the projection heads while retaining the
DIM. During training, sub-volumes are randomly cropped from the volumetric
data. Then, stochastic data augmentations, including random rotation and
cutout, are applied twice to each sub-volume within a mini-batch, resulting in
two different views of each data. All other augmentation parameters align with
those used in Swin UNETR.

For SwinFUSE, we utilize pre-trained weights for both the CT and MRI tasks,
following the official methods outlined in nnUnet [14] and Swin UNETR [8]. To
ensure robustness, we employ a five-fold cross-validation strategy to train
models for BraTS21 and MSD experiments. In each fold, we select the best model
and ensemble their outputs to generate the final segmentation predictions.

Fig. 3: Qualitative visualizations of Swin UNETR and our proposed method.
Colored regions correspond to necrotic tumor core (red), peritumoral edematous
tissue (pink), and enhancing tumor (blue). Dice scores are also given.

##  3 Results

Table 1: Average Dice Score across five folds on MSD for two variants of SwinFUSE: with and without pre-training denoted by the presence and absence of (P) respectively. Organ | Brain | Heart | Hippocampus | Liver | Lung | Pancreas | Prostate | Colon |  | Hepatic
---
Vessel
Spleen
nnUNET [14] | 64.50 | 94.82 | 89.76 | 86.67 | 72.78 | 68.31 | 84.14 | 59.43 | 70.14 | 96.34
Swin UNETR [8] | 66.31 | 94.32 | 88.39 | 87.42 | 76.40 | 72.91 | 81.51 | 60.35 | 70.75 | 95.79
SwinFUSE | 65.17 | 94.67 | 86.76 | 88.31 | 74.39 | 70.95 | 80.63 | 60.42 | 69.78 | 94.54
SwinFUSE (P) | 66.34 | 94.91 | 87.91 | 89.43 | 75.35 | 69.25 | 81.62 | 59.31 | 71.61 | 96.24

#### Quantitative

Average dice scores across five folds for each task in MSD are detailed in
Table. 1 for nnUnet, Swin UNETR, and our model (with and without pre-trained
weights). The proposed method outperforms the current SOTA while segmenting
the Brain, Heart, Liver, and Hepatic Vessels but, on average, is 1-2% worse
than Swin UNETR. This might be because we reuse the same pre-trained weights
to fine-tune SwinFUSE for both CT and MRI tasks, differing from Swin UNETR,
which uses pre-trained weights for only CT tasks.

#### Qualitative

We visualize images from the MSD (heart, hippocampus, prostate, brain) and
BraTS21 datasets in Fig. 3. Segmentation outputs from SwinFUSE are more
concise and perform well in the global context. Moderate improvements are seen
for smaller organs like the hippocampus (second row). In contrast, for organs
like the brain (last row), we notice that the attention mechanism helps locate
the unconnected region in the right hemisphere, which Swim UNETR completely
fails to detect. In Fig. 1 we generate the learned attention weights from
which we can see that the DIM learns similarities between CTs and MRIs and
uses those to anchor itself. It further uses differentiating aspects between
the both to effectively highlight regions before sending the input to the Swin
transformer.

Table 2: Performance on MRI tasks after pre-training on CT organ regions. Dice scores of Swin UNETR (1) and our model (2) are reported in the format 1/2 in each cell Fine- tuned on | Tested on
---|---
Brain | Heart | Prostate |  | Hippo-
---
campus
Liver | 0.22/0.47 | 0.31/0.52 | 0.19/0.32 | 0.18/0.33
Pancreas | 0.21/0.29 | 0.27/0.34 | 0.20/0.40 | 0.16/0.26
Lung | 0.15/0.32 | 0.30/0.46 | 0.22/0.34 | 0.23/0.30
Spleen | 0.25/0.32 | 0.21/0.39 | 0.26/0.44 | 0.25/0.48
Colon | 0.19/0.39 | 0.24/0.53 | 0.29/0.49 | 0.27/0.38
| Hepatic
---
Vessel
0.17/0.35 | 0.22/0.47 | 0.29/0.43 | 0.25/0.41

#### Out-of-distribution

When we fine-tune SwinFUSE and Swin UNETR on CT organ regions in MSD like
Liver, Lung, and Pancreas, and later test on MRI regions, we notice a
significant drop in performance for the latter. For example, when fine-tuned
on the Liver and tested on the Brain, Swin UNETR’s average dice score is 0.22,
whereas SwinFUSE’s is 0.47. In another instance, fine-tuning on the Pancreas
and testing on the Prostate resulted in Swin UNETR achieving a score of 0.16
compared to 0.26 for ours. Multiple other experiments like these are shown in
Table. 2 where we showcase a minimum improvement of 7% and a maximum of 27% in
dice scores. In conclusion, due to the pre-training of our model on a varied
collection of human body compositions and its ability to acquire a versatile
representation from data obtained from various institutions, we assert that
our model is more suitable for clinical applications than existing single-
modality models.

##  4 Discussion and Conclusion

In our research, we have shown noteworthy improvements in the field of medical
imaging by using unsupervised pre-training. However, we acknowledge a
challenge known as domain shift, where the data used for pre-training differs
from the data used for fine-tuning. To address this, we extended the Swin
Transformer framework to pre-train SwinFUSE on two distinct medical imaging
modalities, CT and MRI. This extension offers three key advantages:

  * •

Complementary Feature Representations: By training on both CT and MRI data,
SwinFUSE learns diverse feature representations, making it more adaptable and
robust.

  * •

Domain-Invariance Module: Our DIM helps SwinFUSE adapt to domain shifts by
emphasizing important regions.

  * •

Remarkable Generalizability: SwinFUSE can perform well on tasks it wasn’t
initially trained for, making it highly relevant in real-world applications.

Our experiments show that our approach performs slightly worse than single-
modality models on in-distribution tasks but significantly outperforms them on
out-of-distribution modalities, highlighting its practical applicability.
Quantitatively, our method surpassed state-of-the-art models in segmenting the
Brain, Heart, Liver, and Hepatic Vessels in the MSD dataset. We emphasize that
our diverse pre-training data and versatile representations make SwinFUSE more
suitable for clinical use than single-modality models.

We’ve demonstrated the advantages of pre-training SwinFUSE using our domain-
invariance module and its superior performance on various tasks. We’ve also
highlighted the potential for future research in addressing domain gaps and
applying our framework to other medical imaging modalities like PET and
X-rays. This work represents a significant advancement in the efficiency and
accuracy of medical image analysis, with promising prospects for future
research.

##  5 Compliance with ethical standards

This research study was conducted retrospectively using human subject data
available in open access by RSNA ASNR-MICCAI BraTS 2021, SynthRAD 2023, and
Medical Segmentation Decathlon challenges. Ethical approval was not required,
as confirmed by the license attached to the open access data.

##  6 Acknowledgments

We thank IHub-Data, International Institute of Information Technology,
Hyderabad, for their support.

## References

  * [1] Qi Dou, Lequan Yu, Hao Chen, Yueming Jin, Xin Yang, Jing Qin, and Pheng-Ann Heng,  “3d deeply supervised network for automated segmentation of volumetric medical images,”  Medical image analysis, vol. 41, pp. 40–54, 2017.
  * [2] Zhenghua Xu, Shijie Liu, Di Yuan, Lei Wang, Junyang Chen, Thomas Lukasiewicz, Zhigang Fu, and Rui Zhang,  “$\omega$-net: Dual supervised medical image segmentation with multi-dimensional self-attention and diversely-connected multi-scale convolution,”  Neurocomputing, vol. 500, pp. 177–190, 2022.
  * [3] Tuan Truong, Sadegh Mohammadi, and Matthias Lenga,  “How transferable are self-supervised features in medical image classification tasks?,”  in Machine Learning for Health. PMLR, 2021, pp. 54–74.
  * [4] Xinrui Zhuang, Yuexiang Li, Yifan Hu, Kai Ma, Yujiu Yang, and Yefeng Zheng,  “Self-supervised feature learning for 3d medical images by playing a rubik’s cube,”  in Medical Image Computing and Computer Assisted Intervention–MICCAI 2019: 22nd International Conference, Shenzhen, China, October 13–17, 2019, Proceedings, Part IV 22. Springer, 2019, pp. 420–428.
  * [5] Ting Chen, Simon Kornblith, Kevin Swersky, Mohammad Norouzi, and Geoffrey E Hinton,  “Big self-supervised models are strong semi-supervised learners,”  Advances in neural information processing systems, vol. 33, pp. 22243–22255, 2020.
  * [6] Jiaojiao Zhang, Shuo Zhang, Xiaoqian Shen, Thomas Lukasiewicz, and Zhenghua Xu,  “Multi-condos: Multimodal contrastive domain sharing generative adversarial networks for self-supervised medical image segmentation,”  IEEE Transactions on Medical Imaging, 2023.
  * [7] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo,  “Swin transformer: Hierarchical vision transformer using shifted windows,”  in Proceedings of the IEEE/CVF international conference on computer vision, 2021, pp. 10012–10022.
  * [8] Ali Hatamizadeh, Vishwesh Nath, Yucheng Tang, Dong Yang, Holger R Roth, and Daguang Xu,  “Swin unetr: Swin transformers for semantic segmentation of brain tumors in mri images,”  in International MICCAI Brainlesion Workshop. Springer, 2021, pp. 272–284.
  * [9] Michela Antonelli, Annika Reinke, Spyridon Bakas, Keyvan Farahani, Annette Kopp-Schneider, Bennett A Landman, Geert Litjens, Bjoern Menze, Olaf Ronneberger, Ronald M Summers, et al.,  “The medical segmentation decathlon,”  Nature communications, vol. 13, no. 1, pp. 4128, 2022.
  * [10] Adrian Thummerer, Erik van der Bijl, Arthur Galapon Jr, Joost J. C. Verhoeff, Johannes A. Langendijk, Stefan Both, Cornelis (Nico) A. T. van den Berg, and Matteo Maspero,  “Synthrad2023 grand challenge dataset: Generating synthetic ct for radiotherapy,”  Medical Physics, vol. 50, no. 7, pp. 4664–4674, 2023.
  * [11] Ujjwal Baid, Satyam Ghodasara, Suyash Mohan, Michel Bilello, Evan Calabrese, Errol Colak, Keyvan Farahani, Jayashree Kalpathy-Cramer, Felipe C Kitamura, Sarthak Pati, et al.,  “The rsna-asnr-miccai brats 2021 benchmark on brain tumor segmentation and radiogenomic classification,”  arXiv preprint arXiv:2107.02314, 2021.
  * [12] Mennatullah Siam, Naren Doraiswamy, Boris N Oreshkin, Hengshuai Yao, and Martin Jagersand,  “Weakly supervised few-shot object segmentation using co-attention with visual and semantic embeddings,”  arXiv preprint arXiv:2001.09540, 2020.
  * [13] Qingsong Xie, Yuexiang Li, Nanjun He, Munan Ning, Kai Ma, Guoxing Wang, Yong Lian, and Yefeng Zheng,  “Unsupervised domain adaptation for medical image segmentation by disentanglement learning and self-training,”  IEEE Transactions on Medical Imaging, 2022.
  * [14] Fabian Isensee, Jens Petersen, Simon AA Kohl, Paul F Jäger, and Klaus H Maier-Hein,  “nnu-net: Breaking the spell on successful medical image segmentation,”  arXiv preprint arXiv:1904.08128, vol. 1, no. 1-8, pp. 2, 2019.