2.1. The BK system

When the BK system is used to model water waves over a flat bottom it has the following form [Reference Broer7, Eqs. (3.2)–(3.3)], [Reference Nabelek and Zakharov29, Eqs. (1)–(2)].

(2.1)\begin{gather} \begin{bmatrix} \eta\\\upsilon \end{bmatrix}_{\tilde{t}} = -\begin{bmatrix} (h_0+\eta)\upsilon + \alpha h_0 \upsilon_{\tilde{x}\tilde{x}}\\ g\eta + \tfrac 12 \upsilon^2 \end{bmatrix}_{\tilde{x}} , \end{gather}

where $\tilde{x}$ [m] is the horizontal position, $\tilde{t}$ [s] is the time, ${\mathop{{\eta}\mathopen{\left(\mathord{\tilde{x},\tilde{t}}\right)}}} $ [m] is the free surface elevation with respect to the still water level, ${\mathop{{\upsilon}\mathopen{\left(\mathord{{\tilde{x},\tilde{t}}}\right)}}}$ [m s⁻¹] is the horizontal particle velocity at the free surface, h ₀ [m] is the still water level, τ [N m⁻¹] is the surface tension of the water, ρ [kg m⁻³] is the water density, g [m s⁻²] is the gravitational acceleration, and α [m²] $:= h_0^2/ 3 - \tau/(\rho g).$

Note that subscripts containing function arguments indicate partial derivatives. For example, $\upsilon_{\tilde{x}\tilde{x}}=\frac{\partial^2 \upsilon}{\partial \tilde{x}^2}$. Since the discrete spectrum of capillary waves differs from ordinary, non-capillary waves [Reference Kaup18, §2], we will assume non-capillary waves throughout this paper. That is, $\alpha \gt 0\,\textrm{m}^2 $.

The normalization

(2.2)\begin{align} w & := (h_0+\eta)/ h_0 \geq 0, \end{align}

(2.3)\begin{align} u & := \upsilon / \sqrt{gh_0} \in\mathbb{R}, \end{align}

(2.4)\begin{align} x & := \tilde{x} / \sqrt{\alpha} \in\mathbb{R}, \end{align}

(2.5)\begin{align} t & := \tilde{t} / \sqrt{\alpha/(gh_0)} \in\mathbb{R} \end{align}

results in the normalized BK system [Reference Sachs39, Eq. (2)]

(2.6)\begin{gather} \begin{bmatrix} w\\u \end{bmatrix}_{t} = -\begin{bmatrix} wu + u_{xx}\\ w + \tfrac 12 u^2 \end{bmatrix}_{x} . \end{gather}

The same normalized equation can be obtained for the formulation of the BK system in [Reference Kaup18] by letting $\pi =w-1$ and $\Phi_x=u$ before reversing the sign of x and choosing $\epsilon=-1$ and β = 1.

The BK system is invariant under the following scale transformation. If ${\mathop{{w}\mathopen{\left(\mathord{{x,t}}\right)}}}={\mathop{{w_0}\mathopen{\left(\mathord{{x,t}}\right)}}}$, ${\mathop{{u}\mathopen{\left(\mathord{{x,t}}\right)}}}={\mathop{{u_0}\mathopen{\left(\mathord{{x,t}}\right)}}}$ is any solution of Eq. (2.6) and d is any non-zero real constant, then ${\mathop{{w}\mathopen{\left(\mathord{{x,t}}\right)}}}=d^2{\mathop{{w_0}\mathopen{\left(\mathord{{dx,d^2t}}\right)}}}$, ${\mathop{{u}\mathopen{\left(\mathord{{x,t}}\right)}}}=d{\mathop{{u_0}\mathopen{\left(\mathord{{dx,d^2t}}\right)}}}$ is also a solution. In particular, by choosing $d=-1$ it follows that the BK system is symmetric with respect to the wave direction.

2.2. The soliton solution of the BK equation

The soliton solution of the normalized BK system in Eq. (2.6) can be written as

(2.7)\begin{gather} \left\{ \begin{aligned} {\mathop{{w}\mathopen{\left(\mathord{{x,t}}\right)}}}&= 1 +2\kappa^2{\mathop{{\operatorname{sech}^2}\mathopen{\left(\mathord{{{\theta}\left({x,t}\right)-\beta}}\right)}}} {}+2\kappa^2{\mathop{{\operatorname{sech}^2}\mathopen{\left(\mathord{{{\theta}\left({x,t}\right)+\beta}}\right)}}} \\ {\mathop{{u}\mathopen{\left(\mathord{{x,t}}\right)}}}&=\frac{4\kappa^2{\mathop{{\operatorname{sign}}\mathopen{\left(\mathord{{c}}\right)}}}}{{\mathop{{\cosh^2}\mathopen{\left(\mathord{{{\theta}\left({x,t}\right)}}\right)}}}+\tfrac 12 (\lvert{c}\rvert-1)} \end{aligned} , \right. \end{gather}

where

(2.8)\begin{align} {\mathop{{\theta}\mathopen{\left(\mathord{{x,t}}\right)}}}&:= \kappa (x-ct) - \varphi, \end{align}

(2.9)\begin{align} \kappa&=\tfrac 12 \sqrt{c^2-1}, \end{align}

(2.10)\begin{align} \beta&:=\tfrac 12\ln{(\lvert{c}\rvert-2\kappa)}, \end{align}

(2.11)\begin{align} c&\in\left\{c\in\mathbb{R}\mid\lvert{c}\rvert\geq 1\right\}, \end{align}

(2.12)\begin{align} \varphi&\in\mathbb{R} \end{align}

and ${\mathop{{\operatorname{sign}}\mathopen{\left(\mathord{{c}}\right)}}}=-1$ if c < 0, ${\mathop{{\operatorname{sign}}\mathopen{\left(\mathord{{0}}\right)}}}=0$ and ${\mathop{{\operatorname{sign}}\mathopen{\left(\mathord{{c}}\right)}}}=+1$ if c > 0. The free parameters are the phase offset φ and the normalized celerity c. Here, the celerity is positive for waves that move to the right and negative for waves that move to the left. We remark that the representation of the BK system used in [Reference Kaup18] defines x with a reversed sign compared to the physical BK system and by that switches left and right moving waves. For $\lvert{c}\rvert \gt 2$, Eq. (2.7) becomes a double crested wave [Reference Zhang and Li47, §II], which might be surprising in view of physical considerations. However, if $\lvert{c}\rvert \gt 2$, the amplitude of the free surface elevation is greater than twice the still water level ( $w \gt 3\Leftrightarrow\eta \gt 2h_0$). That amplitude is too large for a shallow water wave and hence the BK system is not an adequate model for such waves in the first place. For shallow water waves $\lvert{c}\rvert \lt 2$, and Eq. (2.7) describes a single crested wave.

2.3. The eigenvalues of the BK-NFT spectrum

In the inverse scattering method, the NFT associated to an integrable equation is derived from the spectrum of the first operator in the corresponding Lax pair. For the BK system, we need to solve the spectral problem [Reference Kaup18, Eq. (2.1)]

(2.13)\begin{gather} {\mathop{{f_{xx}}\mathopen{\left(\mathord{{x,t,c}}\right)}}} = -{\mathop{{Q}\mathopen{\left(\mathord{{x,t,c}}\right)}}}{\mathop{{f}\mathopen{\left(\mathord{{x,t,c}}\right)}}}, \end{gather}

where

(2.14)\begin{gather} {\mathop{{Q}\mathopen{\left(\mathord{{x,t,c}}\right)}}}:= -\frac 14 \left(c - \frac{{\mathop{{u}\mathopen{\left(\mathord{{x,t}}\right)}}}}{2} - \sqrt{{\mathop{{w}\mathopen{\left(\mathord{{x,t}}\right)}}}}\right)\left(c - \frac{{\mathop{{u}\mathopen{\left(\mathord{{x,t}}\right)}}}}{2} + \sqrt{{\mathop{{w}\mathopen{\left(\mathord{{x,t}}\right)}}}}\right) \in\mathbb{R} . \end{gather}

Here we have chosen to use the normalized celerity c as the spectral parameter. However, we will write κ meaning ${\mathop{{\kappa}\mathopen{\left(\mathord{{c}}\right)}}}=\tfrac 12\sqrt{c^2-1}$ whenever it allows us to write equations more compactly. Since the spectrum does not depend on the fixed value of the time t, we omit the dependence on t from here on. The eigenvalues are the c for which the spectral problem admits a square integrable solution f. To localize the eigenvalues, typically scattering theory is utilized. We now outline the procedure for the BK-NFT. We assume that both ${\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}}$ and ${\mathop{{w}\mathopen{\left(\mathord{{x,t}}\right)}}}-1$ satisfy vanishing boundary conditions, Footnote ² to wit

\begin{align*} \int_{-\infty}^\infty \lvert{{\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}}}\rvert(1+\lvert{x}\rvert)\mathord{\textrm{d}} x& \lt \infty, & \int_{-\infty}^\infty \lvert{{\mathop{{w}\mathopen{\left(\mathord{{x}}\right)}}}-1}\rvert(1+\lvert{x}\rvert)\mathord{\textrm{d}} x& \lt \infty. \end{align*}

Then we can define two sets of linearly independent Jost solutions of Eq. (2.13) by the boundary conditions [Reference Kaup18, Eq. (2.9)]

(2.15)\begin{gather} \left\{ \begin{aligned} {\mathop{{\phi}\mathopen{\left(\mathord{{x,c}}\right)}}}&\to{\mathop{{\exp}\mathopen{\left(\mathord{{\kappa x}}\right)}}}\\ {\mathop{{\bar{\phi}}\mathopen{\left(\mathord{{x,c}}\right)}}}&\to{\mathop{{\exp}\mathopen{\left(\mathord{{-\kappa x}}\right)}}} \end{aligned} \right. \quad\text{as }x\to-\infty, \end{gather}

(2.16)\begin{gather} \left\{ \begin{aligned} {\mathop{{\bar{\psi}}\mathopen{\left(\mathord{{x,c}}\right)}}}&\to{\mathop{{\exp}\mathopen{\left(\mathord{{\kappa x}}\right)}}}\\ {\mathop{{\psi}\mathopen{\left(\mathord{{x,c}}\right)}}}&\to{\mathop{{\exp}\mathopen{\left(\mathord{{-\kappa x}}\right)}}} \end{aligned} \right. \quad\text{as }x\to+\infty. \end{gather}

The coefficients that express one set of Jost solutions as a linear combination of the other are known as the scattering parameters. These are implicitly defined by [Reference Kaup18, Eq. (2.11)]

(2.17)\begin{gather} \begin{bmatrix} {\mathop{{\phi}\mathopen{\left(\mathord{{x,c}}\right)}}} & {\mathop{{\bar{\phi}}\mathopen{\left(\mathord{{x,c}}\right)}}} \end{bmatrix} = \begin{bmatrix} {\mathop{{\bar{\psi}}\mathopen{\left(\mathord{{x,c}}\right)}}} & {\mathop{{\psi}\mathopen{\left(\mathord{{x,c}}\right)}}} \end{bmatrix} \underbrace{ \begin{bmatrix} {\mathop{{a}\mathopen{\left(\mathord{{c}}\right)}}} &{\mathop{{\bar{b}}\mathopen{\left(\mathord{{c}}\right)}}} \\ {\mathop{{b}\mathopen{\left(\mathord{{c}}\right)}}} &{\mathop{{\bar{a}}\mathopen{\left(\mathord{{c}}\right)}}} \end{bmatrix} }_{=:{\mathop{{{{S}}}\mathopen{\left(\mathord{{c}}\right)}}}} . \end{gather}

Finally, the eigenvalues of the BK system are the values $c=c_n$ for which ${\mathop{{a}\mathopen{\left(\mathord{{c_n}}\right)}}}=0$, or equivalently ${\mathop{{\phi}\mathopen{\left(\mathord{{x,c_n}}\right)}}}\propto{\mathop{{\psi}\mathopen{\left(\mathord{{x,c_n}}\right)}}}$. This can be seen as follows. We can write any solution of Eq. (2.13) as a linear combination of any set of Jost solutions: ${\mathop{{f}\mathopen{\left(\mathord{{x,c}}\right)}}}=\alpha{\mathop{{\phi}\mathopen{\left(\mathord{{x,c}}\right)}}}+ \beta{\mathop{{\bar{\phi}}\mathopen{\left(\mathord{{x,c}}\right)}}} = \gamma{\mathop{{\bar{\psi}}\mathopen{\left(\mathord{{x,c}}\right)}}} + \delta{\mathop{{\psi}\mathopen{\left(\mathord{{x,c}}\right)}}}$. A solution f is an eigenfunction if and only if it is square-integrable. In light of (2.15)–(2.16), this is equivalent to

(2.18)\begin{align} {\mathop{{f}\mathopen{\left(\mathord{{x,c}}\right)}}}=\alpha{\mathop{{\phi}\mathopen{\left(\mathord{{x,c}}\right)}}} = \delta{\mathop{{\psi}\mathopen{\left(\mathord{{x,c}}\right)}}}, \quad \alpha,\delta \ne 0. \end{align}

Eigenvalues in the interval $(1,\infty)$ each correspond to a soliton that moves to the right, whereas eigenvalues in the interval $(-\infty,-1)$ each correspond to a soliton that moves to the left [Reference Kaup18, p. 406]. The single soliton solution (2.7) leads to exactly one eigenvalue. N-soliton solutions are similarly characterized by N eigenvalues [Reference Zhang and Li47]. The continuous part of the NFT is zero for N-solitons.

3.1. Applicability conditions

The spectral problem of the BK system Eq. (2.13) is not a classic SL equation, because ${\mathop{{Q}\mathopen{\left(\mathord{{x,c}}\right)}}}$ is not an affine function of any suitably chosen spectral parameter. However, we will show that under the assumption $\lvert{{\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}}}\rvert \lt 2$ for all x, a result from generalized SL theory can be applied. This condition is unlikely to be violated while using the BK system to model shallow water waves, since $\lvert{u}\rvert\geq 2$ physically corresponds to a particle velocity of $\upsilon\geq2\sqrt{gh_0}$, so at least twice the still water celerity. The single soliton solution Eq. (2.7) breaks this bound if and only if it is double-crested. As discussed earlier, the double crested soliton lies beyond the regime in which the BK system is an adequate model for shallow water waves. We shall furthermore assume that ${\mathop{{w}\mathopen{\left(\mathord{{x}}\right)}}}$ is bounded, which is certainly true for any real-world wave data. Both conditions can be verified from the data before proceeding.

3.2. Bounds on the eigenvalues

Regardless of the data ${\mathop{{w}\mathopen{\left(\mathord{{x}}\right)}}}$, ${\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}}$, every eigenvalue of the BK system Eq. (2.6) satisfies $\lvert{c_n}\rvert \gt 1$ [Reference Kaup18, p. 406]. For given data, we can furthermore compute upper bounds on $\lvert{c}\rvert$ for the left and right moving part of the discrete spectrum as follows. Recall that eigenfunctions satisfy Eq. (2.18). Without loss of generality, we assume α = 1, so that ${\mathop{{f}\mathopen{\left(\mathord{{x,c}}\right)}}} \gt 0$ and ${\mathop{{f_x}\mathopen{\left(\mathord{{x,c}}\right)}}} \gt 0$ for x sufficiently close to $-\infty$ by (2.15). If now ${\mathop{{Q}\mathopen{\left(\mathord{{x,c}}\right)}}}$ were negative for all x, then by Eq. (2.13) ${\mathop{{f_{xx}}\mathopen{\left(\mathord{{x,c}}\right)}}} \gt 0$ for all x. Therefore the function ${\mathop{{f}\mathopen{\left(\mathord{{x,c}}\right)}}}$ would be monotonously increasing for all x, which is incompatible with square-integrability. Hence c cannot be an eigenvalue if ${\mathop{{Q}\mathopen{\left(\mathord{{x,c}}\right)}}} \lt 0$ for all x. This, together with Eq. (2.14), implies that $c\in(c_\text{lb},-1)$ for the left moving part and $c\in(1,c_\text{ub})$ for the right moving part of the discrete spectrum, where

(3.1)\begin{align} c_\text{lb}&=\inf_x \tfrac 12{\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}} -\sqrt{{\mathop{{w}\mathopen{\left(\mathord{{x}}\right)}}}}, \end{align}

(3.2)\begin{align} c_\text{ub}&=\sup_x \tfrac 12{\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}} +\sqrt{{\mathop{{w}\mathopen{\left(\mathord{{x}}\right)}}}}. \end{align}

These bounds can be computed from the data and reduce the search space from two half lines to two finite (or even empty) intervals.

3.3. Oscillation theory for the BK system

The accounting function is defined as

(3.3)\begin{align} s(c) &:= \text{number of zero-crossings in }\phi(x,c)\,\text{for }-\infty \lt x \lt \infty \end{align}

and satisfies

(3.4)\begin{align} s(c) = \text{number of eigenvalues in }\begin{cases} [1,c), & c \in (1,c_\text{ub}) \\ (c,-1], & c \in (c_\text{lb},-1) \end{cases}. \end{align}

By using bracketing to localize the jump points of s(c), we can find arbitrarily small intervals around each eigenvalue. We now discuss how Eq. (3.4) follows from [Reference Greenberg and Babuška17, Thm. 2.1].

Let us assume that u(x) has compact support, which is an implicit assumption underlying the numerical method anyway. Without loss of generality, we assume that $u(x)=0$ and $w(x)=1$ for $x\notin[0,1]$. Due to Eq. (2.18), any eigenfunction must fulfill the boundary condition for ϕ in (2.15) at the left end of the support, and the boundary condition for ψ in (2.16) at the right end of the support. That is,

(3.5)\begin{align} {\mathop{{f}\mathopen{\left(\mathord{{0,c}}\right)}}} = \alpha, \quad {\mathop{{f_x}\mathopen{\left(\mathord{{0,c}}\right)}}} = \alpha\kappa, \quad {\mathop{{f}\mathopen{\left(\mathord{{1,c}}\right)}}} = \beta, \quad {\mathop{{f_x}\mathopen{\left(\mathord{{1,c}}\right)}}} &= -\beta \kappa. \end{align}

The spectral problem for the BK system can thus be written in the form

(3.6)\begin{align} [{\mathop{{P}\mathopen{\left(\mathord{{x,c}}\right)}}} f_x]_x +{\mathop{{Q}\mathopen{\left(\mathord{{x,c}}\right)}}} f &= 0, \quad x\in[0,1], \end{align}

(3.7)\begin{align} {\mathop{{\alpha_0}\mathopen{\left(\mathord{{c}}\right)}}}{\mathop{{f}\mathopen{\left(\mathord{{0,c}}\right)}}} +{\mathop{{\beta_0}\mathopen{\left(\mathord{{c}}\right)}}}{\mathop{{f_x}\mathopen{\left(\mathord{{0,c}}\right)}}} &= 0, \end{align}

(3.8)\begin{align} {\mathop{{\alpha_1}\mathopen{\left(\mathord{{c}}\right)}}}{\mathop{{f}\mathopen{\left(\mathord{{1,c}}\right)}}} +{\mathop{{\beta_1}\mathopen{\left(\mathord{{c}}\right)}}}{\mathop{{f_x}\mathopen{\left(\mathord{{1,c}}\right)}}} &= 0, \end{align}

where ${\mathop{{P}\mathopen{\left(\mathord{{x,c}}\right)}}}:=1$, ${\mathop{{\alpha_0}\mathopen{\left(\mathord{{c}}\right)}}}:=-\kappa$, ${\mathop{{\beta_0}\mathopen{\left(\mathord{{c}}\right)}}}:=1$, ${\mathop{{\alpha_1}\mathopen{\left(\mathord{{c}}\right)}}}:=\kappa$ and ${\mathop{{\beta_1}\mathopen{\left(\mathord{{c}}\right)}}}:=1$. This kind of spectral problem has been considered in [Reference Greenberg and Babuška17]. We now verify that the relevant assumptions in [Reference Greenberg and Babuška17] are fulfilled for the reformulated BK spectral problem if $c\in (1, c_\text{ub})$.

• • Standard assumptions (S1)–(S5): These are fulfilled as soon as ${\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}}$ and ${\mathop{{w}\mathopen{\left(\mathord{{x}}\right)}}}$ are piecewise continuous. (Note that Remark 1 in [Reference Greenberg and Babuška17] clarifies that S1 can be weakened.)

• • Monotonicity assumptions (M1)-(M4): The applicability conditions in Section 3.1 ensure that $Q(x,c)$ is strictly increasing in c on $(1,c_\text{ub})$. The monotonicity conditions are therefore fulfilled.

• • Limit assumption (L3): We have $\alpha_0(c_\text{ub})\beta_0(c_\text{ub})=-\kappa\le0$ and $\alpha_1(c_\text{ub})\beta_1(c_\text{ub})=\kappa\ge0$. The discussion in Section 3.2 furthermore showed that $Q(x,c_\text{ub})\le 0$ for all $x\in[0,1]$.

• • We have ${\mathop{{\bar{f}}\mathopen{\left(\mathord{{c}}\right)}}} :={\mathop{{\alpha_1}\mathopen{\left(\mathord{{c}}\right)}}}{\mathop{{f}\mathopen{\left(\mathord{{1,c}}\right)}}} +{\mathop{{\beta_1}\mathopen{\left(\mathord{{c}}\right)}}}{\mathop{{f_x}\mathopen{\left(\mathord{{1,c}}\right)}}} = \kappa \beta + 1 (-\beta \kappa) = 0$.

The first case in Eq. (3.4) now follows from Part 3 of Theorem 2.1 in [Reference Greenberg and Babuška17] for $c\in (1, c_\text{ub})$. The case $c\in(c_\text{lb},-1)$ follows similarly using the remark on p. 928 of that paper.

3.4. Numerical computation of the accounting function

Numerically, the accounting function is attractive as it can be computed exactly subject to the numerical discretization in a simple way. The numerically computed accounting function increases by one exactly at the “numerical” eigenvalues and is constant otherwise. This property ensures that we localize all eigenvalues. In contrast, algorithms that are based on finding the zeros of a(c) directly can miss closely spaced eigenvalues. This has been discussed in detail in [Reference Prins and Wahls36] for the spectral problem of the KdV equation. Conveniently, the only change that is needed to adapt the discussion and computation of the accounting function to the BK case is to replace $q_d-\kappa^2$, where $q_d:=q(x_d)$, throughout by ${\mathop{{Q}\mathopen{\left(\mathord{{x_d,c}}\right)}}}$, in particular $\gamma:=\sqrt{{\mathop{{Q}\mathopen{\left(\mathord{{x_d,c}}\right)}}}}$. To see this, note that the first-order representation

(3.9)\begin{align} \begin{bmatrix} \phi(x,c) \\ \phi_x(x,c) \end{bmatrix}_x = \begin{bmatrix} 0 & 1 \\ \kappa^2-q({x}) & 0 \end{bmatrix} \begin{bmatrix} \phi(x,c) \\ \phi_x(x,c) \end{bmatrix} , \end{align}

of the spectral problem for the KdV equation in [Reference Prins and Wahls36, Eq. (2)] turns into a first-order representation of the BK spectral problem in Eq. (2.13) under this transformation.

3.5. Numerical computation of the Jost functions and the scattering matrix

We now discuss how the scattering matrix in Eq. (2.17) and its derivative can be computed numerically. Recall that the zeros of the top-left element of the scattering matrix are exactly the eigenvalues, which is why we can use the scattering matrix and its derivative to perform NR refinements of approximate eigenvalues. For numerical computations it is convenient to rewrite Eq. (2.13) as a system of first-order equations. Furthermore we augment the system to also obtain the derivative of the scattering matrix with respect to c, cf. [Reference Boffetta and Osborne6]. To that end we define an operator

(3.10)\begin{gather} {\mathsf{\widetilde V}}_\text{C}:=\begin{bmatrix}1&\frac{\partial}{\partial x}&\frac{\partial}{\partial c}&\frac{\partial^2}{\partial x\partial c}\end{bmatrix}^\top. \end{gather}

Then Eq. (2.13) can be written as

(3.11)\begin{gather} \frac{\partial}{\partial x}{\mathop{{{{\widetilde f}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}={\mathop{{{{\widetilde A}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}{\mathop{{{{\widetilde f}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}, \end{gather}

where ${\mathop{{{{f}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}:= {\mathsf{V}}_\text{C}{\mathop{{f}\mathopen{\left(\mathord{{x,c}}\right)}}}$ and

(3.12)\begin{gather} {\mathop{{{{\widetilde A}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} := \begin{bmatrix} 0 & 1 & 0 & 0 \\ -{\mathop{{Q}\mathopen{\left(\mathord{{x,c}}\right)}}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ -{\mathop{{Q_c}\mathopen{\left(\mathord{{x,c}}\right)}}} & 0 & -{\mathop{{Q}\mathopen{\left(\mathord{{x,c}}\right)}}} & 0 \end{bmatrix}, \end{gather}

where ${\mathop{{Q_c}\mathopen{\left(\mathord{{x,c}}\right)}}}=\tfrac 14 (2c -{\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}})$. In the representation of Eq. (3.11), the boundary conditions of the Jost solutions can be written as

(3.13)\begin{align} {\mathop{{{{\widetilde\varPhi}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} &:= \begin{bmatrix} {\mathop{{\phi}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\bar\phi}\mathopen{\left(\mathord{{x,c}}\right)}}} & 0 & 0 \\ {\mathop{{\phi_x}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\bar\phi_x}\mathopen{\left(\mathord{{x,c}}\right)}}} & 0 & 0 \\ {\mathop{{\phi_c}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\bar\phi_c}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\phi}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\bar\phi}\mathopen{\left(\mathord{{x,c}}\right)}}} \\ {\mathop{{\phi_{xc}}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\bar\phi_{xc}}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\phi_x}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\bar\phi_x}\mathopen{\left(\mathord{{x,c}}\right)}}} \end{bmatrix} \to {\mathop{{{{\widetilde{T}}}_\text{E}^\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} \text{ as } x\to-\infty, \end{align}

(3.14)\begin{align} {\mathop{{{{\widetilde\varPsi}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} &:= \begin{bmatrix} {\mathop{{\bar\psi}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\psi}\mathopen{\left(\mathord{{x,c}}\right)}}} & 0 & 0 \\ {\mathop{{\bar\psi_x}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\psi_x}\mathopen{\left(\mathord{{x,c}}\right)}}} & 0 & 0 \\ {\mathop{{\bar\psi_c}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\psi_c}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\bar\psi}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\psi}\mathopen{\left(\mathord{{x,c}}\right)}}} \\ {\mathop{{\bar\psi_{xc}}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\psi_{xc}}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\bar\psi_x}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{\psi_x}\mathopen{\left(\mathord{{x,c}}\right)}}} \end{bmatrix} \to {\mathop{{{{\widetilde{T}}}_\text{E}^\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} \text{ as } x\to+\infty, \end{align}

where

(3.15)\begin{equation} \begin{aligned} {\mathop{{{{\widetilde{T}}}_\text{E}^\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} := \begin{bmatrix} {\mathop{{{{T}}_\text{S}^\text{C}}\mathopen{\left(\mathord{{c}}\right)}}} & \begin{bmatrix}0&0\\0&0\end{bmatrix} \\ \frac{\mathord{\textrm{d}}}{\mathord{\textrm{d}} c}{\mathop{{{{T}}_\text{S}^\text{C}}\mathopen{\left(\mathord{{c}}\right)}}} &{\mathop{{{{T}}_\text{S}^\text{C}}\mathopen{\left(\mathord{{c}}\right)}}} \end{bmatrix} \begin{bmatrix} {\mathop{{{{T}}_\text{E}^\text{S}}\mathopen{\left(\mathord{{x,c}}\right)}}} & \begin{bmatrix}0&0\\0&0\end{bmatrix}\\ \frac{\partial}{\partial c}{\mathop{{{{T}}_\text{E}^\text{S}}\mathopen{\left(\mathord{{x,c}}\right)}}} &{\mathop{{{{T}}_\text{E}^\text{S}}\mathopen{\left(\mathord{{x,c}}\right)}}} \end{bmatrix}, \end{aligned} \end{equation}

\begin{equation*} \begin{aligned} {\mathop{{{{T}}_\text{S}^\text{C}}\mathopen{\left(\mathord{{c}}\right)}}} &= \begin{bmatrix} 1&1\\\kappa&-\kappa \end{bmatrix} , & {\mathop{{{{T}}_\text{E}^\text{S}}\mathopen{\left(\mathord{{x,c}}\right)}}} &= \begin{bmatrix} {\mathop{{\exp}\mathopen{\left(\mathord{{\kappa x}}\right)}}} & 0 \\ 0 &{\mathop{{\exp}\mathopen{\left(\mathord{{-\kappa x}}\right)}}} \end{bmatrix} , \\ \frac{\mathord{\textrm{d}}}{\mathord{\textrm{d}} c}{\mathop{{{{T}}_\text{S}^\text{C}}\mathopen{\left(\mathord{{c}}\right)}}} &= \frac{c}{4\kappa} \begin{bmatrix} 0&0\\1&-1 \end{bmatrix} , & \frac{\partial}{\partial c}{\mathop{{{{T}}_\text{E}^\text{S}}\mathopen{\left(\mathord{{x,c}}\right)}}} &= \frac{cx}{4\kappa} \begin{bmatrix} 1&0\\0&-1 \end{bmatrix} {\mathop{{{{T}}_\text{E}^\text{S}}\mathopen{\left(\mathord{{x,c}}\right)}}} . \end{aligned} \end{equation*}

Since ${\mathop{{{{{\widetilde{\Psi}}}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}$ is invertible for $\lvert{c}\rvert \gt 1$ by the linear independence of ${\mathop{{\bar\psi}\mathopen{\left(\mathord{{x,c}}\right)}}}$ and ${\mathop{{\psi}\mathopen{\left(\mathord{{x,c}}\right)}}}$, once ${\mathop{{{{{\widetilde{\Phi}}}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}$ and ${\mathop{{{{{\widetilde{\Psi}}}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}$ are known for the same position x, the scattering matrix and its derivative with respect to c can be computed as

(3.16)\begin{equation} \begin{aligned} \begin{bmatrix} {\mathop{{{{S}}}\mathopen{\left(\mathord{{c}}\right)}}} & \begin{bmatrix}0&0\\0&0\end{bmatrix}\\ \tfrac{\mathord{\textrm{d}}}{\mathord{\textrm{d}} c}{\mathop{{{{S}}}\mathopen{\left(\mathord{{c}}\right)}}}&{\mathop{{{{S}}}\mathopen{\left(\mathord{{c}}\right)}}} \end{bmatrix} \equiv {\mathop{{{{\widetilde\varPsi}}^{-1}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} {\mathop{{{{\widetilde\varPhi}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} {.} \end{aligned} \end{equation}

We remark that computation time can be saved by pre- and post-multiplying the expression above by suitable selection matrices to compute only the required values, e.g.

(3.17)\begin{align} \begin{bmatrix} {\mathop{{a}\mathopen{\left(\mathord{{c}}\right)}}}\\ \frac{\mathord{\textrm{d}}}{\mathord{\textrm{d}} c}{\mathop{{a}\mathopen{\left(\mathord{{c}}\right)}}} \end{bmatrix} &= \left( \begin{bmatrix} 1&0&0&0\\ 0&0&1&0\\ \end{bmatrix} {\mathop{{{{\widetilde\varPsi}}^{-1}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} \right) \left( {\mathop{{{{\widetilde\varPhi}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} \begin{bmatrix} 1&0&0&0 \end{bmatrix}^\top \right), \end{align}

(3.18)\begin{align} {\mathop{{{{S}}}\mathopen{\left(\mathord{{c}}\right)}}} &= \left( \begin{bmatrix} 1&0&0&0\\ 0&1&0&0 \end{bmatrix} {\mathop{{{{\widetilde\varPsi}}^{-1}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} \right) \left( {\mathop{{{{\widetilde\varPhi}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}} \begin{bmatrix} 1&0&0&0\\ 0&1&0&0 \end{bmatrix}^\top \right). \end{align}

At this point it remains to propagate ${\mathop{{{{\widetilde{\Phi}}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}$ and ${\mathop{{{{\widetilde{\Psi}}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}$ to the same position x. Firstly we truncate the data to a window $x\in[X_-,X_+]$. That is, we assume ${\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}}=0$ and ${\mathop{{w}\mathopen{\left(\mathord{{x}}\right)}}}=1$ for all $x\notin[X_-,X_+]$, such that ${\mathop{{{{{\widetilde{\Phi}}}}_\text{C}}\mathopen{\left(\mathord{{X_-,c}}\right)}}}={\mathop{{{{\widetilde{T}}}_\text{E}^\text{C}}\mathopen{\left(\mathord{{X_-,c}}\right)}}}$ and ${\mathop{{{{{\widetilde{\Psi}}}}_\text{C}}\mathopen{\left(\mathord{{X_+,c}}\right)}}}={\mathop{{{{\widetilde{T}}}_\text{E}^\text{C}}\mathopen{\left(\mathord{{X_+,c}}\right)}}}$. Secondly we solve Eq. (3.11) by the exponential midpoint rule. That is, we approximate ${\mathop{{u}\mathopen{\left(\mathord{{x}}\right)}}}$ and ${\mathop{{w}\mathopen{\left(\mathord{{x}}\right)}}}$ by piecewise constant functions. Let D be the number of piecewise constant steps, then the stepsize is $\varepsilon:= (X_+-X_-)/ D$ and the midpoints are $x_d:= X_-+(d-\tfrac 12)\varepsilon$ for $d=1,2,\ldots,D$. In each step Eq. (3.11) can be solved exactly:

(3.19)\begin{gather} {\mathop{{{{\widetilde{f}}}_\text{C}}\mathopen{\left(\mathord{{x_d+\tfrac{\varepsilon}{2},c}}\right)}}} = {\mathop{{{{\widetilde{H}}}_\text{C}}\mathopen{\left(\mathord{{x_d,c}}\right)}}} {\mathop{{{{\widetilde{f}}}_\text{C}}\mathopen{\left(\mathord{{x_d-\tfrac{\varepsilon}{2},c}}\right)}}}, \end{gather}

where

(3.20)\begin{align} {\mathop{{{{\widetilde{H}}}_\text{C}}\mathopen{\left(\mathord{{x_d,c}}\right)}}} &:= \begin{bmatrix} {\mathop{{{{H}}_\text{C}}\mathopen{\left(\mathord{{x_d,c}}\right)}}} &\begin{bmatrix}0&0\\0&0\end{bmatrix} \cr \frac{\partial}{\partial c}{\mathop{{{{H}}_\text{C}}\mathopen{\left(\mathord{{x_d,c}}\right)}}} &{\mathop{{{{H}}_\text{C}}\mathopen{\left(\mathord{{x_d,c}}\right)}}} \end{bmatrix}, \end{align}

(3.21)\begin{align} {\mathop{{{{H}}_\text{C}}\mathopen{\left(\mathord{{x_d,c}}\right)}}} &= {\mathop{{\exp}\mathopen{\left(\mathord{{\begin{bmatrix}0&1\\- {Q}\left({x_d,c}\right) & 0\end{bmatrix}\varepsilon}}\right)}}} = \begin{bmatrix} {\mathop{{\cos}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} & \varepsilon{\mathop{{\operatorname{sinc}}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} \\ -\gamma{\mathop{{\sin}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} &{\mathop{{\cos}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} \end{bmatrix}, \end{align}

(3.22)\begin{align} \frac{\partial}{\partial c}{\mathop{{{{H}}_\text{C}}\mathopen{\left(\mathord{{x_d,c}}\right)}}} &= \bigl(2c-{\mathop{{u}\mathopen{\left(\mathord{{x_d}}\right)}}}\bigr) \frac{\varepsilon}{8} \begin{bmatrix} \varepsilon{\mathop{{\operatorname{sinc}}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} & \varepsilon^2{\mathop{{h}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} \\ {\mathop{{\operatorname{sinc}}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} +{\mathop{{\cos}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} & \varepsilon{\mathop{{\operatorname{sinc}}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}} \end{bmatrix}, \end{align}

where $\gamma:=\sqrt{{\mathop{{Q}\mathopen{\left(\mathord{{x_d,c}}\right)}}}}$, and

(3.23)\begin{align} {\mathop{{\operatorname{sinc}}\mathopen{\left(\mathord{{\vartheta}}\right)}}}&:= \begin{cases} \mathord{{\mathop{{\sin}\mathopen{\left(\mathord{{\vartheta}}\right)}}}}/ \vartheta & \vartheta\neq 0\\ 1 & \vartheta=0 \end{cases}, \end{align}

(3.24)\begin{align} {\mathop{{h}\mathopen{\left(\mathord{{\vartheta}}\right)}}}&:=\begin{cases} ({\mathop{{\operatorname{sinc}}\mathopen{\left(\mathord{{\vartheta}}\right)}}}-{\mathop{{\cos}\mathopen{\left(\mathord{{\vartheta}}\right)}}})/{\vartheta^2}& \vartheta\neq 0,\\ {1}/{3}&\vartheta=0\ \end{cases} \end{align}

(3.25)\begin{align} &=\sum_{m=0}^\infty \frac{(-1)^m\,\vartheta^{2m}}{(2m+1)!\,(2m+3)} = \frac 13 - \frac{\vartheta^2}{30} + \frac{\vartheta^4}{840} - \frac{\vartheta^6}{45360} + \ldots. \end{align}

We will evaluate ${\mathop{{h}\mathopen{\left(\mathord{{\gamma\varepsilon}}\right)}}}$ for $\lvert{\gamma\varepsilon}\rvert \lt 1$ by a truncated Taylor series expansion of order 17 (i.e. 9 non-zero terms) to prevent catastrophic cancellation near $\gamma\varepsilon=0$. Finally we can propagate ${\mathop{{{{\widetilde\varPhi}}_\text{C}}\mathopen{\left(\mathord{{x,c}}\right)}}}$ from $x=X_-$ to $x=X_+$ as

(3.26)\begin{gather} {\mathop{{{{\widetilde\varPhi}}_\text{C}}\mathopen{\left(\mathord{{X_+,c}}\right)}}} = {\mathop{{{{\widetilde{H}}}_\text{C}}\mathopen{\left(\mathord{{x_D,c}}\right)}}} {\mathop{{{{\widetilde{H}}}_\text{C}}\mathopen{\left(\mathord{{x_{D-1},c}}\right)}}} \cdots {\mathop{{{{\widetilde{H}}}_\text{C}}\mathopen{\left(\mathord{{x_2,c}}\right)}}} {\mathop{{{{\widetilde{H}}}_\text{C}}\mathopen{\left(\mathord{{x_1,c}}\right)}}} {\mathop{{{{\widetilde\varPhi}}_\text{C}}\mathopen{\left(\mathord{{X_-,c}}\right)}}}. \end{gather}

5.1. Multisoliton solutions of the BK system

The multi-soliton solutions are defined in terms of real parameters $\{\lambda_j\}_{j=1}^n$ satisfying $|\lambda_j| \gt 1$ and nonzero real parameters $\{b_j\}_{j=1}^n$ such that λ_j and b_j have the same sign. These parameters can be chosen freely subject to the given constraints. The eigenvalues are determined by λ_j according to the relationship $c_j=-\tfrac{1}{2} (\lambda_j + 1/\lambda_j)$ while b_j determines the phase of the jth soliton (and are therefore related to $b(c_j)$). We can then form the matrix-valued function ${{M}}(x,t)$ with entries

(5.1)\begin{align} M_{ij} (x,t) = \delta_{ij} + \frac{b_i \lambda_i}{\lambda_i \lambda_j -1} \exp \left(-(\lambda_i-\lambda_i^{-1})x-\frac{1}{2}(\lambda_i^2-\lambda_i^{-2})t \right), \end{align}

and the vector-valued function $\mathbf{y}(x,t)$ with entries

(5.2)\begin{align} y_j (x,t)= b_j \exp \left(-(\lambda_i-\lambda_i^{-1})x-\frac{1}{2}(\lambda_i^2-\lambda_i^{-2})t \right). \end{align}

Using the entries $f_i(x,t)$ in the unique solution $\mathbf{f}(x,t)$ to

(5.3)\begin{align} {{M}}(x,t) \mathbf{f}(x,t) = \mathbf{y} (x,t) , \end{align}

we can produce solutions to the nondimensional version of the BK system from [Reference Nabelek and Zakharov29] using

(5.4)\begin{align} \varphi(x,t) &= - \log \left( 1- \sum_{j=1}^n \frac{f_j(x,t)}{\lambda_j} \right), \end{align}

(5.5)\begin{align} \eta(x,t) &= - \frac{1}{2} \frac{\partial^2 \varphi}{\partial x^2} (x,t) - \sum_{j=1}^n \frac{\partial f_j}{\partial x} (x,t). \end{align}

The function $\varphi(x,t)$ has the interpretation of a velocity potential, so it makes sense to define the velocity $v(x,t) = \frac{\partial \varphi}{\partial x} (x,t)$. The solutions to the version of the Kaup-Broer solutions $w(x,t)$ and $u(x,t)$ to the BK system considered here can the be produced as

(5.6)\begin{align} w(x,t) = \eta\left(\frac{x}{2},\frac{t}{2} \right) + 1, \quad u(x,t) = v\left( \frac{x}{2}, \frac{t}{2} \right). \end{align}

The above formulas can be implemented symbolically, and this provides exact formulas for a class of solutions depending on the parameters λ_j and b_j for which the eigenvalues c_j are known.

5.2. Validation of a 5-soliton solution

We first validate the numerical computation of the eigenvalues on a 5-soliton solution $w(x,t)$, $u(x,t)$ to the BK system produced in Mathematica using the above approach with parameters λ_n and b_n given in Table 1.

Table 1. Parameters and eigenvalues of the 5-soliton solution.

The corresponding eigenvalues c_n for this solution were evaluated up to the accuracy of double precision floating-point numbers and are provided in in Table 1 as well. The 5-soliton synthetic data are shown in Figure 1 at times $t=-50$, t = 0 and t = 50 at 1001 equally spaced grid-points for $x\in[-200,200]$ with both endpoints included. The five eigenvalues detected by the numerical method are shown in Figure 2 as time evolves from $t=-50$ to t = 50 using 1001 grid-points at each time. We see that the computed eigenvalues are consistent and correct with at least 9 correct digits throughout the whole evolution from time $t=-50$ to time t = 50.

Figure 1. An exact 5-soliton solution to the BK system evaluated at three times

Figure 2. The five eigenvalues of the 5-soliton synthetic data are consistently computed as time evolves

The variations occur due to three error sources. First, there is a discretization error because the signal is assumed to be piecewise constant during the derivation of the algorithm. This error can be reduced by decreasing the step size. Second, there is a truncation error because the signal is assumed to be zero outside the computational domain, which is only approximately true for the considered 5-soliton solutions. This error can be reduced by extending the computational domain (and using the same step size, which means more samples). Third, there is a finite precision error due to the limitations of floating-point arithmetic. This error can be reduced by the use of arbitrary precision floating-point numbers. However, the computational cost of the method will increase when any of these errors is reduced.

The theory discussed in the previous sections indicates that the algorithm should be fourth-order accurate, i.e. the error should be proportional to $(\Delta x)^4$, where $\Delta x$ is the grid-spacing. We verify the accuracy order on the synthetic data produced from the exact 5-soliton solution. At time t = 0, the 5-soliton solution was evaluated on 2¹⁵ uniformly space grid-points on $[-200,200)$ with only the left endpoint included. This was then sub-sampled to produce synthetic data with 2ⁿ grid-points for $n=15, 14, ... , 9$ respectively. The absolute value of the errors in the computation of each of the 5 detected eigenvalues for various step sizes $\Delta x = 400/2^n$ on a log-log plot are shown in blue in Figure 3. For each eigenvalue, a line with the same slope as $(\Delta x)^4$ is also plotted, and we see that the rate of convergence is in good agreement with the predicted fourth-order convergence of the method.

Figure 3. The computation of the five eigenvalues converge to machine precision with fourth-order accuracy

5.3. Validation of a 8-soliton solution

We also validated the numerical method for the eigenvalues using an 8-soliton solution with the parameters λ_n and b_n given in Table 2. The 8 eigenvalues c_n for this solution up to the accuracy of double precision floating-point numbers are also provided in Table 2. The 8-soliton synthetic data are shown in Figure 4 at times $t=-20$, t = 0, and t = 20 sampled on 1001 equally space grid-points for $x \in [-200,200]$ with both endpoints. The eight eigenvalues detected by the algorithm as time evolves from $t=-20$ to t = 20 using 1001 grid-points at each time are shown in Figure 5. We see that the computed eigenvalues are correct with at least 8 correct digits throughout the whole evolution from time $t=-20$ to time t = 20. This again indicates the algorithm is accurate and can consistently detect the correct eigenvalues with good accuracy throughout the time evolution from time $t=-20$ to time t = 20. Notice that the errors are most severe near the earliest time for the 8-soliton solution, which is likely due to the edge of a soliton starting to pass out of the interval $[-200,200]$. We did not conduct the fourth-order convergence tests, because when evaluating the 8-soliton solution in Mathematica there appeared to be random noise on the order of 10⁻⁸ introduced in the evaluation of the exact 8-soliton solution, likely due to cut-off error in the calculation of the solution to the system of linear equations. This made it so the smallest possible error in the numerical calculation due to the noise was achieved before the trend could be observed clearly on the log-log plot. However, we note that the error obtained using only a grid spacing of 0.4 has many correct digits, and is likely accurate enough for any practical application in which possible measurement errors could be assumed to be much larger than the error introduced in evaluation of the 8-soliton solution.

Figure 4. An exact 8-soliton solution to the BK system evaluated at three times

Figure 5. The eight eigenvalues of the 8-soliton synthetic data are consistently computed as time evolves

Table 2. Parameters and eigenvalues of the 8-soliton solution.