2. Pöppe’s bracket for semi-additive operators

Our method extends that of Pöppe [Reference Pöppe59], who realized that solutions of $KP$ are given by tau functions, namely Fredholm determinants of suitable families of integral operators. Pöppe considers semi-additive operators, which one can define as families of integral operators $\Phi_{(t)}$ on $L^2((0, \infty ), {\mathbb{C}})$ that depend upon complex parameters $(t)=(t_j)_{j=1}^\infty$ so that $\Phi$ has kernel $\phi (x+t_1, y +t_1; t_2, t_3, \dots )$ and

In this section, we choose $t_1=t \gt 0$, fix and suppress $t_2, t_3, \dots $ and regard $(x,y)$ as the primary variables.

Pöppe [Reference Pöppe59] used Fredholm determinants of semi-additive operators to solve the $KP$ equation and additive operators to solve the $KdV$ equation. There were steps towards a universal tau function theorem that would generate solutions to PDE from the tau functions of suitable integral operators. The $\tau$ function satisfies differential equations relating to (2.3). In this section, we make the basic results precise.

Pöppe’s calculations involve certain operator formulas which reduce integral operators to finite-rank operators, and can conveniently be expressed in terminology from [Reference Cuntz and Quillen21, Section 5] or [Reference Bezrukavnikov and Kaledin2, Lemma 3.6] used to describe the Fedosov product. The origins of this idea in geometric quantization date back to Lichnérowicz. Let $\mathscr{L}$ be a unital complex algebra with ideal $\mathscr{F}$; introduce the algebra $\mathscr{M}$ with ideal $\mathscr{J}$ by

so that $\mathscr{J}^2=0$. Let $\pi : \mathscr{M}\rightarrow\mathscr{M}/\mathscr{J}$ be the canonical homomorphism.

Let $(\Omega^{ev}\mathscr{L} , \circ)$ be the space of noncommutative differential forms on $\mathscr{L}$ that have even order, which forms an algebra for Fedosov’s product $a\circ b=ab-dadb$; for a detailed discussion, see [Reference Cuntz and Quillen21]. The following statements should be self-explanatory, with indexing to match [Reference Cuntz and Quillen22, (15)].

Proposition 2.2(i) applies in Theorem 2.4(iii). Clearly, $\partial :\mathscr{L}\rightarrow\mathscr{L}$ is a derivation if and only if $\rho$ is an algebra homomorphism, as in Theorem 2.4(i). Also, $\partial (\mathscr{L} )\subset \mathscr{F} $ if and only if $\pi\circ\rho =I$; this situation may be compared with Theorem 2.4(iv).

The following includes the fundamental identities of Pöppe in terms of translation semigroups; see [Reference Pöppe60, p. 622].

Suppose that $\psi \in L^2 ((0, \infty ); H')$ and $\xi \in L^2 ((0, \infty );H)$, and $\phi (z, \zeta )=\psi (z)\xi (\zeta )$. Then with $R_x=\int_x^\infty \xi (z)\psi (z)dz$, we have a family of bounded linear operators $I+R_x$ $(x \gt 0)$ on $H$ which are invertible for all $x \gt x_0$ for some $x_0 \gt 0$. Then we can introduce a kernel

In previous papers [Reference Blower6], [Reference Blower and Newsham13], we have developed Pöppe’s approach, starting from linear systems and forming differential rings. In Definition 4.9 and Proposition 7.1, we introduce another bracket operation which facilitates calculation of $[K]_{x,x}$ and its derivatives, so we can compute $\tau$. Our results on differential rings are equivalent in some cases to Pöppe’s results on semi-additive operators, although we regard our approach as more natural, more closely aligned with algebraic formalism that is used in differential Galois theory and elsewhere, and less dependent on isolated ingenious identities.

By [Reference Blower6, Lemma 5.1], the solution $T(x,y)$ to (2.3) satisfies

which is a wave equation, albeit with a potential $-2{\frac{d}{dx}}T(x,x)$ that is typically a matrix, not necessarily self-adjoint. In the next Section 3, we develop a functional calculus for this context.

3. Spectra of linear systems

Let $q\in C_b([0, \infty );{\mathbb R} )$ with $\int_0^\infty x\vert q(x)\vert dx \lt \infty$ and consider the differential operator $L=-{\frac{d^2}{dx^2}}+q$, which is self-adjoint for suitable boundary conditions. Hence, for $f\in L^2(0, \infty );{\mathbb C} )$, we can define $\cos (t\sqrt{L})$ so that $u(t,x)=\cos (t\sqrt{L})f(x)$ satisfies the wave equation

Gelfand and Levitan used this as the foundation of their spectral theory of second-order differential operators. In the current paper, we make a modest extension of their theory to deal with potential $q$ that are not necessarily real, so $L$ is not necessarily self-adjoint, but the wave equation is still useful. In Example 3.7, we consider examples of differential equations which in Section 4 we analyse in terms of linear systems.

We review some notions of spectral theory.

Consider the operator

and suppose that $\Delta$ generates a strongly continuous semigroup $e^{t\Delta}$ such that $\Vert e^{t\Delta}\Vert \leq Me^{\omega_0t}$ for all $t\geq 0$. Then $s(\Delta )\leq \omega_0$. The topological boundary of $\sigma (\Delta )$ and the approximate point spectrum of $\Delta$ are related by $\partial \sigma (\Delta )\subseteq \sigma_{ap} (\Delta )$, and $e^{t\sigma_{ap}(\Delta )}\subseteq \sigma_{ap} (e^{t\Delta })$ for all $t\geq 0$. We defer discussion of the point spectrum until Proposition 3.6.

Let $\mathscr{R}_0$ the subalgebra of $\mathscr{L} (H)$ that generated by the set of $r(\Delta )$, where $r$ is a complex proper rational function which is holomorphic on $\sigma (\Delta )$, and let $\mathscr{R}$ be the norm closure of $\mathscr{R}_0$ in $\mathscr{L} (H)$. Given a complex and commutative Banach algebra $\mathscr{A}$ and a homomorphism $\theta:\mathscr{A}\rightarrow \mathscr{R}$, there is an induced map $\mathop{\mathrm{spec}} m(\mathscr{R})\rightarrow \mathop{\mathrm{spec}} m(\cal A)$ between the corresponding maximal ideal spaces given by $\varphi\mapsto \varphi\circ \theta$ for every multiplicative linear functional $\varphi :\mathscr{R}\rightarrow{\mathbb C}$ that corresponds to a maximal ideal of $\mathscr{R}$. A useful choice of $\mathscr{A}$ is obtained from cosine families.

We have two cases to consider for $\Delta$ and $\omega_0$ as in (3.2). First, suppose $\omega_0=0$, so the cosine family $(\cos (t\sqrt{L_1})$ is uniformly bounded on $L^2$. Then there exist an invertible $U\in {\mathscr{L}}(L^2, H)$ and a self-adjoint and non-negative $K\in \mathscr{L}(H)$ such that $\cos (t\sqrt{L_1})=U\cos (t\sqrt {K})U^{-1}$, and in particular $\sigma (L)=\sigma (K)\subseteq [0, \infty )$.

Now consider $\omega_0 \gt 0$. For $\omega \gt \omega_0 \gt 0$, we introduce the horizontal strip $\mathscr{S}_{\omega}=\{\lambda\in {\mathbb C} :\vert \mathop{{\rm Im}} \lambda \vert \lt \omega\}$ with closure $cl(\mathscr{S}_\omega )$.

Let $L_0$ be the self-adjoint differential operator which is densely defined in $H$ with $L_0\geq 0$. Then we consider the Cauchy problem for the symmetric hyperbolic system

with initial condition

Then there exists a unique solution

To determine the eigenvalues of $I-{\frac{d^2}{dx^2}}+U^\dagger$, it is convenient to work with its inverse operator. In Proposition 3.6, we obtain a criterion for eigenvalues involving Pincus’s principal function [Reference Carey and Pincus19], which requires the following Lemma.

From Lemma 3.5, we deduce that $G_U=X+iY$ where $X,Y$ are bounded and self-adjoint operators in $L^2$ such that $2i[X,Y]=[G_U^\dagger ,G_U]$ is trace class; one says that $(X,Y)$ are an almost commuting pair of self-adjoint operators. With nonzero $\lambda\in \rho (X)$ and $\ell \in \rho (Y)$, we have $W=\log (I-X/\lambda )$ and $Z=\log (I-Y/\ell )$ such that $[W,Z]\in \mathscr{L}^1(H)$; then $\det (e^We^Ze^{-W}e^{-Z})=\exp(\mathop{\mathrm{trace}} [W,Z])$. Pincus has shown that there exists a compactly supported $g\in L^1({\mathbb R}\times {\mathbb R} ; {\mathbb R} )$ such that

We introduce the numerical range of $G_U^\dagger$ as

Since $U-U^\dagger$ is typically non-zero and skew, there is no reason to presume that $\mathscr{W}(G_U^\dagger )$ is contained in ${\mathbb R}$.

Let $\mathscr{R}_0$ be the subalgebra of $\mathscr{L} (H)$ that is generated by the set of $r(\Delta^\dagger)$, where $\Delta^\dagger =I-{\frac{d^2}{dx^2}}+U^\dagger$ and $r$ is a complex and proper rational function which is holomorphic on $\sigma (\Delta )$; then let $\mathscr{R}$ be the norm closure of $\mathscr{R}_0$ in $\mathscr{L} (H)$. Also, if $\lambda $ is a non-zero eigenvalue of $G_U^\dagger$, then $\lambda^{-1}$ is an eigenvalue of $\Delta^\dagger$, so $\lambda^{-1}\in \mathop{\mathrm{spec}} m(\mathscr{R} )$.

Let $\phi : (0, \infty )\rightarrow M_{n\times n}({\mathbb C} )$. We say that $\phi$ is of Floquet type if there exist $\zeta \in {\mathbb C}$ with $\mathop{{\rm Re}}\zeta \gt 0$ and $F: {\mathbb R} \rightarrow M_{n\times n}({\mathbb C} ) $ which is once continuously differentiable and $2\pi$-periodic such that $\phi (t)=e^{-\zeta t}F(t)$.

4. Zakharov–Shabat system

The Zakharov–Shabat system was introduced in [Reference Zakharov and Shabat70], and expressed in terms of operators in [Reference Bauhardt and Pöppe1]. The following definition is suggested by [Reference Jacob, Partington and Pott39].

The hypothesis for Lemma 4.2(iii) is somewhat coarse. In Proposition 4.5(iv) and Proposition 6.1, we only need $Q_x$ and $M_x$ to be Hilbert–Schmidt, and in Section 9, we give examples where this occurs under milder hypotheses. In the abstract, we express our results in terms of Hankels, instead of the less familiar $R_x$.

Now let $H_0={\mathbb C}$. From the linear system

with input and output space $\mathbb{C}^2$ and state space $H\oplus H$, we introduce $\phi (x)=Ce^{-xA}B$ and

so

We note that $M_x, Q_x\geq 0$ as operators, so $\sigma (M_xQ_x)\setminus \{0\}=\sigma (M_x^{1/2}Q_xM_x^{1/2})\setminus \{0\}\subset (0, \infty )$, hence, $I-\lambda M_xQ_x$ is invertible for all $\lambda\in {\mathbb C} \setminus (0, \infty )$

This bracket operation is related to that of Pöppe [Reference Pöppe59, p 622] and Definition 2.3 (ii), and our identities that correspond to his Theorem 3.1 are summarized in the following Proposition.

Now we let $T(x,y)=-\hat Ce^{-x\hat A}(I+\hat R_x)^{-1}e^{-y\hat A}\hat B$, or more explicitly

There exists $x_0$ such that $\Vert Q_xM_x\Vert_{\mathscr{L}(H)} \lt 1$ for all $x \gt x_0$ since $Q_x, M_x\rightarrow 0$ as $x\rightarrow\infty $, so $T$ is well defined.

Note that we can write

where $X(x,x)=\overline{V(x,x)}$. Within $M_{2\times 2} (\mathbb{C})$ we note two real-linear subspaces that have real dimension four, namely the quaternions ${\mathbb H}$ and the $2\times 2$ complex Hermitian matrices. These arise as follows.

5. Spectra of Hankel operators associated with ZS

Under the hypotheses of (ii) and (iii) of Proposition 4.6, we consider

By the spectral theorem [Reference Birman and Solomjak3, p 177], a bounded and self-adjoint operator $\Gamma$ on a separable Hilbert space $H$ has a resolution $H=\int^\oplus H(t)\mu (dt)$ where $\mu$ is a positive Radon measure on ${\mathbb R}$ called the scalar spectral measure or maximal spectral type and $\Gamma h(t)=th(t)$ where $h(t)\in H(t)$. The spectral multiplicity is $\nu :{\mathbb R}\rightarrow {\mathbb N}\cup \{\infty \}$ with $\nu (t)=\dim H(t)$. In particular, let $\mathop{\mathrm{null}} (\Gamma )=\{\xi\in H:\Gamma \xi =0\}$.

In the case of a nonnegative and self-adjoint compact Hankel operator, nonzero eigenvalues are all simple. This is in contrast to Proposition 5.1, as we have not proposed a bound on $\nu (t)$. The decay rate of singular numbers is reflected in the asymptotics of the Fredholm determinant, as follows.

6. Integral equations relating to KP

We start by introducing families of linear systems and related operators, and obtain a determinant formula which we then express in more classical terms. Let $(-A, B(\zeta ),C)$ be a one-parameter family of continuous-time linear systems with state space $H$ and input and output space $H_1$. Let $\Phi (z,\zeta )=Ce^{-zA}B(\zeta )$ be the impulse response function and $R_x=\int_x^\infty B(\zeta )Ce^{-\zeta A}\, d\zeta$ as in (1.3) these are our basic operator functions, from which we introduce various auxiliary functions. Let $\Theta_x:L^2((0,\infty );H_1)\rightarrow H$ and $\Xi_x: L^2((0,\infty );H_1)\rightarrow H$ be defined by

and

and one computes $\Theta^\dagger_x :H\rightarrow L^2((0, \infty ); H_1)$

Hence

so $\Theta^\dagger_0\Xi_0$ is the integral operator on $L^2((0, \infty );H_1)$ that has kernel $\Phi (z, \zeta )$; while

which is our other basic operator. Next, we let

which determines a Volterra-type integral operator

Squares of Hankel operators are themselves associated with tau functions, and on account of Theorem 2.4(iv) can be simpler to use in computations than Hankel operators. For a specific example relating to Painlevé equation and concentric $KdV$, see [Reference Drazin and Johnson27, p. 172]. In [Reference Blower and Doust9, Section 4], we discussed a differential ring related to the Hastings–McLeod solution of Painlevé II and obtained the following determinant from a Hankel square operator in [Reference Blower and Doust9, Proposition 4.1(ii)] and Example 9.2. The values of $\det (I\pm R_x)$ are equivalent data to the values of $\det (I-R_x^2)$ and $\det ((I+R_x)(I-R_x)^{-1})$ for $R_x\in\mathscr{L}^1$ such that $I-R_x$ is invertible.

Next we provide an explicit expression for this solution in classical style. We use spatial coordinates $(z, \zeta )\in \mathbb{R}^2$, variables $x,y\in\mathbb{R}$, a spectral parameter $\kappa\in {\mathbb{R}}\cup i{\mathbb{R}}$ for an operator in $(z, \zeta )$, and a Fourier transform variable $\omega$ which is dual to $y$.

We show how to obtain $\psi (\zeta ,z;\kappa )$ that is a solution of the system

Let $J_n$ be Bessel’s function of the first kind of order $n\in {\mathbb Z}$, defined by

and let $K_0$ be the modified Bessel function of the second kind

such that $K_0(w)=J_0(iw)$.

We return to general $\beta\in \mathbb{R}$, and observe that

as a function of $y$ is the Fourier transform in the $\omega$-variable of

Then the Gelfand–Levitan equation

has Fourier transform

and we can solve this as above. Suppose that $c:(a_1, a_2)\rightarrow \mathbb{C}$ is an integrable function and

Here, we choose the continuous-time linear systems $(-A,B(\zeta ),C)$ with state space $H=L^2((0, \infty )\times (a_1,a_2); \mathbb{C}^{2\times 1})$ and input and output space $H_1=\mathbb{C}$ with operators

In particular, if

then

7. A differential ring for KP

From $\Phi$, we now obtain a solution of the linearized $KP$ equation in the standard form. We consider the semi-additive family of kernels

such that

and

For each real $x$, there is an operator $P_x$ given by

which is a multiplication operator as a function in the $y$ variable and an integral operator in $\zeta$. Consider a family of linear systems $\Sigma_{(y,t)}=(-A, B(y), C(t))$ with impulse response function $\Phi (z, \zeta ; y,t)=C(t)e^{-(z+\zeta )A}B(y)$; then

gives

which satisfies

We now introduce a version of the Pöppe bracket operation that is suited to the $KP$ equation. The starting point is Lyapunov’s identity

For the family of linear systems $(-A, e^{-y\Delta }B,C)$, we have

then we define $F_{z,y}=(I+R_{z,y})^{-1}$. Also, we have

and

Hence, the second-order partial derivatives of $\log\tau$ satisfy

by calculations as in Proposition 7.2. With $K(x,z)=-Ce^{-xA}(I+R_x)^{-1}e^{-zA}B$, we reconcile the brackets $[\cdot ]$ with $\lfloor\cdot \rfloor$, as in

8. Solution of KP

Suppose that there is another parameter $t$ such that

For the $KP$ equation, we adjust the choice of $L_u$ to

to accommodate the extra variable $t$, nevertheless, we have $L_u(I+K)=(I+K)L_0$, as in Lemma 6.6. As in (7.20), we define

and recall that

Let

9. Examples and Remarks

In Propositions 4.5 and 8.1, we use a trace-class $R_x$ that arises as the product of Hilbert–Schmidt operators. This is a more stringent hypothesis than admissibility, as in 4.1, so in this section, we give conditions for various linear systems to produce operators in trace ideals. Let $H=L^2(\mathbb{R}; \mathbb{C})$ and $\mathscr{D}(A)=\{f\in H: vf(v)\in H\}$, let ${\mathscr{D}}(\Delta )=\{f\in H: v^3f(v)\in H\}$ which are themselves Hilbert spaces for the appropriate graph norms; let $b,c\in H\cap L^\infty (\mathbb{R}; \mathbb{C})$; then introduce bounded linear operators

This example is not covered by Propositions 2.2 and 2.3 of [Reference Blower6], so we give a special argument to show that $R_x$ exists.

10. Numerical simulations

We present numerical simulations of solutions to the KP equation. We use four different numerical approaches as follows:

(1) GLM solution using Riemann Rule approximation (GLM-RR): We solve the linear integral GLM equation, with the coefficients given by semi-additive scattering data representing solutions to the linearized $KP$ equation. The solution to the linearized equations (10.2) can be analytically advanced to any time $t \gt 0$ and substituted into the GLM equation, which is then solved. For this method we use the left-hand Riemann Rule to approximate the integral in the GLM equation, which is then solved as a large linear system of equations at that time $t \gt 0$;

(2) GLM solution using Clenshaw–Curtis quadrature (GLM-CC): This is similar to the last method except that we use Clenshaw–Curtis quadrature to approximate the integral in the GLM equation. Clenshaw–Curtis quadrature is based on Chebyshev polynomial approximation and its use here is inspired by the approximation method for computing Fredholm determinants developed by Bornemann [Reference Bornemann15], as we outline next;

(3) Fredholm determinant using Nyström–Clenshaw–Curtis method (Det-CC): The solution to the $KP$ equation is given in terms of the second derivative, with respect to $x$, of the logarithm of the $\tau$-function, which can be expressed as the Fredholm determinant of the scattering data. We use the Nyström–Clenshaw–Curtis method developed by Bornemann [Reference Bornemann15] to approximate Fredholm determinants to very high accuracy;

(4) Direct pseudo-spectral time-stepping approximation (FFT2-exp): This is a direct, exponential split-step, pseudo-spectral method that utilizes the fast Fourier transform (FFT) in both the $x$- and the $y$- directions. It also utilizes the window method outlined by Kao and Kodama [Reference Kao and Kodama40], to deal with non-periodic boundary conditions. This method is outlined in detail in Blower and Malham [Reference Blower and Malham11, App. B].

We provide further details of these numerical approaches presently. Before doing so, let us outline the direct linearization approach we have outlined in Sections 6–8, and connect that to the direct linearization approach given in Blower and Malham [Reference Blower and Malham11], as well as the expression for the solution via the $\tau$-function. The $KP$ equation for the field $g=g(x,y;t)$ in potential form is given by

This corresponds to the KP equation (1.1) in the case $\alpha=-1/4$, $\beta=\pm1$ and $u=2{\partial}_xg$. The linearized form of the $KP$ equation for $p=p(t)$ is given by

The linear integral GLM equation is of the form,

for the solution operator $G$, or equivalently its kernel, $g$. Here, $P$ is the scattering operator associated with the kernel function solution $p$ solving (10.2). We assume that $P$ is a Hilbert–Schmidt valued integral operator on $(-\infty,0]$ with kernel of the form,

This is the semi-additive form first introduced by Pöppe [Reference Pöppe60] with $z,\zeta\in(-\infty,0]$ the primary variables parametrizing the operator $P$, while $x,y\in\mathbb R$ and $t\geqslant0$ are regarded as additional parameters. This form guarantees that the first of the following two constraints on $p$ is automatically satisfied:

These two constraints arise as the first two equations in the $KP$ hierarchy; see for example [Reference Blower and Malham12]. While our semi-additive assumption for $p$ means that (10.5a) is satisfied, we henceforth assume that $p$ satisfies (10.5b) as well. Using the constraints (10.5), the linearized $KP$ equation (10.2), has the alternative formulation,

Note that if $P$ is Hilbert–Schmidt valued on $(-\infty,0]$, then the solution operator $G$ to the GLM equation (10.3) is Hilbert–Schmidt valued on $(-\infty,0]$ as well—see Blower and Malham [Reference Blower and Malham11, Lemma 7]—and the kernel $g=g(z,\zeta;x,y,t)$ of $G$ is square-integrable. Furthermore, $g$ automatically adopts any regularity that $p$ possesses. The solution $g=g(z,\zeta;x,y,t)$ to the GLM equation evaluated at $z=\zeta=0$, i.e. $g=g(0,0;x,y,t)$, satisfies the $KP$ equation (10.1). See, for example, [Reference Blower and Malham11]. Solutions to the $KP$ equation (10.1) can thus be generated by solving the linearized $KP$ equation (10.2), or equivalently here (10.6), and the following linear integral equation for $g=g(0,\zeta;x,y,t)$,

Furthermore, we also know that $g=g(0,0;x,y,t)$ is given by the trace formula,

where $V:=({\mathrm{id}}-P)^{-1}$; see [Reference Blower and Malham11, Cor. 12]. Here, ${\partial}_{\mathfrak{l}}G$ and ${\partial}_{\mathfrak{r}}G$ represent the trace-class operators with the respective kernels ${\partial}_zg$ and ${\partial}_\zeta g$. In other words we have, $g(0,0;x,y,t)=\mathrm{trace}\,\bigl(({\partial}_x P)V\bigr)$ since ${\partial}_xP={\partial}_{\mathfrak{l}}P+{\partial}_{\mathfrak{r}}P$ using (10.5a). By a standard calculation, we thus have,

The quantity $\tau:=\mathrm{det}({\mathrm{id}}-P)$ is known as the $\tau$-function, as in Definition 2.7(ii). Thus, in our numerical methods GLM-RR and GLM-CC outlined above, we solve the linear integral equation (10.7), respectively, using the left-hand Riemann Rule and Clenshaw–Curtis quadrature for the integral on the right-hand side. Utilizing either of these quadrature approximations generates a linear algebraic system of equations for the solution $g$ at the nodal points $\zeta_n$. We give more details on this procedure presently. Furthermore, for our numerical method Det-CC, we use the Nyström–Clenshaw–Curtis method developed by Bornemann [Reference Bornemann15] to approximate the $\tau$-function Fredholm determinant. Lastly, in this context, since $g$ solves the potential form of the $KP$ equation (10.1), the solution to the $KP$ equation itself is ${\partial}_xg(0,0;x,y,t)$.

Before we discuss the implementation of our numerical methods, let us relate the quantities above to those in Sections 6–8. Suppose the operator $R=R(x,y,t)$, see for example (1.3) or (7.5), has the form,

Here, we have assumed that, the operator $A$ simply represents real multiplication by the constant $a$, and the operator $B=B(y,\zeta)$ has the form,

where $b$ is a real constant. Furthermore, suppose we define the operator from $\mathbb C$ to itself, or function, $\hat{K}=\hat{K}(z,\zeta;x,y,t)$ by,

It is then straightforward to verify that $\hat{p}(z+x,\zeta+x;y,t):= C(t)\mathrm{e}^{-a(z+x)-b(\zeta+x)}\hat{B}(y)$ and $\hat{K}=\hat{K}(z,\zeta;x,y,t)$ satisfy the linear integral equation,

If we make the change of variables,

and set,

and

then we see that $p$ and $g$ satisfy the GLM equation $P=G({\mathrm{id}}-P)$, i.e.

This generates (10.7) when we set $z=0$.

Let us now outline in detail the four numerical algorithms we used to compute solutions to the $KP$ equation. For all four numerical methods, we truncate the $(x,y)\in\mathbb R^2$ domain to $[-L_x/2,L_x/2]\times[-L_y/2,L_y/2]$ for sufficiently large domain lengths $L_x \gt 0$ and $L_y \gt 0$. First, we outline the simple solution method, GLM-RR, that solves the linear integral GLM equation, by using the left-hand Riemann Rule to approximate the integral therein. Our goal is to compute $g=g(0,0;x,y,t)$. To achieve this, for given scattering data $p$, and given $x\in[-L_x/2,L_x/2]$, $y\in[-L_y/2,L_y/2]$ and $t\geqslant0$, we numerically solve the GLM equation (10.7) for $g=g(0,\zeta;x,y,t)$ and then set $\zeta=0$. In practice we use $N_x$ nodal points $x_n$ in the truncated $x$-domain and $N_y$ nodal points $y_{n^\prime}$ in the truncated $y$-domain. For each nodal point-pair $(x_n,y_{n^\prime})$ and any given fixed $t\geqslant0$, we numerically solve (10.7) as follows. Note that $\zeta,\xi\in[-L_x/2,0]$. We use $M/2+1$ nodal points $\zeta_m$ and $\xi_{m^\prime}$ for these variables in $[-L_x/2,0]$, with separation $h=L_x/M$. We always take $M$ to be even. Using each of the nodal point $\zeta_m$, we generate the row vector $\widehat{P}$ of length ‘ $M/2+1$’ containing the values $p(x_n,\zeta_m+x_n;y_{n^\prime},t)$ for $m\in\{0,1,\ldots,M/2\}$. Note that $\zeta_{M/2+1}=0$. Furthermore, for each of the nodal point pairs $(\zeta_m,\xi_{m^\prime})$, we generate the matrix $\widehat{Q}$ of size $(M/2+1)\times(M/2+1)$ containing the values $p(\xi_{m^\prime}+x_n,\zeta_m+x_n;y_{n^\prime},t)$ for $m^\prime,m\in(0,1,\ldots,M/2)$. Suppose that $\widehat{G}$ is the row vector of unknown values $g(0,\zeta_m;x_n,y_{n^\prime},t)$ for $m\in(0,1,\ldots,M/2)$. We then solve the linear algebraic system,

for $\widehat{G}$, where $I=I_{M/2+1}$ is the $(M/2+1)\times(M/2+1)$ identity matrix. The matrix product ‘ $h\widehat{G}\widehat{Q}$’ with the multiplicative scaling $h$, naturally implements the left-hand Riemann Rule. We solve the linear system in (10.18) using Gaussian elimination. For a given $t\geqslant0$, the procedure just outlined, is carried out for each $(x_n,y_{n^\prime})$ with $n\in\{0,1,\ldots,N_x\}$ and $n^\prime\in\{0,1.\ldots,N_y\}$. For each $(x_n,y_{n^\prime})$ we extract the final $(M/2+1)$th component of $\widehat{G}$ which represents an approximation for $g(0,0;x_n,y_{n^\prime},t)$. This outlines the GLM-RR method. We compute the solution to the $KP$ equation, namely ${\partial}_xg(0,0;x_n,y_{n^\prime},t)$, by approximating the derivative via a finite difference.

Second, we outline the GLM-CC method, inspired by the method for computing Fredholm determinants developed by Bornemann [Reference Bornemann15]. To begin with, this method follows that of the GLM-RR method, we truncate the $(x,y)\in\mathbb R^2$ domain in precisely the same way and utilize the nodes $x_n$ and $y_{n^\prime}$ as outlined above. The difference comes into play with the choice of the $M/2+1$ nodal points $\zeta_m$ and $\xi_{m^\prime}$ in $[-L_x/2,0]$. Here, we choose these nodal points according to the Clenshaw–Curtis quadrature rule approximation for the integral term in (10.7). Indeed, using the notation,

and then suppressing the implicit $x_n$, $y_{n^\prime}$ and $t$-dependence in $\gamma$ and $\rho$, we approximate,

Here, $\zeta_m$ and $\xi_{m^\prime}$ are chosen to be the Clenshaw–Curtis nodal points and the $w_{m}$ are the Clenshaw–Curtis quadrature weights. Clenshaw–Curtis quadrature is based on the expansion of the integrand using Chebyshev polynomials of the first kind. The nodal points and weights are given explicitly in Bornemann [Reference Bornemann15, p. 909], including a Matlab code for generating them, which in fact, we utilized directly. Bornemann points out that alternatively, Gauss–Legendre quadrature could also be used; however, Clenshaw–Curtis quadrature is more efficient. Thus, for each nodal point pair $(x_n,y_{n^\prime})$ and any given fixed $t\geqslant0$, we construct the row vectors $\widehat{P}$ and $\widehat{G}$ and the matrix $\widehat{Q}$, as above, except now based on the Clenshaw–Curtis nodal points $\zeta_m$ and $\xi_{m^\prime}$, which are not uniformly distributed. Let $W$ denote the diagonal matrix of Clenshaw–Curtis quadrature weights. We then solve the linear algebraic system,

for $\widehat{G}$, again using Gaussian elimination. The Clenshaw–Curtis quadrature approximation is implicit in the matrix product $\widehat{G}W\widehat{Q}$. As above, for a given $t\geqslant0$ and for each $(x_n,y_{n^\prime})$ with $n\in\{0,1,\ldots,N_x\}$ and $n^\prime\in\{0,1.\ldots,N_y\}$, we solve the linear system (10.20) and extract the final $(M/2+1)$th component of $\widehat{G}$ which represents our approximation for $g(0,0;x_n,y_{n^\prime},t)$. This outlines the GLM-CC method. We again use a finite difference approximation to compute ${\partial}_xg(0,0;x_n,y_{n^\prime},t)$.

Third, we now outline the Det-CC solution method based on the method for computing Fredholm determinants developed by Bornemann [Reference Bornemann15]. Again, initially this follows the set-up of the GLM-CC method, we truncate the $(x,y)\in\mathbb R^2$ domain as above and utilize the uniformly distributed nodes $x_n$ and $y_{n^\prime}$. Our goal here is to compute the quantity ‘ $\mathrm{det}({\mathrm{id}}-P)$’ in the formula (10.9). Bornemann [Reference Bornemann15] provides a simple and accurate approximation formula for computing such a determinant based on Clenshaw–Curtis quadrature. Using the same notation to that outlined above for the GLM-CC method, Bornemann [Reference Bornemann15, p. 890, 894] suggests we compute,

where $W^{1/2}$ is the diagonal matrix of entries consisting of the square roots of the Clenshaw–Curtis weights. We compute this approximation for any given $t\geqslant0$ and for each $(x_n,y_{n^\prime})$ with $n\in\{0,1,\ldots,N_x\}$ and $n^\prime\in\{0,1.\ldots,N_y\}$. This outlines the Det-CC method. The approximate solution to the $KP$ equation can be computed by approximating the partial derivative ${\partial}_x^2$ of this determinant approximation using a second-order central difference scheme. A comprehensive error and performance analysis is provided in Bornemann [Reference Bornemann15]. We remark on this in our implementation in Example 10.3 below. Fourth, we now outline the exponential split-step pseudo-spectral algorithm we used to directly integrate the $KP$ equations. For convenience, suppose $A$ denotes the linear $KP$ operator, $A:={\partial}_x^3+3{\partial}_x^{-1}{\partial}_y^2$. Then the exponential split-step we use to integrate the $KP$ equation 10.1 is given by,

where $\mathrm{N}(v)=6{\partial}_x(\mathscr{F}^{-1}v)^2$. Here, $\mathscr{F}$ represents the two-dimensional Fourier transform, and $\hat{u}_\ell$ is the two-dimensional Fourier transform of the approximate solution $u$ to the $KP$ equation at time $t_\ell\in\{0\}\cup\mathbb N$. The quantity $\mathscr{F}(A)$ represents the Fourier transform of the operator $A$. In practice, if $2\pi\mathrm{i}k_x/L_x$ and $2\pi\mathrm{i}k_y/L_y$ are the wavenumbers, respectively, in the $x$- and $y$- directions, we set,

Here, following Klein and Roidot [Reference Klein and Roidot41, p. 3341], we have approximated the Fourier transform of ${\partial}_x^{-1}$ by $1/(2\pi\mathrm{i}k_x/L_x+2\pi\delta)$, where $\delta=2^{-52}$. The initial data is generated by numerically solving the GLM equation, using the GLM-CC method, for the given scattering data $p$ at time $t=0$, as outlined above. We use a pseudo-spectral algorithm due to its simplicity and efficiency, see, for example, Klein and Saut [Reference Klein and Saut42] and Grava, Klein, and Pitton [Reference Grava, Klein and Pitton36]. To deal with the fact that the solutions we compute are not periodic, we use the ‘window method’ employed by Kao and Kodama [Reference Kao and Kodama40]. Precise details on this method can be found in Blower and Malham [Reference Blower and Malham11, App. B]. The FFT2-exp method we have employed thus combines the exponential split-step pseudo-spectral algorithm with the ‘window method’.

Figure 1.

We plot the two-soliton interaction solution outlined in Example 10.3 at $t=0$. The top panel set shows the solution computed by numerically solving the GLM equation (10.7) using Clenshaw–Curtis quadrature, i.e. the GLM-CC method. The bottom panel set shows the solution computed using the $\tau$-function Fredholm determinant, i.e. using the Nyström–Clenshaw–Curtis method Det-CC. The right-hand panels show the corresponding contour plots.

1 Long description

The image A showing a three-dimensional surface plot titled GLM Clenshaw-Curtis method T=0. The x-axis is labeled x with tick labels minus 5, 0, 5, 10, 15. The y-axis is labeled y with tick labels minus 5, 0, 5, 10, 15. The vertical axis is labeled g subscript x with tick labels 0.5, 0, minus 0.5, minus 1, minus 1.5, minus 2, minus 2.5. No unit text is shown for x, y, or g subscript x. The surface is near 0 over much of the domain, then drops sharply for x values near 10 to 15, reaching values below minus 2. Several repeating trough-like bands appear along the steep drop region. The image B showing a contour plot titled GLM Clenshaw-Curtis method, partial subscript x g: T=0. The x-axis is labeled x with tick labels minus 6, minus 4, minus 2, 0, 2, 4, 6, 8, 10, 12, 14. The y-axis is labeled y with tick labels minus 6, minus 4, minus 2, 0, 2, 4, 6, 8, 10, 12, 14. No unit text is shown for x or y. Contours form a narrow, near-vertical high-gradient band for x values around 11 to 13 across most y values. Additional contours approach the band along diagonal paths from lower x values at higher y values and the contour spacing tightens around x about 11 to 12 and y about 6 to 8. The image C showing a three-dimensional surface plot titled Determinant Nystrom method T=0. The x-axis is labeled x with tick labels minus 5, 0, 5, 10, 15. The y-axis is labeled y with tick labels minus 5, 0, 5, 10, 15. The vertical axis is labeled partial subscript x squared log tau with tick labels 0.5, 0, minus 0.5, minus 1, minus 1.5, minus 2, minus 2.5. No unit text is shown for x, y, or partial subscript x squared log tau. The surface is near 0 over much of the domain and drops sharply for x values near 10 to 15, reaching values below minus 2, with repeating trough-like bands along the steep drop region. The image D showing a contour plot titled Determinant Nystrom method, partial subscript x squared log tau: T=0. The x-axis is labeled x with tick labels minus 6, minus 4, minus 2, 0, 2, 4, 6, 8, 10, 12, 14. The y-axis is labeled y with tick labels minus 6, minus 4, minus 2, 0, 2, 4, 6, 8, 10, 12, 14. No unit text is shown for x or y. Contours form a narrow, near-vertical high-gradient band for x values around 11 to 13 across most y values. Diagonal contour paths approach the band from lower x values at higher y values and contour spacing tightens around x about 11 to 12 and y about 6 to 8.

Figure 2.

We plot the two-soliton interaction solution outlined in Example 10.3 at $t=0.25$. The top panel set shows the solution computed by numerically solving the GLM equation (10.7) using the GLM-CC method. The middle panel set shows the solution computed using the Det-CC method. The bottom panel set shows the solution computed by direct numerical integration using the FFT2-exp method. The right-hand panels show the corresponding contour plots.

Figure 2 Long description

The image A showing a 3D surface plot titled “GLM Clenshaw Curtis method t=0.25”. The horizontal axes are labeled x and y. The vertical axis is labeled partial derivative x g. The x-axis shows tick labels from minus 5 to 15. The y-axis shows tick labels from minus 5 to 15. The partial derivative x g axis shows tick labels from minus 0.5 to 0.5. The surface includes a broad, nearly flat region and a steep ridge like structure that bends across the plot, with a deep trough adjacent to the ridge. The image B showing a contour plot titled “GLM Clenshaw Curtis method, partial derivative x g: T=0.25”. The x-axis is labeled x with tick labels from minus 6 to 14. The y-axis is labeled y with tick labels from minus 6 to 14. Contour lines form a narrow, nearly vertical band centered near x around 9 to 10 and several diagonal contour bands that intersect this vertical band. The image C showing a 3D surface plot titled “Determinant Nystrom method t=0.25”. The horizontal axes are labeled x and y. The vertical axis is labeled partial derivative x g. The x-axis shows tick labels from minus 5 to 15. The y-axis shows tick labels from minus 5 to 15. The partial derivative x g axis shows tick labels from minus 0.5 to 0.5. The surface shows a broad, nearly flat region and a steep ridge like structure with an adjacent trough, similar in placement to image A. The image D showing a contour plot titled “Determinant Nystrom method, partial derivative x g: T=0.25”. The x-axis is labeled x with tick labels from minus 6 to 14. The y-axis is labeled y with tick labels from minus 6 to 14. Contour lines form a narrow, nearly vertical band near x around 9 to 10, with multiple diagonal contour bands crossing it. The image E showing a 3D surface plot titled “Direct pseudo spectral (FFT2) method t=0.25”. The horizontal axes are labeled x and y. The vertical axis is labeled partial derivative x g. The x-axis shows tick labels from minus 5 to 15. The y-axis shows tick labels from minus 5 to 15. The partial derivative x g axis shows tick labels from minus 0.5 to 0.5. The surface appears more rippled, with multiple smaller undulations across the flatter region and a steep ridge and trough structure toward the right side. The image F showing a contour plot titled “Direct pseudo spectral (FFT2) method, partial derivative x g: T=0.25”. The x-axis is labeled x with tick labels from minus 6 to 14. The y-axis is labeled y with tick labels from minus 6 to 14. Contour lines include a narrow, nearly vertical band near x around 9 to 10, intersected by diagonal contour bands, with additional irregular contour loops and waviness compared with images B and D.