\documentclass[reqno]{amsart} 

\AtBeginDocument{{\noindent\small 
{\em Electronic Journal of Differential Equations},
Vol. 2003(2003), No. 90, pp. 1--27.\newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu  (login: ftp)}
\thanks{\copyright 2003 Texas State University-San Marcos.} 
\vspace{9mm}}

\begin{document} 

\title[\hfilneg EJDE--2003/90\hfil Weak asymptotic method]
{Weak asymptotic method for the study of propagation
and interaction of infinitely narrow $\delta$-solitons} 

\author[Vladimir G. Danilov \& Georgii A. Omel'yanov\hfil EJDE--2003/90\hfilneg]
{Vladimir G. Danilov \& Georgii A. Omel'yanov}  % in alphabetical order

\address{Vladimir G. Danilov\newline
Moscow Technical University of Communication and Informatics}
\email{danilov@miem.edu.ru}

\address{Georgii A. Omel'yanov\newline
Universidad de Sonora,  Rosales y Blvd. Luis Encinas, Hermosillo, Mexico\newline
Permanent address: Moscow State Institute of Electronics
and Mathematics}
\email{omel@hades.mat.uson.mx \quad omel@miem.edu.ru}
\date{}


\thanks{Submitted April 26, 2002. Published September 5, 2003.}
\thanks{Partially supported  by grant F-41421 from  CONACYT Mexico.}
\subjclass[2000]{35Q53, 35Q51, 35C20}
\keywords{KdV type equations, soliton, soliton interaction}


\begin{abstract}
 We present a new method for studying the interaction of solitons
 for non-integrable Korteweg-de Vries (KdV) type equations with small
 dispersion and test this method for the KdV equation.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}

\section*{Introduction}

We present a method for studying the interaction of solitons
for Korteweg-de Vries (KdV) type equations with small dispersion.
It is well known that the KdV equation
$$
u_t+(u^2)_x+\epsilon^2u_{xxx}=0,
$$
where $\epsilon\to+0$ is a small parameter,
has a soliton solution in the form of a traveling wave
\begin{equation}
u=\frac{A}{2}\cosh^{-2}\Big(\beta\frac{x-Vt}{\epsilon}\Big),\quad
\beta=\Big(\frac{A}{12}\Big)^{1/2},\quad
V=4\beta^2.
\tag {0.1}
\end{equation}
If we consider this solution as a generalized function and
calculate its asymptotic behavior in the weak sense,
then we obtain
\begin{equation}
u\simeq\epsilon\frac{A}{\beta}\delta(x-Vt)\frac1{2}
\int_{\mathbb{R}^1} \cosh^{-2}(\eta)\,d\eta.
\tag {0.2}
\end{equation}
Due to this formula we call the solution (0.1) $\epsilon$--$\delta$ solution.

Similarly, if we consider a KdV type equation in the form
\begin{equation}
u_t+(u^m)_x+\epsilon^2u_{xxx}=0,
\tag {0.3}
\end{equation}
where $m$ is integer, $m>2$,
then we obtain a soliton solution in the form
\begin{equation}
u=A\omega\Big(\beta\frac{x-Vt}{\epsilon}\Big),
\tag {0.4}
\end{equation}
where $A>0$ is a constant,
$$
\beta^2=\gamma^{m-1}A^{m-1},\quad
V=a_mA^{m-1},\quad
\gamma=\Big(\frac{m-1}{m+3}\frac{a_ma_2}{a'_2}\Big)^{1/(m-1)}.
$$
Here and in the following,
we denote
$a_k=\int_{\mathbb{R}^1}\omega^k(\eta)\,d\eta$ for any $k\geq1$,
$a'_2=\int_{\mathbb{R}^1}(d\omega(\eta)/d\eta)^2\,d\eta$,
and $\omega(\eta)$ is the exact solution of the model equation
corresponding to (0.3):
\begin{equation}
\frac{d\omega}{d\eta}+\frac{d\omega^m}{d\eta}
+\frac{d^3\omega}{d\eta^3}=0,
\quad\omega\to0\quad\text{as }\eta\to\pm\infty.
\tag {0.5}
\end{equation}
It is clear that this solution belongs to the Schwartz space
$S(\mathbb{R}^1)$.
Again, in the weak sense, we have
$$
u\simeq \epsilon\frac{A}{\beta}\delta(x-Vt)\int_{\mathbb{R}^1}
\omega(\eta)\,d\eta.
$$


It is well known that the interaction of soliton travelling waves
(0.1) can be described by the famous inverse
scattering transform method \cite{AblS}, which cannot be
used for describing the interaction of waves for $m>3$.
However, one can easily see that, in a weak sense,
the waves (0.1) and (0.4) are similar.
So one can put forward the hypothesis that
there must be a procedure for describing
the wave interaction in both the integrable ($m=2,3$) and
non-integrable ($m>3$) cases.

The goal of the present paper is to confirm this hypothesis
and to propose some general procedure
(the {\it weak asymptotic method\/})
for describing the interaction of nonlinear waves.


Note that, since (0.5) has a solution
rapidly decreasing as $|\eta|\to\infty$,
(0.3) has a special class of solutions \cite{MasOm1,MasOm2}, namely,
distorted infinitely narrow solitons of the form
$$
u=u_0(x,t,\epsilon)+W\Big(\beta\frac{x-v(t)}{\epsilon},t,x,\epsilon\Big)+O(\epsilon^N),
$$
where $N\gg1$, $u_0(x,t,\epsilon)$ and $W(\eta,t,x,\epsilon)$ are smooth
functions, and
$$
W(\eta,t,x,0)=A(t)\omega(\eta).
$$
Such solutions correspond to special initial data
(of the same form as the solution) and are very sensitive
to any general small perturbations \cite{MasOm1,MasOm2,IlKal}.
These perturbations need not necessarily lead to instability,
but they take the solutions out of the class of infinitely
narrow solitons.
In this paper we restrict ourselves to considering propagation
and interaction of nonlinear waves corresponding to solutions
from the class of infinitely narrow solitons, and we do not
consider the problem of their stability \cite{IlKal}.

We briefly consider the weak asymptotic method.
The main notation of the method is the following.

1) By $f_\epsilon(x)=O_{\mathcal{D}'}(\epsilon^\alpha)$ we denote
all generalized functions
(depending on $\epsilon$ as on a parameter)
that for any test function satisfy the relation
$$
\langle f_\epsilon(x),\psi(x)\rangle=O(\epsilon^\alpha),
$$
where $\langle,\rangle$ denotes the action of a generalized function
$f_\epsilon(x)$ on a test function $\psi(x)$
and the estimation on the right-hand side is treated
in the usual sense.

2) Let a function $u_\epsilon(x,t)$ belong to
$C^\infty([0,T]\times \mathbb{R}^1_x)$, $T>0$, for $\epsilon>0$ and
belong to $C(0,T;\mathcal{D}'(\mathbb{R}^1_x))$ uniformly in $\epsilon\geq0$.
We say that $u_\epsilon(x,t)$ is a {\it weak asymptotic solution\/}
($\bmod O_{\mathcal{D}'}(\epsilon^\alpha)$) of the equation
$$
Lu=0
$$
if $Lu_\epsilon=O_{\mathcal{D}'}(\epsilon^\alpha)$,
where, uniformly with respect to $\epsilon$,
the right-hand side is a weakly continuous and
weakly piecewise smooth function of $t$ (for more detail, see
below).

In fact, this means that the difference between the
method of weak asymptotic and the method of ordinary asymptotic
expansions is that the smallness of the remainder is understood
in a different way.
Usually, the remainder is assumed to be small in some uniform sense
with sufficient accuracy. Here we assume exactly the same but in the
sense of $O_{\mathcal{D}'}$.

These two remarks allow us to construct asymptotic formulas
that describe the propagation and interaction of nonlinear waves.
Some results obtained by this method
(for instance, results concerning the collision, formation, and
destruction of shock waves)
were published recently in the papers \cite{Dan,DaSh1,DaSh2,DaSh3}.

In the present paper we use this method for solving
the problem of solitary waves (solitons) interaction
for KdV type equations with small dispersion,
including the integrable case.
Here we encounter a problem of constructing
a suitable definition of the weak asymptotic solution
in the case where the dispersion tends to zero,
which is well known
as the zero limit dispersion problem \cite{LaxLev1,LaxLev2}.
However,
in contrast to \cite{LaxLev1,LaxLev2}, we deal not with oscillating solutions
of the KdV equation or even with more general solutions, but we
consider a special class of solutions, i.e., solitons.
Therefore, in our case the zero dispersion limit leads to a
system of differential equations, instead of
integro-differential equations, obtained in \cite{LaxLev1,LaxLev2}.

Consider in detail problems related to the definition
of the weak asymptotic
$\pmod{O_{\mathcal{D}'}(\epsilon^2)}$ solution $u^*_\epsilon$ for the KdV
equation. According to the notation in item 2), we define it so that
\begin{equation}
L_{\text{KdV}}u^*_\epsilon=O_{\mathcal{D}'}(\epsilon^2).
\tag {0.6}
\end{equation}
However, we note that (0.6) implies the loss of
distinction between the weak asymptotic solutions of the KdV and
Hopf equations. Indeed, since for any generalized functions $f$,
we obviously have the estimate
$$
\epsilon^2 f_{xxx}=O_{\mathcal{D}'}(\epsilon^2),
$$
we have
$$
L_{\rm KdV} u^*_\epsilon
=\frac{\partial u^*_\epsilon}{\partial t}+\frac{\partial(u^*_\epsilon)^2}{\partial x}
+O_{\mathcal{D}'}(\epsilon^2)=: L_H u^*_\epsilon+O_{\mathcal{D}'}(\epsilon^2),
$$
where
$$
L_Hu^*_\epsilon=\frac{\partial u^*_\epsilon}{\partial t}+\frac{\partial (u^*_\epsilon)^2}{\partial x}.
$$
Next, let us try to construct a weak asymptotic solution in the
sense of (0.6) in the simplest one-soliton situation.
We write a smooth ansatz in the form
\begin{equation}
\label{9}
u^*_\varepsilon= u_0(x,t)
+g(t)\omega\Big(\beta\frac{x-\phi}{\varepsilon} \Big)
+ e(x,t){\varepsilon }
\omega_0\Big(\beta\frac{x-\phi}{\varepsilon} \Big), \quad \varepsilon >0,
\tag {0.7}
\end{equation}
where $u_0(x,t)$, \ $g(t)$, \ $e(x,t)$, \ $\phi(t)$ are the desired
smooth functions, and, according to the results obtained in \cite{MasOm1,MasOm2},
$\beta=\sqrt{g/6}$.
The function $\omega(z)$ is a solution of the problem (0.5),
and
$\omega_0(z) \in C^{\infty}(\mathbb{R})$
satisfies the conditions
$\lim_{z \to +\infty}\omega_0(z)= 1$,
$\lim_{z \to -\infty}\omega_0(z)= 0$,
and $\big|\frac{d^\alpha\omega_0}{dz^\alpha}\big|
\leq C_\alpha(1+|z|)^{-3}$.
Then we have in the sense of ${\mathcal{D}}'(\mathbb{R})$
\begin{gather*}
\omega\Big(\beta\frac{x-\phi}{\varepsilon} \Big)
=\frac{\varepsilon}{\beta} \delta(x-\phi)
+ O_{\mathcal{D}'}(\varepsilon ^2), \\
\omega_0\Big(\beta\frac{x-\phi}{\varepsilon} \Big)
= \theta(x-\phi)
+ O_{\mathcal{D}'}(\varepsilon), \quad \varepsilon  \to +0,
\end{gather*}
where $\theta (z)$ is the Heaviside function.

  The system of equations for the functions $u_0(x,t)$,  $g(t)$,  $e(x,t)$,
$\phi(t)$ follows from (0.6) (this system is derived
in detail in \cite{DaSh3}) and has the following form:
(a comparative analysis of systems (0.8), (0.10), and (0.11)
was performed in \cite{DaSh3})
\begin{equation}
\begin{gathered}
u_{0t}+(u^2_0)_x =0,\\
\phi_{t} - 2u_0(\phi(t),t) - \frac{2}{3}g(t) =0, \\
e(\phi(t),t)-\frac{3\sqrt{6}}{2}g_t(t)/g^{3/2}(t) = 0,\\
\big(e_{t}(x,t) + 2(u_0(x,t)e(x,t))_x \big)\Big|_{x > \phi(t)} = 0.
\end{gathered}
\tag {0.8}
\end{equation}

    It is easy to verify that under the condition $g>0$ (which is an
analog of the admissibility condition in the theory of shock waves)
the solution of system (0.8) exists on any interval $t\in[0, T]$
such that the smooth solution $u_0$ of the Hopf equations exists on
this interval.

    System (0.8) can be solved in the following way: first, one finds
the smooth solution of the Hopf equation, next, one finds the function
$e(x,t)$ from the last equations (which is uniquely solvable in view
of the inequality $2u_0(\phi,t)<\phi_t$), then one finds the (positive)
function $g(t)$ from the next to the last equation, and finally, one
finds the function $\phi(t)$.

    Note that system (0.8) contains no obstacles to setting
$e(x,t)= 0$. If so, $g(t)={\mbox{const}}$ in the case of an
arbitrary (nonconstant) background function $u_0(x,t)$. But this
conclusion is contrary to well known properties of soliton
solutions of the KdV equation \cite{MasOm1,MasOm2}.

    Moreover, under our notation, the weak asymptotic of the
asymptotic one-soliton solution to the KdV equation, constructed
by  Maslov and Omelyanov \cite{MasOm1,MasOm2}, has the form
\begin{equation}
\begin{aligned}
u^*_{1,\varepsilon }(x,t)
=& u_{01}(x,t)+g_1(t) \omega\big(\beta\frac{x-\phi_1}{\varepsilon}\big)\\
&+ e_1(x,t){\varepsilon }
\Big[1-\omega_0\big(\beta\frac{x-\phi_1}{\varepsilon} \big)\Big]
+O_{\mathcal{D}'}(\epsilon^2), \quad \varepsilon  \to+0.
\end{aligned} \tag {0.9}
\end{equation}
In other words, in the case (0.7), the ``shock wave" with a
small amplitude \break ${\varepsilon }e(x,t)\theta(x-\phi_1(t))$
propagates {\it in front of the soliton\/}
${\varepsilon }\delta(x-\phi_1(t))$, but in the asymptotic one-soliton
solution constructed in \cite{MasOm1,MasOm2} the small shock wave
${\varepsilon }e_1(x,t)[1-\theta(x-\phi_1(t))]$ arises
{\it behind the soliton}.

When we apply (0.6) to the
asymptotic solution obtained in \cite{MasOm1,MasOm2}, whose weak asymptotic
yields (0.9), we obtain the following system of equations \cite{DaSh3}
\begin{equation}
\begin{gathered}
u_{01t}+(u^2_{01})_x =0,\\
\phi_{1t} - 2u_{01}(\phi_1(t),t) - \frac{2}{3}g_1(t) =0, \\
e_1(\phi(t),t)+\frac{3\sqrt{6}}{2}g_{1t}(t)/g^{3/2}_1(t) =0,\\
\big(e_{1t}(x,t) + 2(u_{01}(x,t)e_1(x,t))_x \big)\Big|_{x < \phi_1(t)}=0.
\end{gathered} \tag {0.10}
\end{equation}
The solution of this  system for $g_{1t}(t)\ne 0$ is not
uniquely determined by the initial conditions $e_1(x,0)$ for
$x \le \phi_1(0)$, since the velocity along the characteristic
($\dot x= 2u_{01}(x(t),t)$) is less (for $g_{1}(t)>0$) than that
of the soliton $\phi_{1t}=2u_{01}(\phi_1(t),t)+\frac{2}{3}g_1(t)$.

    Thus, the assumption that the structure of the solution to the KdV
equation is specified by (0.9) due to (0.6)
leads to an ill-posed Cauchy problem (with a non-unique
solution) for the functions $u_{01}(x,t)$,  $g_1(t)$,  $e_1(x,t)$,
 $\phi_1(t)$.

    On the other hand, the complete system of equations obtained
in \cite{MasOm1,MasOm2}  for these functions has the form
\begin{equation}
\begin{gathered}
u_{01t}+(u^2_{01})_x = 0,\\
\phi_{1t} - 2u_{01}(\phi_1(t),t) - \frac{2}{3}g_1(t) =0, \\
e_1(\phi(t),t)+\frac{3\sqrt{6}}{2}g_{1t}(t)/g_1^{3/2}(t) =0,\\
\big(e_{1t}(x,t) + 2(u_{01}(x,t)e_1(x,t))_x \big)\Big|_{x<\phi_1(t)}=0,\\
g_{1}(t) + 2u_{01}(\phi_1(t),t) = {\mbox{const}}.
\end{gathered}
\tag {0.11}
\end{equation}
It is evident that this system differs from system (0.10)
obtained from (0.6) by the additional equation
$g_{1}(t)+2u_{01}(\phi_1(t),t)=g_{1}(0)+2u_{01}(\phi_1(0),0)$.
The presence of this equation implies that system (0.11) splits
into the two systems
\begin{equation}
\begin{gathered}
u_{01t}+(u^2_{01})_x = 0,\\
\phi_{1t} - 2u_{01}(\phi_1(t),t) - \frac{2}{3}g_1(t) = 0, \\
g_{1}(t) + 2u_{01}(\phi_1(t),t) = {\mbox{const}},
\end{gathered}
\tag {0.12}
\end{equation}
and
\begin{gather*}
\big(e_{1t}(x,t) + 2(u_{01}(x,t)e_1(x,t))_x \big)\Big|_{x<\phi_1(t)}=0,
\tag {0.13}\\
e_1(\phi(t),t)+\frac{3\sqrt{6}}{2}g_{1t}(t)/g_1^{3/2}(t) = 0,
\tag {0.14}
\end{gather*}
and equality (0.14) is the boundary condition for equation (0.13),
which turns the Cauchy problem for equation (0.13) into the well-posed
one (the Cauchy condition, in view of (0.9),
has the form
$e_1(x,0)= e_1^0(x)[1-\theta(x-\phi_1(0))]$).

    Moreover, if equation (0.13) is considered formally in the
domain $x>\phi_1(t)$, which corresponds to the solution structure given by
formula (0.7), then the ``redundant" condition
$$
e(\phi(t),t)-\frac{3\sqrt{6}}{2}g_{t}(t)/g^{3/2}(t) = 0,
$$
analogous to (0.14), overdetermines the problem.

     Thus, the weak asymptotic corresponding to the asymptotic solution
of the Cauchy problem for the KdV equation constructed in \cite{MasOm1,MasOm2} 
cannot be derived from the solution to the KdV equation with the help of
(0.6), and vice versa.

 Why is it so? The essence of the matter lies in the definition of
weak (generalized) solution to nonlinear equation. It turns out that
the definition of the weak (generalized) solution to nonlinear equation
depends on the structure of the kernel of the operator adjoint to
the linearized operator of the initial differential equation which
arises when constructing the smooth asymptotic. This construction
of the definition of weak solutions was previously discussed
in \cite{DaSh3,DaOmRa}.

 In the present paper we do not come into details of construction of
the definition of the weak solution to our problem. We just point out
that, in terms of this construction, the KdV equation is analogous to
the phase field system discussed in \cite{DaOmRa}.
The difference is that for the KdV equation the kernel of
the adjoint operator mentioned above is {\it two-dimensional}.
Therefore, relation (0.6) is not sufficient for the correctness
and it has to be supplied with another condition.

 At least for special soliton type initial data for the KdV
equation, we can give a definition of a weak solution
admitting the zero dispersion limit.

This definition was first presented in \cite{DaSh3}.

\begin{definition} \rm  \label{def0.1}
A function $u_\epsilon(x,t)$ belonging to
$C^\infty([0,T]\times \mathbb{R}^1_x)$ for $\epsilon>0$
and to $C(0,T;\mathcal{D}'(\mathbb{R}^1_x))$ uniformly in $\epsilon\geq0$
is called a {\it weak asymptotic solution\/} of the KdV type
equation $Lu=0$ if the following two relations are satisfied:
\begin{gather*}
Lu_\epsilon=O_{\mathcal{D}'}(\epsilon^2),\\
u_\epsilon Lu_\epsilon=O_{\mathcal{D}'}(\epsilon^2),
\end{gather*}
where, uniformly with respect to $\epsilon$,
the right-hand sides are weakly continuous and
weakly piecewise smooth functions of $t$.
\end{definition}

\begin{remark} \rm  \label{rmk1.1}
In fact, this definition means that we do not refuse to use
relation (0.6), but we impose additional requirements on the
right-hand side of (0.6).
\end{remark}

Of course, these relations can be written in the form of
usual integral identities (see below).
The condition that the remainders $O_{\mathcal{D}'}(\epsilon^2)$ are weakly
continuous in $t$  allows one to pass from Definition \ref{def0.1} to
integral identities in the usual sense.
This theme was discussed in \cite{DaSh3} in detail.
In particular, in \cite{DaSh3} it was shown that the passage to the limit
within the framework of this definition for single exact and
distorted solitons
leads to the same system of ordinary differential equations,
which was derived by using the classical asymptotic method \cite{MasOm1,MasOm2}.
The consequences of the use of this definition
in the case of two solitons (the case of soliton interaction)
are the main point of the present paper.

It should be noted that we do not study some weak solution of
(0.3) in the Banach space but construct a solution of a specific
structure and describe some fine properties of its dynamics.
This is impossible by using methods traditional to the theory
of nonlinear PDE (e.g., see \cite{BonaSa,LinSc}).
Earlier, such results were obtained only for integrable problems
by the inverse scattering transform method
or for problems similar in a sense to integrable ones
and by using the inverse scattering transform method
again \cite{AblS,MasOm1,MasOm2,McSc}.

We also note that, even in the case of a solitary soliton,
our approach allows obtaining results that are new
as compared to the well-known results.
Namely, the asymptotic methods known today \cite{AblS,MasOm1,MasOm2,McSc,MolVak}
allows calculating the dynamics of a distorted soliton
if there is a perturbed right-hand side and/or a variable
background.
However, in this case there arise very rigid restrictions
on the Cauchy data.
In the integrable case, in particular, for the KdV equation,
it is possible to get rid of the restrictions on lower-order
terms (of order $O(\epsilon)$) in the initial data only  by using
very complicated constructions \cite{IlKal} based on the inverse
scattering transform method.
In our approach, only the leading term of the asymptotic
expansion is fixed and the construction is
reduced to simple algebraic calculations.

The text is organized as follows.
First, we present some auxiliary results,
which are useful for calculations in the weak asymptotic
method.
In the next section, we study the two-soliton solution of the KdV equation.
In the final section,
we consider non-integrable KdV type equations.

\section{Auxiliary formulas of the weak asymptotic method}%1.

\subsection*{Weak asymptotic expansions}%1.1
Let $\omega(\eta)$ be a continuous function decreasing sufficiently
fast as $|\eta|\to\infty$
(so that the integrals below were meaningful),
and let $\beta$ be independent of $x$.
Let us consider the expression $\omega(\beta x/\epsilon)$
as an element of $\mathcal{D}'$.
For any test function $\psi(x)$ we have
\begin{align*}
\big\langle\omega\big(\beta\frac{x}\epsilon\big),\psi(x)\big\rangle
&=\int_{\mathbb{R}^1}\omega\big(\beta\frac{x}\epsilon\big)\psi(x)\,dx
=\frac{\epsilon}{\beta}\int_{\mathbb{R}^1} \omega(\eta)\psi\big(\epsilon\frac{\eta}\beta\big)\,d\eta
\\
&=\sum^{n}_{k=0}\frac{\epsilon^{k+1}}{\beta^{k+1}k!}\psi^{(k)}(0)
\int_{\mathbb{R}^1}\omega(\eta)\eta^k\,d\eta +O(\epsilon^{n+2}).
\end{align*}
This formula is well known in the theory of algebras of
Colombeau generalized functions
as the momentum decomposition.

Using the notation $O_{\mathcal{D}'}(\epsilon^\alpha)$ introduced above,
we can rewrite this relation as follows:
\begin{equation}
\omega\big(\beta\frac{x}\epsilon\big)
=\sum^{n}_{k=0}\frac{\epsilon^{k+1}}{\beta^{k+1}k!}(-1)^k
\Omega_k\delta^{(k)}(x)+O_{\mathcal{D}'}(\epsilon^{n+2}),
\tag {1.1}
\end{equation}
where
$\Omega_k=\int \eta^k\omega(\eta)\,d\eta$ and $\delta(x)$ is the Dirac
delta function.
In what follows, we restrict ourselves to studying the leading
term in (1.1).
Then we have
$$
\omega\Big(\beta\frac{x}\epsilon\Big)
=\epsilon\frac{\Omega_0}{\beta}\delta(x)+O_{\mathcal{D}'}(\epsilon^2)
$$
in the general case and the last relation can be made more
precise by specifying the
properties of the function $\omega(\eta)$.
Namely, the term $O_{\mathcal{D}'}(\epsilon^2)$ has the form
$$
O_{\mathcal{D}'}(\epsilon^2)=-\epsilon^2\delta'(x)\beta^{-2}
\int_{\mathbb{R}^1} \eta\omega(\eta)\,d\eta +O_{\mathcal{D}'}(\epsilon^3).
$$
Thus if $\omega(\eta)$ is an even function, then
$$
\omega\Big(\beta\frac{x}\epsilon\Big)=\epsilon\frac{\Omega_0}{\beta}\delta(x)
+O_{\mathcal{D}'}(\epsilon^3).
$$
Now let us consider a more complicated expression
$$
f\Big(\omega_1\Big(\beta_1\frac{x-\varphi_1}\epsilon\Big)
+\omega_2\Big(\beta_2\frac{x-\varphi_2}\epsilon\Big)\Big),
$$
where $f(\eta)$ is a smooth function
and $\omega_i(\eta)$, $i=1,2$, possess the same properties
as $\omega(\eta)$,
$\beta_1,\beta_2=\mathop{\rm const}>0$, and $\varphi_1,\varphi_2$ are constant.
We have
\begin{align*}
&\big\langle
f\Big(\omega_1\big(\beta_1\frac{x-\varphi_1}{\epsilon}\big)
+\omega_2\big(\beta_2\frac{x-\varphi_2}{\epsilon}\big)\Big),\psi \big\rangle\\
&\quad
=\Big\{\beta_1(x-\varphi_1)=\epsilon \eta\Big\}\\
&\quad
=\big\langle\frac{\epsilon}{\beta_1}\int
f_1\Big(\omega_1(\eta)+\omega_2\big(\frac{\beta_2}{\beta_1}\eta
+\frac{\beta_2}\epsilon\Delta \varphi \big)\Big)\,d\eta \delta(x-\varphi_1)
+f(0),\psi \big\rangle+O(\epsilon^2),
\end{align*}
where $f_1(\eta)=f(\eta)-f(0)$ and $\Delta\varphi=\varphi_1-\varphi_2$.

On the other hand, introducing the variable $\eta$ by the formula
$\beta_2(x-\varphi_2)=\epsilon\eta$, we obtain
\begin{align*}
&f\Big(\omega_1\big(\beta_1\frac{x-\varphi_1}{\epsilon}\big)
+\omega_2\big(\beta_2\frac{x-\varphi_2}{\epsilon}\big)\Big)\\
&=f(0)+\frac{\epsilon}{\beta_2}\int_{\mathbb{R}^1}
f_1 \Big(\omega_1\big(\frac{\beta_1}{\beta_2}\eta-
\frac{\beta_1}\epsilon\Delta \varphi\big)
+\omega_2(\eta)\Big)\,d\eta\delta(x-\varphi_2)
+O_{\mathcal{D}'}(\epsilon^2).
\end{align*}
It is easy to see that the coefficients
of $\delta(x-\varphi_1)$ and $\delta(x-\varphi_2)$
in the last formulas are equal to each other,
\begin{align*}
&\beta^{-1}_1\int_{\mathbb{R}^1}  f_1
\Big(\omega_1(\eta)
+\omega_2\big(\frac{\beta_2}{\beta_1}\eta+
\frac{\beta_2}\epsilon\Delta \varphi\big)\Big)\,d\eta\\
&=\beta^{-1}_2\int_{\mathbb{R}^1}  f_1
\Big(\omega_1\big(\frac{\beta_1}{\beta_2}\eta-\frac{\beta_1}\epsilon
\Delta \varphi\big)
+\omega_2(\eta)\Big)\,d\eta := \Omega_{\Delta \varphi}.
\end{align*}
This implies the relation
\begin{align*}
&f\Big(\omega_1\big(\beta_1\frac{x-\varphi_1}{\epsilon}\big)
+\omega_2\big(\beta_2\frac{x-\varphi_2}{\epsilon}\big)\Big)\\
&=f(0)+\epsilon\Omega_{\Delta \varphi}
(\lambda\delta(x-\varphi_1)+\nu\delta(x-\varphi_2))
+O_{\mathcal{D}'}(\epsilon^2),
\end{align*}
where $\lambda$ and $\nu$ are arbitrary constants,
$\lambda+\nu=1$.

\subsection*{Asymptotic linear independence}%1.2
If we consider linear combinations of generalized functions with
accuracy up to $O_{\mathcal{D}'}(\epsilon^\alpha)$,
then we need to introduce the notion of linear independence.
Consider the relation
$$
g_1\delta(x-\varphi_1)+g_2\delta(x-\varphi_2)
=O_{\mathcal{D}'}(\epsilon^\alpha),\quad \alpha>0,\quad \varphi_1\ne \varphi_2,
$$
where $g_i$ are independent of $\epsilon$.
Obviously, we obtain the relations
$$
g_i=O_{\mathcal{D}'}(\epsilon^\alpha),\quad i=1,2,
$$
which, by virtue of our assumption, imply
$$
g_i=0,\quad i=1,2.
$$

Everything is different
if we assume that the coefficients $g_i$ can depend on $\epsilon$.
Here we consider only a special case of such dependence, which
we will use later.
Namely, let
$$
g_i=A_i+S_i\Big(\frac{\Delta \varphi}{\epsilon}\Big), \quad i=1,2,
$$
where $A_i$ are independent of $\epsilon$ and continuous
in $\varphi_1$ and $\varphi_2$ and $S_i(\tau)$ decrease
sufficiently fast as $|\tau|\to\infty$.

Let us find out what properties of the coefficients $g_i$
follow from the relation
$$
g_1\delta(x-\varphi_1)+g_2\delta(x-\varphi_2)=O_{\mathcal{D}'}(\epsilon).
$$
Applying both sides of the equality
to a test function $\psi$,
we obtain
$$
g_1\psi(\varphi_1)+g_2\psi(\varphi_2)=O(\epsilon)
$$
or, which is the same,
\begin{equation}
[A_1\psi(\varphi_1)+A_2\psi(\varphi_2)]
+[S_1\psi(\varphi_1)+S_2\psi(\varphi_2)]
=O(\epsilon).
\tag {1.2}
\end{equation}
Let us consider the expression in the second brackets.
Using Taylor's formula, we obtain
$$
[S_1\psi(\varphi_1)+S_2\psi(\varphi_2)]
= S_1\psi(\varphi_1)+S_2\psi(\varphi_1)
+S_2(\varphi_2-\varphi_1)\psi'(\varphi_1+\theta \varphi_2).
$$
Now we see that
$$
S_2\Big(\frac{\Delta \varphi}{\epsilon}\Big)(\varphi_2-\varphi_1)
=\{-\rho S_2(\rho)\}\big|_{\rho=\Delta \varphi/\epsilon}\cdot \epsilon=O(\epsilon)
$$
since the function $\rho S_2(\rho)$ is assumed to be bounded
uniformly in $\rho\in\mathbb{R}^1$.
So we can rewrite relation (1.2) as
$$
A_1\psi(\varphi_1)+A_2\psi(\varphi_2)+(S_1+S_2)\psi(\varphi_1)=O(\epsilon).
$$

Our next goal is to obtain from this relation
the condition for the solution to be {\it uniform\/}
under the choice of $\varphi_1$ and $\varphi_2$.
Let $|\Delta\varphi| \geq \mathop{\rm const}$.
Then $|S_1|=O(\epsilon)$ and $|S_2|=O(\epsilon)$,
and in the usual way  we obtain
$$
A_1=0,\quad A_2=0.
$$
Further we obtain the relation
$S_1+S_2=0$.
It follows from the continuity that $A_1=0$ and $A_2=0$
also for $\Delta\varphi=0$. Hence we finally obtain
\begin{equation}
A_1=0,\quad A_2=0,\quad S_1+S_2=0.
\tag {1.3}
\end{equation}
Another method for analyzing relation (1.2) is the following.
We choose a point $x^*$. Then for any test function $\psi(x)$ we
have
\begin{align*}
\langle S_1\delta(x-\varphi_1),\psi\rangle
+\langle S_2\delta(x-\varphi_2),\psi\rangle
&=S_1\psi(x^*)+S_2\psi(x^*)+O(\epsilon)\\
&=\langle(S_1+S_2)\delta(x-x^*),\psi\rangle+O(\epsilon),
\end{align*}
which implies
\begin{align*}
&(A_1+S_1)\delta(x-\varphi_1)+(A_2+S_2)\delta(x-\varphi_2)\\
&=A_1\delta(x-\varphi_1)+A_2\delta(x-\varphi_1)
+(S_1+S_2)\delta(x-x^*)+O_{\mathcal{D}'}(\epsilon).
\end{align*}
Following the above argument, we again obtain relations (1.3).

\section{$\epsilon$--$\delta$ soliton interaction in the KdV model}%2.

In this and in the following sections,
by using the ideas briefly considered above,
we qualitatively describe
the interaction of $\epsilon$--$\delta$-solitons
(of narrow solitary waves)
in models governed by KdV type equations.
>From methodological considerations,
we first study the simplest case,
i.e.,
the KdV equation with small dispersion
\begin{equation}
u_t+(u^2)_x+\epsilon^2u_{xxx}=0
\tag {2.1}
\end{equation}
supplemented with the initial data,
which are
the superposition of two solitary waves
\begin{equation}
u\big|_{t=0}=u^0(x,\epsilon),\quad
u^0=A_1\omega\big(\beta_1\frac{x-x^0_1}{\epsilon}\big)
+A_2\omega\big(\beta_2\frac{x-x^0_2}{\epsilon}\big).
\tag {2.2}
\end{equation}
We set $\beta_i=\sqrt{A_i/12}$,
$\omega(\eta)=1/2\cosh^2(\eta)$, and assume that
$$
x^0_2<x^0_1,\quad A_2>A_1>0
$$
so that there exist a time moment $t^*$
at which the trajectories of the solitary waves
corresponding to \thetag{2.2}
intersect:
$V_1t^* + x^0_1=V_2t^* + x^0_2$, $V_i=4\beta^2_i$.

Needless to say, for the KdV equation,
as well as for the modified KdV (MKdV) equation,
the behavior of the solution is well known:
the solitons interact passing through each other
without changing their shapes and, after the interaction,
they are shifted by some distance.
This assertion can be obtained directly, for instance,
from the formula for the exact two-soliton solution
of problem \thetag{2.1}, \thetag{2.2}
derived by the inverse scattering transform method
(e.g., see \cite{AblS})
\begin{equation} \tag {2.3}
\begin{gathered}
u_{\text{sol}}=6\beta^2_1(\beta^2_2-\beta^2_1)v_1(\eta_1,\rho)
+6\beta^2_2(\beta^2_2-\beta^2_1)v_2(\eta_2,\rho),\\
v_1=\big\{\beta_1\sinh(\beta_1\eta_1-\mu)
-\beta_2\coth\big(\beta_2(\eta_1-\rho)+\mu\big)
\cosh(\beta_1\eta-\mu)\big\}^{-2}, \\
v_2=\big\{\beta_2\cosh(\beta_2\eta_2+\mu)
-\beta_1\tanh\big(\beta_1(\eta_2+\rho)-\mu\big)
\sinh(\beta_2\eta_2+\mu)\big\}^{-2},
\end{gathered}
\end{equation}
where
$\eta_i=\big(x-\varphi_{i0}(t)\big)/\epsilon$,
$\varphi_{i0}=V_it+x^0_i$,
$\rho=\big(\varphi_{20}(t)-\varphi_{10}(t)\big)/\epsilon$, and
\begin{equation}
\mu=\frac12\ln\frac{\beta_2-\beta_1}{\beta_2+\beta_1}
\tag {2.4}
\end{equation}
is the displacement of the soliton trajectories
after the interaction.

It is well known that as $\rho\to-\infty$
(that is, for $\varphi_{20}-\varphi_{10}<0$, $\epsilon\to0$,
i.e., prior to the interaction)
and as $\rho\to\infty$
(that is, for $\varphi_{20}-\varphi_{10}>0$, $\epsilon\to0$, i.e.,
after the interaction),
formula (2.3)
is the sum of isolated solitons
with accuracy up to terms admitting the estimate $O(\epsilon^N)$ for
any $N>0$.
However,
after the interaction, some constants are added
to the arguments of the solutions.

In fact, the problem of describing the soliton interaction
(as soon as the qualitative mechanism of interaction becomes
clear) is just the calculation of these constants.

As mentioned above, the weak asymptotic method
permits us to obtain a description of the single soliton
dynamics in an extremely simple way.
There is a natural question of whether we can describe
the interaction of solitons
passing  to the limit in the weak sense as $\epsilon\to0$.

In what follows, we present some results obtained in this way.
According to our approach,
for problem (2.1), (2.2) we present the two-soliton
weak asymptotic solution $\bmod O_{\mathcal{D}'}(\epsilon^2)$
in the form
\begin{equation} \tag {2.5}
\begin{gathered}
u=g_1(\tau)\omega\Big(\beta_1\frac{x-\phi_1(t,\tau,\epsilon)}{\epsilon}\Big)
+g_2(\tau)\omega\Big(\beta_2\frac{x-\phi_2(t,\tau,\epsilon)}{\epsilon}\Big),\\
g_i=A_i+S_i(\tau),\quad
\phi_i=\varphi_{i0}(t)+\epsilon\varphi_{i1}(\tau),\quad
\tau=\beta_1\big(\varphi_{20}(t)-\varphi_{10}(t)\big)/\epsilon,
\end{gathered}
\end{equation}
where $\varphi_{i0}(t)=V_it+x^0_i$
so that $x=\varphi_{i0}(t)$ is the trajectory of motion
of the solitary wave with constant amplitude $A_i$,
the function $\tau=\tau(t,\epsilon)$ has the meaning of ``fast'' time,
$\tau(t^*,\epsilon)=0$, $S_i$ and $\varphi_{i1}$
are corrections to the amplitudes and the phases
rapidly varying during the time of interaction
with exponential velocities:
\begin{equation} \tag {2.6}
\begin{gathered}
S_i(\tau)\to0\quad\text{as } \tau\to\pm\infty,\\
\varphi_{i1}(\tau)\to0\quad\text{as } \tau\to-\infty,
\quad
\varphi_{i1}(\tau)\to\varphi^\infty_{i1}=\mathop{\rm const}
\quad\text{as } \tau\to+\infty.
\end{gathered}
\end{equation}
It turned out that
(in contrast to the inverse problem method)
the construction of the solution significantly
depends on the ratio of the parameters $\beta_1$ and $\beta_2$,
i.e., on  $\theta=\beta_1/\beta_2$, and this ratio can be easily
recalculated for the ratio of soliton amplitudes.
Namely, it turned out that the solution in the form of (2.5) can
be constructed for $0<\theta<1/2$.
For $\theta=1/2$ the formulas degenerate,
and this case must be studied separately.
For $\theta\in(\frac12,1)$, the functions contained in the
expression on the right-hand side of (2.5) become complex-valued,
and, in order to deal with real functions,
we must seek the solution in the form
$$
u=g_1(\tau)\omega\Big(\beta_1\frac{x-\phi_1(t,\tau,\epsilon)}{\epsilon}\Big)
+g_2(\tau)\omega\Big(\beta_2\frac{x-\phi_2(t,\tau,\epsilon)}{\epsilon}\Big)
+\text{c.c.}.
$$

It is clear that these additional terms make our calculations
much more complicate. So, in what follows,
we restrict our study of interaction to the case
$0<\theta<1/2$.
Strictly speaking, we can prove that there exists a
$\theta^*<1/2$  such that our statements are true for
$0<\theta<\theta^*$.
But numerical calculations show that $\theta^*=1/2$.


Under this restriction on the initial amplitudes of solitons,
we construct an asymptotic solution of the KdV equation
and, and in the next section, for the KdV type equations.

Our aim is to  reconstruct all the qualitative properties
of soliton interaction obtained by the inverse problem method:
\\
1) solitons pass through each other without changing
their structure,
\\
2) the result of interaction is a phase shift.

In what follows, such a mechanism of interaction of
solitary waves  will be called the {\it soliton scenario of
interaction}.
The description of interaction is reduced to the study of
autonomous first-order differential equations.
Moreover, these equations are of similar form
for both the KdV equation and the non-integrable KdV type equations.

One cannot solve these equations explicitly, however, for small
values of the parameter $\theta$, it is possible to obtain
enough information about the solution and to calculate the phase
shifts. The formulas thus obtained differ from (2.4).
The matter is that, as will be shown below,
the KdV type equations with small dispersion are not well
defined in the asymptotic sense. Namely, an arbitrary  small
perturbation of the solution can lead to
a contribution to the solution of order $O_{\mathcal{D}'}(\epsilon^2)$ and
thus to change the phase shift.

The main result of this section,
which we obtain without using the inverse scattering transform
method,  is the following theorem.

\begin{theorem} \label{thm2.1}
Let $\theta=\beta_1/\beta_2\in(0,1/2)$, then
the interaction of $\epsilon$--$\delta$-solitons of the KdV equation
follows the soliton scenario.
The weak asymptotic of the solution
of problem \thetag{2.1}, \thetag{2.2} is asymptotically
not unique in the terms of order $O_{\mathcal{D}'}(\epsilon^2)$.
\end{theorem}

We prove this assertion in two steps.
First, we construct the weak asymptotic
modulo $O_{\mathcal{D}'}(\epsilon^2)$
of the solution of problem \thetag{2.1}, \thetag{2.2}
in the form  (2.5), (2.6)
At the second stage of the proof, we show that a perturbation
(arbitrary small in the sense of $\mathcal{D}'$) of the leading term of
the asymptotic of (2.5) changes the constants
$\varphi^{\infty}_{i1}$.

Prior to proving the theorem, let us discuss the leading term of
the weak asymptotic of \thetag{2.5} and compare it with the exact
solution \thetag{2.3}.
First, recall that we understand
the solution of problem \thetag{2.1}, \thetag{2.2}
in the weak sense (see Introduction),
that is,
by a weak asymptotic modulo $O_{\mathcal{D}'}(\epsilon^2)$ of the solution
we mean a function $u=u(x,t,\epsilon)$ such that
for any test function $\psi=\psi(x)$
the following relations are satisfied:
\begin{equation} \tag {2.7}
\begin{gathered}
\frac{d}{dt}\int u\psi\,dx-\int u^2\frac{\partial \psi}{\partial x}\,dx
=O(\epsilon^2),\\
\frac{d}{dt}\int u^2\psi\,dx
-\frac43\int u^3\frac{\partial \psi}{\partial x}\,dx
+3\int\Big(\epsilon\frac{\partial u}{\partial x}\Big)^2 \frac{\partial \psi}{\partial x}\,dx
=O(\epsilon^2), \\
\int \big(u\big|_{t=0}-u^0(x,\epsilon)\big)\psi\,dx
=O(\epsilon^2).
\end{gathered}
\end{equation}
Here and in the following, the notation $\int\dots\,dx$
means the integration over $\mathbb{R}^1$.

The most unexpected result obtained
while proving the theorem is that the following fact.

\begin{theorem} \label{thm2.2}
The condition $\theta\in(0,1/2)$, the conservation laws
\begin{gather}
\frac{d}{dt}\int u\,dx =0, \tag {2.8}\\
\frac{d}{dt}\int u^2\,dx =0, \tag {2.9}
\end{gather}
and the energy relations
\begin{gather}
\frac{d}{dt}\int xu\,dx-\int u^2\,dx =0, \tag {2.10}\\
\frac{d}{dt}\int xu^2\,dx -\frac 43\int u^3\,dx
+3\int\Big(\epsilon\frac{\partial u}{\partial x}\Big)^2\,dx =0 \tag {2.11}
\end{gather}
are necessary and sufficient conditions
for the function \thetag{2.5} to be a weak asymptotic
of the solution of problem \thetag{2.1}, \thetag{2.2}
in the sense of \thetag{2.7}.
\end{theorem}

Let us discuss the structure of the asymptotic solution  of \thetag{2.5}.
Obviously, the above formulas imply that
the function \thetag{2.5} is
the sum of solitary waves until they interact
$$
u=\sum^{2}_{i=1}A_i\omega\Big(\beta_i\frac{x-\varphi_{i0}(t)}{\epsilon}\Big)
+O(\epsilon^\infty),\quad t<t^*,
$$
and the sum of solitary waves with phase shift
$$
u=\sum^{2}_{i=1}A_i\omega\Big(\beta_i\frac{x-\varphi_{i0}(t)}{\epsilon}
-\beta_i\varphi^\infty_{i1}\Big)
+O(\epsilon^\infty),\quad t>t^*,
$$
after the interaction.
We stress
that the resulting phase shifts $\beta_i\varphi^\infty_{i1}$
satisfy the same relation
\begin{equation}
\beta_1\varphi^\infty_{11}+\beta_2\varphi^\infty_{21}=0
\tag {2.12}
\end{equation}
as those for the exact two-soliton solution (2.3).
In the sense of $\mathcal{D}'$,
for the function \thetag{2.5} uniformly in $t$
we have
\begin{equation}
u=\epsilon\sum^{2}_{i=1}\frac{g_i(\tau)}{\beta_i}\delta(x-\phi_i)
+O_{\mathcal{D}'}(\epsilon^3).
\tag {2.13}
\end{equation}
However,
$\epsilon\delta(x-\phi_i)=\epsilon\delta(x-\varphi_{i0})+O_{\mathcal{D}'}(\epsilon^2)$
and we see that only soliton components of the amplitude $S_i$
are essential, while the role of phase shifts $\epsilon\varphi_{i1}$
seems to be unessential.
Moreover,
a similar conclusion can readily be obtained by calculating
the weak asymptotic of the exact solution \thetag{2.3}:
\begin{equation}
u_{\text{sol}}=\epsilon\sum^2_{i=1}\frac{P_i(\tau)}{\beta_i}
\delta(x-\varphi_{i0}(t))+O_{\mathcal{D}'}(\epsilon^2),
\tag {2.14}
\end{equation}
where
$$
P_i=6\beta^2_i(\beta^2_2-\beta^2_1)
\int v_i(\eta,\tau/\beta_1)\,d\eta\to A_i\quad
\text{as } \tau\to\pm\infty.
$$
Nevertheless, the phase shifts play an important role.
To verify this, we calculate the weak asymptotic for the
derivatives of the function \thetag{2.5}.
Taking into account that fact that $S_i$ and $\varphi_{i1}$ vary
rapidly with time, we obtain
\begin{align*}
\frac{\partial u}{\partial t}
&=\sum^2_{i=1}\Big\{\frac{c}{\beta_i}
\frac{\partial S_i}{\partial \tau}\delta(x-x^*)
-\epsilon\frac{A_i}{\beta_i}V_i\delta'(x-\varphi_{i0})\\
&\quad
-c\epsilon\frac{\partial}{\partial\tau}
\Big[V_i\frac{\tau}{c}\frac{S_i}{\beta_i}
+\frac{g_i}{\beta_i}\varphi_{i1}\Big]\delta'(x-x^*)\Big\}
+O_{\mathcal{D}'}(\epsilon^2),
\end{align*}
where
$c=\beta_1(V_2-V_1)$ and $x^*=\varphi_{i0}(t^*)$ is the point of
intersection of the soliton trajectories.
Similar calculations for the exact solution \thetag{2.3} lead to
the expansion
\begin{align*}
\frac{\partial u_{\text{sol}}}{\partial t}
&=\sum^2_{i=1}\Big\{\frac{c}{\beta_i}
\frac{\partial P^s_i}{\partial \tau}\delta(x-x^*)
-\epsilon\frac{A_i}{\beta_i}V_i\delta'(x-\varphi_{i0})\\
&\quad
-c\epsilon\frac{\partial}{\partial\tau}
\big[V_i\frac{\tau}{c}\frac{P^s_i}{\beta_i}
+\frac{D_i}{\beta_i}\big]\delta'(x-x^*)\Big\}
+O_{\mathcal{D}'}(\epsilon^2),
\end{align*}
where $D_i=6\beta^2_i(\beta^2_2-\beta^2_1)
\int\eta v_i(\eta,\tau/\beta_1)\,d\eta$
and we set $P_i=A_i+P^s_i$
so that $P^s_i\to0$ as $\tau\to\pm\infty$.
Comparing these two formulas,
as well as (2.13) and (2.14),
we readily see that
the varying phase shifts of the function \thetag{2.5} play the
same role as the odd (with respect to $\eta$) parts
of the functions $v_i(\eta,\tau/\beta_i)$ in the exact solution.

If we discuss the interaction process and compare
the functions \thetag{2.3} and \thetag{2.5} in the pointwise sense,
then we see that these functions are, of course, different.
Nevertheless, they possess an important common property:
the soliton type (with respect to time) components of the amplitudes
compensate one another.
We mean that uniformly in $\tau\in\mathbb{R}^1$, the relation
\begin{equation}
S_1(\tau)/\beta_1+S_2(\tau)/\beta_2=0
\tag {2.15}
\end{equation}
holds for the weak asymptotic \thetag{2.5}
(we prove this later).
This coincides with the identity
$$
P_1(\tau)/\beta_1+P_2(\tau)/\beta_2=A_1/\beta_1+A_2/\beta_2,
$$
which can be derived directly for the exact solution
by rewriting formula (2.3) in the form
standard for the inverse problem method
or by using the conservation law \thetag{2.8}.

Prior to proving Theorem \ref{thm2.1},
we introduce some auxiliary functions
we need to describe the weak asymptotic \thetag{2.5}
exactly.
By $\lambda_j=\lambda_j(\sigma)$ and $L_i=L_i(\sigma)$
we denote the following convolutions depending on the parameter
$\theta=\beta_1/\beta_2$:
\begin{gather}
\begin{gathered}
\lambda_0=\frac1{a_2}\int\omega(\eta)\omega(\theta\eta+\sigma)\,d\eta,
\quad
\lambda_1=\frac1{a_2}\int\eta\omega(\eta)\omega(\theta\eta+\sigma)\,d\eta,
\\
\lambda_2=\frac1{a'_2}\int\omega'(\eta)\omega'(\theta\eta+\sigma)\,d\eta,
\quad
\omega'(z)=\partial\omega(z)/\partial z,
\end{gathered} \tag {2.16}\\
\begin{gathered}
L_1=\frac1{a_3}\int\Big\{
\theta(\theta+k)\omega(\theta\eta+\sigma)
+(1-k)\omega(\eta)\Big\}^3\,d\eta-(1+\theta^5), \\
L_2=\theta^3k(2\theta+k)+k^2-2k +2\theta^2(\theta+k)(1-k)\lambda_2,
\end{gathered} \tag {2.17}
\end{gather}
where $a_2, a_3, a'_2$ are the numbers
\begin{equation}
\begin{gathered}
a_2=\int\omega^2(\eta)\,d\eta=\frac13,\quad
a_3=\int\omega^3(\eta)\,d\eta=\frac2{15},\\
a'_2=\int\big(\omega'(\eta)\big)^2\,d\eta=\frac4{15},
\end{gathered}
\tag {2.18}
\end{equation}
and $k=k(\sigma)$ is the function determined for $\theta\in(0,1/2)$
by the expression
\begin{equation}
\begin{gathered}
k=\big\{(1-\theta)(1+\theta-\lambda_0\theta)-\sqrt{I}\big\}
(1+\theta-2\lambda_0\theta)^{-1},\\
I=(1+\theta)^2(1-\theta-\lambda_0\theta)^2
+2\lambda_0\theta^2(1+\theta)(1-2\theta).
\end{gathered} \tag {2.19}
\end{equation}
For $\theta=1/2$ we set
\begin{equation}
k=\lambda_0/(3-2\lambda_0).
\tag {2.20}
\end{equation}
Next, we introduce a function $\mathcal{F}=\mathcal{F}(\sigma)$ by the formula
\begin{equation}
\mathcal{F}=F/(1-\theta+\mathcal{J}'),
\tag {2.21}
\end{equation}
where
\begin{align*}
F&=\Big(\frac43\frac{a_3}{a^2_2}(L_1-L_2)+L_2\Big)
\big(\theta(1-\theta)^2\big)^{-1}-\mathcal{J}',
\\
\mathcal{J}'&=k\big(\frac{1-\theta-\theta^2}{\theta}-k\big)
+\big\{\big(\frac{1-\theta-\theta^2}{\theta}-2k\big)\sigma
+2\theta(1-\theta-2k)\lambda_1\big\}\frac{\partial k}{\partial\sigma}\\
&\quad
+2\theta(\theta+k)(1-k)\frac{\partial\lambda_1}{\partial\sigma}.
\end{align*}
Let us define a function $\sigma_0=\sigma_0(\tau)$ depending
also on the parameter $\theta$ as a solution of the problem
\begin{equation}
\frac{d\sigma_0}{d\tau}=\mathcal{F}(\tau+\sigma_0),\quad
\sigma_0\to0\quad\text{as } \tau\to-\infty.
\tag {2.22}
\end{equation}
Next, we show that problem \thetag{2.22} has
solution for all $\tau\in\mathbb{R}^1$.
Moreover,
\begin{equation}
\sigma_0\to\sigma^\infty_0(\theta)\quad\text{as } \tau\to+\infty
\tag {2.23}
\end{equation}
as a function of exponential type.


\begin{theorem} \label{thm2.3}
Let $\theta\in(0,1/2)$, then
the weak asymptotic solution of problem \thetag{2.1},
\thetag{2.2}
has the form \thetag{2.5},
where
\begin{gather}
S_1=\sqrt{A_1A_2}k(\tau+\sigma_0),\quad
S_2=-A_2k(\tau+\sigma_0),
\tag {2.24}\\
\begin{aligned}
\varphi_{11}&=\frac1{\beta_1(1+\theta)}
\big\{\tau k(\tau+\sigma_0)-\sigma_0(1-k(\tau+\sigma_0))\big\},\\
\varphi_{21}&=\frac1{\beta_2\theta(1+\theta)}
\big\{\theta \sigma_0+(\tau+\sigma_0)k(\tau+\sigma_0)\big\}.
\end{aligned}
\tag {2.25}
\end{gather}
\end{theorem}

\begin{remark} \label{rmk2.1} \rm
The above formulas readily imply that
\begin{gather*}
k(\sigma)\to0\quad\text{as }\sigma\to\pm\infty,
\quad
\varphi_{i1}(\tau)\to0\quad\text{as }\tau\to-\infty,
\quad i=1,2,
\\
\beta_1\varphi_{11}(\tau)\to-\sigma^\infty_0/(1+\theta),\quad
\beta_2\varphi_{21}(\tau)\to\sigma^\infty_0/(1+\theta)
\quad\text{as }\tau\to+\infty
\end{gather*}
as functions of exponential type.
\end{remark}

\begin{proof}[Proof of Theorem  \ref{thm2.3}]
We present the leading term of the weak asymptotic
of the solution of problem \thetag{2.1}, \thetag{2.2}
in the form \thetag{2.5},
where $S_i=S_i(\tau)$ and
$\phi_i=\varphi_{i0}(t)+\epsilon\varphi_{i1}(\tau)$
are the desired smooth functions.
We assume that
\begin{gather*}
S_i(\tau)\to0\quad\text{as }\tau\to\pm\infty,
\quad i=1,2,
\\
\varphi_{i1}(\tau)\to0\quad\text{as }\tau\to-\infty,
\quad
\varphi_{i1}(\tau)\to\varphi^\infty_{i1}\quad\text{as }\tau\to+\infty
\end{gather*}
as functions of exponential type.
Now we apply formulas from Sec. 1 and calculate the weak asymptotic
\[
u^n =\epsilon a_n\sum^2_{i=1}K^{(n)}_{i0}\delta(x-\phi_i)
+\epsilon a_n\Big\{\sum^2_{i=1}K^{(n)}_{i1}+R_n\Big\}\delta(x-x^*)
+O_{\mathcal{D}'}(\epsilon^2),
\]
\begin{equation} \tag {2.26}
(\epsilon u_x)^2
=\epsilon a'_2\sum^2_{i=1}\beta^2_iK^{(2)}_{i0}\delta(x-\phi_i)
+\epsilon a'_2\Big\{\sum^2_{i=1}\beta^2_iK^{(2)}_{i1}+Q_2\Big\}\delta(x-x^*)
+O_{\mathcal{D}'}(\epsilon^2),
\end{equation}
\begin{align*}
\frac{\partial u}{\partial t}
&=
c\frac{\partial}{\partial\tau}\sum^2_{i=1}K^{(1)}_{i1}\delta(x-x^*)
-\epsilon\sum^2_{i=1}V_iK^{(1)}_{i0}\delta'(x-\phi_i)\\
&\quad -\epsilon c\frac{\partial}{\partial\tau}\Big\{\sum^2_{i=1}
K^{(1)}_{i}\varphi_{i1}+\frac{V_i}{c}\tau K^{(1)}_{i1}\Big\}
\delta'(x-x^*)+O_{\mathcal{D}'}(\epsilon^2),
\end{align*}
\begin{align}
\frac{\partial u^2}{\partial t}
&=c a_2\frac{\partial}{\partial\tau}
\Big\{\sum^2_{i=1}K^{(2)}_{i1}+R_2\Big\}\delta(x-x^*)
-\epsilon a_2\sum^2_{i=1}V_iK^{(2)}_{i0}\delta'(x-\phi_i)
\tag {2.27}\\
&\quad
-\epsilon c a_2\frac{\partial}{\partial\tau}\Big\{\sum^2_{i=1}
\Big(K^{(2)}_{i}\varphi_{i1}+\frac{V_i}{c}\tau K^{(2)}_{i1}\Big)
+R^{(1)}_{2}\Big\}\delta'(x-x^*)+O_{\mathcal{D}'}(\epsilon^2).
\nonumber
\end{align}
Here $K^{(n)}_{i}=g^n_i/\beta_i$,
$K^{(n)}_{i0}=A^n_i/\beta_i$,
$K^{(n)}_{i1}=K^{(n)}_{i}-K^{(n)}_{i0}$,
and the numbers $a_n$ for $n=2,3$ and $a'_2$
are calculated by formulas \thetag{2.18},
while $a_1=1$.
We also set $V_i=\varphi'_{i0t}(t^*)$,
$c=\beta_1(V_2-V_1)$, and
$$
R_n=\frac1{\beta_2 a_n}\int
\big(g_1\omega(\theta\eta+\sigma)+g_2\omega(\eta)\big)^n\,d\eta
-\sum^{2}_{i=1}K^{(n)}_i.
$$
Moreover, using the notation \thetag{2.16},
we obtain
\begin{equation} \tag {2.28}
R_2=2\frac{g_1g_2}{\beta_2}\lambda_0,\quad
R^{(1)}_2=2\frac{g_1g_2}{\beta_2}
\Big\{\big(\frac{V_2}{c}\tau+\varphi_{21}\big)\lambda_0
+\frac1{\beta_2}\lambda_1\Big\},
\quad
Q_2=2g_1g_2\beta_1\lambda_2.
\end{equation}
By definition, the weak asymptotic of the solution
must satisfy the relations
\begin{equation}
\frac{\partial u}{\partial t}+\frac{\partial u^2}{\partial x}=O_{\mathcal{D}'}(\epsilon^2),
\quad
\frac{\partial u^2}{\partial t}+\frac43\frac{\partial u^3}{\partial x}
-3\frac{\partial}{\partial x}(\epsilon u_x)^2=O_{\mathcal{D}'}(\epsilon^2).
\tag {2.29}
\end{equation}
We substitute expansions \thetag{2.26}, \thetag{2.27}
into \thetag{2.29} and collect the term containing
$\delta(x-x^*)$, $\epsilon\delta'(x-\phi_i)$,
and $\epsilon\delta'(x-x^*)$.
By equating the obtained relations with zero,
we arrive at the relations
\begin{gather}
\frac{d}{d\tau}\sum^{2}_{i=1}K^{(1)}_{i1}=0,
\tag {2.30}\\
\frac{d}{d\tau}\Big\{\sum^{2}_{i=1}K^{(2)}_{i1}+R_2\Big\}=0,
\tag {2.31}\\
\frac{d}{d\tau}\sum^{2}_{i=1}\Big\{K^{(1)}_{i}\varphi_{i1}
+\tau\frac{V_i}{c}K^{(1)}_{i1}\Big\}=0,
\tag {2.32}\\
\begin{aligned}
&c\frac{d}{d\tau}\Big\{\sum^{2}_{i=1}
\Big(K^{(2)}_{i}\varphi_{i1}+\tau\frac{V_i}{c}K^{(2)}_{i1}\Big)
+R^{(1)}_2\Big\}\\
&=\frac43\frac{a_3}{a_2}\Big\{\sum^{2}_{i=1}K^{(3)}_{i1}+R_3\Big\}
-3\frac{a'_2}{a_2}\Big\{\sum^{2}_{i=1}\beta^2_{i}K^{(2)}_{i1}+Q_2\Big\},
\end{aligned} \tag {2.33}
\end{gather}
as well as the relations
\begin{equation}
-V_i K^{(1)}_{i0}+a_2K^{(2)}_{i0}=0,\;
-a_2 V_i K^{(2)}_{i0}+\frac43 a_3 K^{(3)}_{i0}
-3a'_2\beta^2_i K^{(2)}_{i0}=0,\; i=1,2,
\tag {2.34}
\end{equation}
describing the motion of solitary waves without
taking into account their interaction.
Relations \thetag{2.34} readily imply
\begin{equation}
V_i=a_2A_i,\quad \beta^2_i=\gamma A_i,\quad
\gamma=(\tfrac43 a_3-a^2_2)/3a'_2\,.
\tag {2.35}
\end{equation}

Let us transform system \thetag{2.30}--\thetag{2.33}.
Integrating the first three relations
and performing simple transformations,
we obtain
\begin{gather}
K^{(1)}_{11}+K^{(1)}_{21}=0,
\tag {2.36}\\
\sum^2_{i=1}K^{(2)}_{i1}+R_2=0,
\tag {2.37}\\
\frac{\sigma}{\beta_1}K^{(1)}_{11}
-\sum^2_{i=1}K^{(1)}_{i0}\varphi_{i1}=0.
\tag {2.38}
\end{gather}
 Next, we use \thetag{2.28}, \thetag{2.36} and rewrite
\thetag{2.37} as follows:
\begin{equation}
-2(A_2-A_1)K^{(1)}_{11}
+(\beta_1+\beta_2)(K^{(1)}_{11})^2
+\frac{2}{\beta_2}(A_1+\beta_1 K^{(1)}_{11})
(A_2-\beta_2 K^{(1)}_{11})\lambda_0\!=\!0.
\tag {2.39}
\end{equation}
Solving \thetag{2.39} for $K^{(1)}_{11}$
and choosing the root from the condition
that $K^{(1)}_{11}\to0$ as $\tau\to-\infty$,
we obtain
\begin{equation}
K^{(1)}_{11}=\frac{\beta_2}{\gamma}k(\sigma),
\tag {2.40}
\end{equation}
where the function $k(\sigma)$ is defined in \thetag{2.19}
for $\theta\ne1/2$ and in \thetag{2.20} for $\theta=1/2$.
Let us analyze formulas (2.19), (2.20).
The obvious inequality
$\lambda_0(\sigma,\theta)\leq 1/\sqrt{\theta}$ for $\theta<1$
implies
$$
1+\theta-2\lambda_0(\sigma,\theta)\theta\geq(1-\sqrt{\theta})^2\geq\mathop{\rm const}>0.
$$
Next, let us consider the expression $I$ defined in \thetag{2.19}.
It is easy to see that the maximal value
(with respect to $\sigma$) of $\lambda_0$,
which we denote by $\lambda_0(0,\theta)$,
is a function of $\theta\in[0,1]$ monotonically decreasing
from $3/2$ to $1$.
Therefore,
for $\theta\geq2/5$ there exists $\sigma^*=\sigma^*(\theta)$
for which the first term in the formula for $I$ is zero.
Nevertheless,
the second term in this formula is positive for $\theta<1/2$
and negative for $\theta>1/2$.
Hence for $\theta<1/2$ we have the following estimate
uniform with respect to $\sigma$:
\begin{equation}
I\geq\mathop{\rm const}>0,
\tag {2.41}
\end{equation}
and we see that for $\theta>1/2$ this inequality does not hold.

For $\theta=1/2$ we define $k$ by formula \thetag{2.20}.
Since $3-2\lambda_0\big|_{\theta=1/2}\geq\mathop{\rm const}>0$,
this formula determines the smooth function $k(\sigma)$.

Finally, for  $\theta\in(c_1,1/2]$, $c_1>0$,
the function $k(\sigma)$ is smooth.
It is also easy to see that this function is positive,
even with respect to $\sigma$, and tends to zero
as $\sigma\to\pm\infty$ as an exponential function.

Now let us consider relation \thetag{2.38}
together with the relation
\begin{equation}
\sigma=\tau+\beta_1(\varphi_{21}-\varphi_{11}).
\tag {2.42}
\end{equation}
By solving \thetag{2.38}, \thetag{2.42} for $\varphi_{i1}$,
we obtain
\begin{equation}
\begin{gathered}
\varphi_{11}=(\beta_1\sum^{2}_{i=1}K^{(1)}_{i0})^{-1}
(\tau K^{(1)}_{20}-\sigma K^{(1)}_{2}),\\
\varphi_{21}=(\beta_1\sum^{2}_{i=1}K^{(1)}_{i0})^{-1}
(-\tau K^{(1)}_{10}+\sigma K^{(1)}_{1}).
\end{gathered}
\tag {2.43}
\end{equation}
Let us transform  \thetag{2.33}.
After simple calculations,
taking into account \thetag{2.28}, \thetag{2.36},
 \thetag{2.40}, and \thetag{2.43},
we obtain
$$
\sum^2_{i=1}\Big(K^{(2)}_i\varphi_{i1}
+\tau\frac{V_i}{c}K^{(2)}_{i1}\Big)+R^{(1)}_2
=\frac{\beta^2_2}{\gamma^2}
\{-\tau(1-\theta)+\sigma(1-\theta)+\mathcal{J}(\sigma)\},
$$
where
\begin{equation}
\mathcal{J}(\sigma)=\sigma \big(\theta^{-1}(1-\theta-\theta^2)k-k^2\big)
+2\theta(\theta+k)(1-k)\lambda_1.
\tag {2.44}
\end{equation}
Now we note that the following relations hold:
\begin{gather*}
\sum^{2}_{i=1}K^{(3)}_{i1}+R_3
=\frac{\beta^5_2}{\gamma^3}L_1(\sigma),
\quad
\sum^{2}_{i=1}\beta^2_iK^{(2)}_{i1}+Q_2
=\frac{\beta^5_2}{\gamma^2}L_2(\sigma),
\\
c=a_2\frac{\beta^3_2}{\gamma}\theta(1-\theta^2),\quad
3\gamma a'_2=\frac43a_3-a^2_2,
\end{gather*}
where $L_1,L_2$ are defined in \thetag{2.17}.
Hence we can rewrite Eq. \thetag{2.33} in the autonomous form
\begin{equation}
(1-\theta+\mathcal{J}')\frac{d\sigma}{d\tau}
=1-\theta+\frac{F_1(\sigma)}{\theta(1-\theta^2)},
\tag {2.45}
\end{equation}
where
$$
F_1=\frac43\frac{a_3}{a^2_2}(L_1-L_2)+L_2,
\quad
\mathcal{J}'=\frac{\partial \mathcal{J}}{\partial\sigma}.
$$
Formulas  \thetag{2.16}, \thetag{2.44}
and the properties of $k(\sigma)$
readily imply that $\mathcal{J}(\sigma)\in C^\infty(\mathbb{R}^1)$
for $0<\theta\leq 1/2$.
Moreover, $\mathcal{J}(\sigma)$ is odd and tends to zero
as $\sigma\to\pm\infty$ as an exponential function.

Consider $\mathcal{J}(\sigma)$ for small $\theta$.
Expanding $\omega(\theta\eta+\sigma)$ in the Taylor series
around the point $\theta=0$, we obtain
$$
\lambda_0=\frac1{a_2}\omega(\sigma)+O(\theta^2),\quad
\lambda_1=O(\theta).
$$
Thus we have
\begin{equation}
k(\sigma)=\frac{\theta^2}{a_2}\omega(\sigma)+O(\theta^3),\quad
\mathcal{J}(\sigma)=\frac{\theta(1-\theta)}{a_2}\sigma\omega(\sigma)+O(\theta^3),
\tag {2.46}
\end{equation}
and there exists $\theta^*_1\in(0,1/2]$ such that uniformly in $\sigma$
\begin{equation}
1-\theta+\mathcal{J}'(\sigma)\geq\mathop{\rm const}>0\quad\text{for}\quad\theta<\theta^*_1.
\tag {2.47}
\end{equation}
It follows from this relation and the fact that $F_1(\sigma)$ is
smooth that Eq. \thetag{2.45} has a solution for all
$\tau\in\mathbb{R}^1$.
Next, for small $\theta$ we find
$$
\lambda_2=O(\theta),\quad
L_1=-\frac32\theta^2\omega(\sigma) +O(\theta^3),\quad
L_2=-6\theta^2\omega(\sigma)+O(\theta^4).
$$
Now we readily calculate that for small $\theta$
\begin{equation}
1-\theta+\frac{F_1(\sigma)}{\theta(1-\theta^2)}
=1-\theta+\frac65\theta\omega(\sigma)+O(\theta^2)>0.
\tag {2.48}
\end{equation}
Hence there exists $\theta^*_2\in(0,1/2]$ such that uniformly in $\sigma$
\begin{equation}
1-\theta+\frac{F_1(\sigma)}{\theta(1-\theta^2)}\geq\mathop{\rm const}>0
\quad\text{for}\quad\theta<\theta^*_2.
\tag {2.49}
\end{equation}
Inequalities \thetag{2.47} and \thetag{2.49}
and the fact that the function $F_1(\sigma)$ is even in $\sigma$
imply that $\sigma(\tau)$ tends to $\pm\infty$
as $\tau\to\pm\infty$ and is odd in $\tau$
up to an additive constant.

We set $\sigma_0=\sigma-\tau$ and
$\theta^*=\min\{\theta^*_1,\theta^*_2\}$.
As was mentioned above,
we cannot rigorously prove a more precise estimate for $\theta^*$,
but computer calculations give $\theta^*=1/2$.

Obviously, it follows from the above that problem \thetag{2.22}
has a solution smooth in $\tau$ and odd
up to an additive constant $\sigma_0(0)$.
The exponential stabilization rate
of $\sigma_0$ as $\tau\to\infty$
is a direct consequence of the fact
that the function $\mathcal{F}(\sigma)$ decreases
as an exponential function as $\tau\to\infty$.

To complete the proof of Theorem \ref{thm2.3},
it remains to note
that formulas \thetag{2.21}, \thetag{2.22}
readily follow from \thetag{2.33}, \thetag{2.37},
and \thetag{2.40}.
The third relation in \thetag{2.7} for functions
of the form \thetag{2.5} is obvious.
\end{proof}


\begin{proof}{Proof of Theorem \ref{thm2.2}}
Let us consider \thetag{2.27}. Taking into account
the normalization of the function $\omega$, we obtain
\begin{align*}
\sum^2_{i=1}
K^{(1)}_{i1}
&=\sum^2_{i=1}\frac{g_i-A_i}{\beta_i}\int\omega(\eta)\,d\eta
=\sum^2_{i=1}g_i\int\omega(\beta_i z)\,dz+\mathop{\rm const}\\
&=\frac1\epsilon\int\sum^2_{i=1}g_i
\omega\big(\beta_i\frac{x-\phi_i}{\epsilon}\big)\,dx+\mathop{\rm const}.
\end{align*}
Since $cd/d\tau=\epsilon d/dt$, we see that \thetag{2.27} is exactly
the conservation law \thetag{2.8} calculated for functions $u$
of the form \thetag{2.5}.
Consider relation \thetag{2.30}. With the help of the identity
$$
\frac{V_i}{c}\tau=\frac{\varphi_{i0}}{\epsilon}+q,\quad
q=\frac{V_1\varphi_{20}-V_2\varphi_{10}}{\epsilon(V_2-V_1)},
\quad i=1,2,
$$
by setting $\xi=\theta\eta+\sigma$, we calculate
\begin{align*}
&M := a_2 \sum^2_{i=1}
\big\{K^{(2)}_i\varphi_{i1}+\frac{V_i}{c}\tau K^{(2)}_{i1}\big\}
+a_2 R^{(1)}_2\\
&=\frac{g^2_1}{\beta_1}\frac{\phi_1}{\epsilon}\int\omega^2(\xi)\,d\xi
+\frac{g^2_2}{\beta_2}\frac{\phi_2}{\epsilon}\int\omega^2(\eta)\,d\eta
+2\frac{g_1g_2}{\beta_2}
\Big\{\big(\frac{\phi_2}{\epsilon}+q\big)
\int\omega(\eta)\omega(\xi)\,d\eta\\
&\quad
+\frac1{\beta_2}\int\eta\omega(\eta)\omega(\xi)\,d\eta\Big\}
+a_2\sum^2_{i=1}
\big(q K^{(2)}_{i1}-\frac{\varphi_{i0}}{\epsilon}K^{(2)}_{i0}\big).
\end{align*}
Performing the change $\eta=\beta_2(x-\phi_2)/\epsilon$,
we prove the relation
\begin{align*}
&\frac{\phi_2}{\epsilon}\int\omega(\eta)\omega(\xi)\,d\eta
+\frac1{\beta_2}\int\eta\omega(\eta)\omega(\xi)\,d\eta\\
&=\frac{\beta_2}{\epsilon^2}
\int x\omega\big(\beta_1\frac{x-\phi_1}{\epsilon}\big)
\omega\big(\beta_2\frac{x-\phi_2}{\epsilon}\big)\,dx.
\end{align*}
Now it is easy to see that even functions $\omega(\eta)$
satisfy the relation
\begin{align*}
M=&\frac1{\epsilon^2}\int x\Big(
g_1\omega\big(\beta_1\frac{x-\phi_1}{\epsilon}\big)
+g_2\omega\big(\beta_2\frac{x-\phi_2}{\epsilon}\big)\Big)^2\,dx\\
&+a_2 q\Big(\sum^2_{i=1} K^{(2)}_{i1}+R_2\Big)
-\frac{a_2}{\epsilon}\sum^2_{i=1}\varphi_{i0}K^{(2)}_{i0}.
\end{align*}
In a similar way, we obtain
\begin{gather*}
a_3\Big\{\sum^2_{i=1}K^{(3)}_{i1}+R_3\Big\}
=\frac1\epsilon\int\Big(\sum^2_{i=1}g_i
\omega\big(\beta_i\frac{x-\phi_i}{\epsilon}\big)\Big)^3\,dx
-a_3\sum^2_{i=1}K^{(3)}_{i0},
\\
a'_2\Big\{\sum^2_{i=1}\beta_iK^{(2)}_{i1}+Q_2\Big\}
=\frac1\epsilon\int\Big\{\epsilon\frac{\partial}{\partial x}
\Big(\sum^2_{i=1}g_i
\omega\big(\beta_i\frac{x-\phi_i}{\epsilon}\big)\Big)\Big\}^2\,dx\\
-a'_2\sum^2_{i=1}\beta^2_iK^{(2)}_{i0}.
\end{gather*}
We substitute the relations obtained into \thetag{2.33} and
take into account \thetag{2.29}.
Then \thetag{2.33} takes the form
\begin{equation} \tag {2.50}
\begin{aligned}
&\frac1\epsilon\Big\{\frac{d}{dt}\int xu^2\,dx
-\frac43\int u^3\,dx
+3\int(\epsilon u_x)^2\,dx\Big\}\\
&=\sum^2_{i=1}\Big\{a_2\frac{d}{dt}\varphi_{i0}K^{(2)}_{i0}
-\frac43 a_3 K^{(3)}_{i0}+3a'_2\beta^2_iK^{(2)}_{i0}\Big\}.
\end{aligned}
\end{equation}
However, we have $d\varphi_{i0}/dt=V_i$
and thus the right-hand in \thetag{2.50} is zero
by virtue of the second relation in \thetag{2.34}.
Hence Eq. \thetag{2.33} is exactly the energy
relation \thetag{2.11}  calculated for the function $u$
of the form \thetag{2.5}.

In a similar way we can prove that \thetag{2.31} and \thetag{2.32}
coincide with the equalities \thetag{2.9} and \thetag{2.10}
calculated for the function $u$ of the form \thetag{2.5}.

To complete the proof of Theorem \ref{thm2.2},
it remains to recall that relations \thetag{2.30} and \thetag{2.33}
are necessary and sufficient conditions for the function \thetag{2.5}
to satisfy  \thetag{2.7} with any test function $\psi$.
It should be noted
that relations \thetag{2.34} also coincide exactly
with the energy relations  \thetag{2.10} and \thetag{2.11}
whose left-hand sides, however,
are calculated for the solitary wave
$A_i\omega(\beta_i(x-\varphi_{i0})/\epsilon)$.
\end{proof}


\begin{proof}[Proof of Theorem \ref{thm2.1}]
Let $W_i(\eta)$ be some even functions from $C^\infty(\mathbb{R}^1)$
decreasing faster than $|\eta|^{-2}$ as $\eta\to\pm\infty$.
By $W_{i,\ell}(\eta)$ we denote their $2\ell+1$-order derivatives
with respect to $\eta$ for arbitrary $\ell\geq0$
and consider the function
\begin{equation}
U=\sum^2_{i=1}
\Big\{g_i(\tau)\omega\big(\beta_i\frac{x-\phi_i}{\epsilon}\big)
+A_iG_i(\tau)W_{i,\ell}
\big(\frac{\beta_i}{\nu_i(\theta)}\frac{x-\phi_i}{\epsilon}\big)\Big\},
\tag {2.51}
\end{equation}
where $\omega, g_i, \beta_i, \phi_i$ are the same as
in \thetag{2.5},
$G_i(\tau)\in C^\infty$ are arbitrary functions exponentially
decreasing as $\tau\to\pm\infty$,
and $\nu_i>0$ are sufficiently small variables depending
on the parameter $\theta$.

It is easy to prove that in the sense of $\mathcal{D}'$
the function \thetag{2.51} differs from the function \thetag{2.5}
by a value of order $O(\epsilon^{2\ell+2})$,
while their first derivatives differ by a value of order
$O(\epsilon^{2\ell+1})$. Nevertheless,
we have
\begin{gather*}
U^n=(u^n)_0+\epsilon a_n M_{(n)}(\tau,\nu,G)\delta(x-x^*)
+O_{\mathcal{D}'}(\epsilon^2), \\
(\epsilon U_x)^2=(\epsilon u_x)^2_0+\epsilon a'_2 P(\tau,\nu,G)\delta(x-x^*)
+O_{\mathcal{D}'}(\epsilon^2),\\
\begin{aligned}
\frac{\partial U^2}{\partial t}=&\big(\frac{\partial u^2}{\partial t}\big)_0
+a_2 C(\tau,\nu,G)\delta(x-x^*)\\
&+\epsilon a_2 C^{(1)}(\tau,\nu,G)\delta'(x-x^*)+O_{\mathcal{D}'}(\epsilon^2),
\end{aligned}
\end{gather*}
where $(u^n)_0$, $(\epsilon u_x)_0$, and $(\partial u^2/\partial t)_0$
are distributions defined in formulas \thetag{2.26}, \thetag{2.27},
the functions $M_{(n)}$, $P$, $C$, and $C^{(1)}$
tend to zero as $\tau\to\pm\infty$ as exponential functions,
and
$$
M_{(n)}=O(\nu_1+\nu_2), \quad
P=O(1/\nu_1+1/\nu_2),\quad
C=O(\nu_1+\nu_2).
$$

We require that the function $U$ is the weak modulo
$O_{\mathcal{D}'}(\epsilon^2)$ asymptotic of the solution of
problem \thetag{2.1}, \thetag{2.2}.
Just as above,
we obtain a system of model equations from \thetag{2.7}.
Four of these equations coincide with \thetag{2.30},
 \thetag{2.32}, and \thetag{2.34}, and the other two
have the form
\begin{gather}
\frac{d}{d\tau}\Big\{\sum^2_{i=1}K^{(2)}_{i1}+R_2+M_{(2)}\Big\}=0,
\tag {2.52}\\
\begin{aligned}
&c\frac{d}{d\tau}\Big\{\sum^2_{i=1}
\big(K^{(2)}_{i}\varphi_{i1}+\tau\frac{V_i}{c}K^{(2)}_{i1}\big)
+R^{(1)}_2+C^{(1)}\Big\}\\
&=\frac43\frac{a_3}{a_2}
\Big\{\sum^2_{i=1}K^{(2)}_{i1}+R_3+M_{(3)}\Big\}
-3\frac{a'_2}{a_2}
\Big\{\sum^2_{i=1}\beta^2_iK^{(2)}_{i1}+Q_2+P\Big\}.
\end{aligned}  \tag {2.53}
\end{gather}
Repeating the above constructions and using the
notation \thetag{2.40},
we obtain the auxiliary function $k$ from \thetag{2.30}
and \thetag{2.52}.
However, in contrast to \thetag{2.19}, we now have
$k=k(\sigma,\theta,\nu,G)$.
Here it is important that we can choose parameters $\nu_i$
arbitrary small so that the corresponding quadratic equation
has a real root for $\theta<\theta^*_1$.
Further, formulas \thetag{2.52} and \thetag{2.42}
yield
\begin{gather*}
\varphi_{11}=\Big(\beta_1\sum^{2}_{i=1}K^{(1)}_{i0}\Big)^{-1}
(\tau K^{(1)}_{20}-\sigma K^{(1)}_{2}+D_1),\\
\varphi_{21}=\Big(\beta_1\sum^{2}_{i=1}K^{(1)}_{i0}\Big)^{-1}
(-\tau K^{(1)}_{10}+\sigma K^{(1)}_{1}+D_2),
\end{gather*}
which differ from \thetag{2.43} by an additional function
$D_i=D_i(\tau,\sigma,\nu,G)$ such that $D_i\to0$
as $\tau\to\pm\infty$.
Here it should be pointed out that,
since $M_{(2)}$ vanishes as $\tau\to\pm\infty$,
the symmetry law for the resulting phase shifts of
trajectories \thetag{2.12}
remains valid for all perturbations of the function $u$
contained in \thetag{2.51}.

Obviously,
we now can transform \thetag{2.53} to the equation
similar to \thetag{2.45}:
\begin{equation}
\frac{d}{d\tau}\big\{(\sigma-\tau)(1-\theta)+\mathcal{J}(\sigma)
+\mathcal{J}_1(\tau,\sigma,\nu,G)\big\}
=\frac{F_2}{\theta(1-\theta^2)},
\tag {2.54}
\end{equation}
where
$$
F_2=F_1(\sigma)+\widetilde{F}_1(\tau,\sigma,\nu,G),\quad
\widetilde{F}_1=\frac43\frac{a_3}{a^2_2}(M_{(3)}-P)+P,
$$
$\mathcal{J}$ and $F_1$ are the same functions as in \thetag{2.45},
and $\mathcal{J}_1$ is a smooth function tending to zero as
$\tau\to\pm\infty$.
Let us integrate \thetag{2.54} with respect to $\tau$
and let $\tau\to\infty$.
Since $\sigma=\tau+\sigma_0$, we obtain the relation
\begin{equation}
\sigma^\infty_0(1-\theta)
=\Big(\int^{\infty}_{-\infty}F_1\,d\tau
+\int^{\infty}_{-\infty}\widetilde{F}_1\,d\tau\Big)
\Big/\theta(1-\theta^2).
\tag {2.55}
\end{equation}
The first integral on the right-hand side of \thetag{2.55}
corresponds to the fixed value of $\sigma^\infty_0(\theta)$
obtained in the proof of Theorem \ref{thm2.3}.
We see that the perturbations of $u$
contained in \thetag{2.51}
yield a function of soliton type (with respect to $\tau$),
namely the function $\widetilde{F}_1=(1-4a_3/3a^2_2)P+O(\nu_1+\nu_2)$.
Hence, by changing $G_i$ and $\nu_i$,
we can make the limit $\sigma^\infty_0$ to be  an arbitrary
finite number.
This fact means the asymptotic non-uniqueness.
The proof of Theorem \ref{thm2.1} is complete.
\end{proof}

In conclusion, we analyze  \thetag{2.45} for small $\theta$.
As shown in \thetag{2.48},
$F_1(\sigma)=\frac65\theta^2\omega(\sigma)+O(\theta^3)$.
Hence we have $(1-\theta)\sigma^{\infty}_0=\frac65\theta+O(\theta^2)$.
Thus the trajectory of the first soliton is shifted
by $-\epsilon\beta_1\varphi^{\infty}_{11}=\frac65\epsilon\theta+O(\epsilon\theta^2)$,
while according to the exact formula
this shift is equal to $-2\epsilon\mu=\epsilon\theta+O(\epsilon\theta^2)$.
In view of the fact that the weak asymptotic \thetag{2.5} is
asymptotically non-unique, this coincidence is astonishingly good.


\section{$\epsilon$-$\delta$-interaction in KdV type models}%3.

In this section we discuss a natural generalization of
the KdV equation
\begin{equation}
u_t+(u^m)_x+\epsilon u_{xxx}=0.
\tag {3.1}
\end{equation}
Under the assumption that $m\geq2$ is an integer, we study
the problem of interaction of two solitary waves.
Suppose that
\begin{equation}
u\big|_{t=0}=u^0(x,\epsilon),\quad
u^0=A_1\omega\big(\beta_1\frac{x-x^0_1}{\epsilon}\big)+
A_2\omega\big(\beta_2\frac{x-x^0_2}{\epsilon}\big),
\tag {3.2}
\end{equation}
where $\omega(\eta)$ is the exact solution of the model equation
(0.5) corresponding to \thetag{3.1}
(see Introduction)
and
\begin{equation}
\beta^2_i=\gamma^{m-1} A^{m-1}_i,\quad
\gamma=\Big(\frac{m-1}{m+3}\frac{a_m a_2}{a'_2}\Big)^{1/(m-1)}.
\tag {3.3}
\end{equation}
We assume that $A_i>0$ for even $m$
and that the sign of $A_i$ can be arbitrary for odd $m$.
In both cases, we assume that
$$
|A_2|>|A_1|,\quad x^0_1<x^0_1,
$$
which ensures the interaction of solitons at time $t^*$.

Precisely as in Sec. 2,
a weak modulo $O_{\mathcal{D}'}(\epsilon^2)$ asymptotic of the solution of
problem \thetag{3.1}  \thetag{3.2}
is a function $u$ satisfying for all $\psi=\psi(x)\in\mathcal{D}(\mathbb{R}^1)$
the relations
\begin{equation} \tag {3.4}
\begin{gathered}
\frac{d}{dt}\int u\psi\,dz-\int u^m\frac{\partial\psi}{\partial x}\,dx
=O(\epsilon^2),\\
\frac{d}{dt}\int u^2\psi\,dz-\frac{2m}{m+1}
\int u^{m+1}\frac{\partial\psi}{\partial x}\,dx
+3\int\Big(\epsilon\frac{\partial u}{\partial x}\Big)^2\frac{\partial\psi}{\partial x}\,dx
=O(\epsilon^2),\\
\int\big(u\big|_{t=0}-u^0\big)\psi\,dx
=O(\epsilon^2).
\end{gathered}
\end{equation}

The main result of this section is  the following
assertion.

\begin{theorem} \label{thm3.1}
There exists $\theta^*(m)\in(0,1)$ such that
for $\theta:= \beta_1/\beta_2\in\big(0,\theta^*(m)\big)$
the interaction of solitary waves
in the KdV type equation \thetag{3.1} follows the soliton
scenario.
The weak asymptotic of the solution of problem \thetag{3.1},
 \thetag{3.2} is asymptotically non-unique
in the terms of order $O_{\mathcal{D}'}(\epsilon^2)$.
\end{theorem}

The proof of this theorem differs from the proof of a similar statement in Sec. 2
only in small details.
Namely, we represent the leading term of the weak asymptotic
in the form
\begin{equation}
u=\sum^2_{i=1}(A_i+S_i(\tau))
\omega\Big(\beta_i\frac{x-\varphi_{i0}(t)-\epsilon\varphi_{i1}(\tau)}{\epsilon}\Big),
\tag {3.5}
\end{equation}
where $S_i$ and $\varphi_{i1}$ are assumed
to satisfy the same properties as in Sec. 2,
while the function $\omega(\eta)$,
precisely as in  \thetag{3.2}, depends on the parameter $m$.
As in Sec. 2, we establish the following assertion.

\begin{theorem} \label{thm3.2}
The condition $\theta\in(0,\theta^*(m))$,
the conservation laws
$$
\frac{d}{dt}\int u\,dx=0,\quad \frac{d}{dt}\int u^2\,dx=0,
$$
and the energy relations
\begin{gather*}
\frac{d}{dt}\int xu\,dx-\int u^m\,dx=0,\\
\frac{d}{dt}\int xu^2\,dx-\frac{2m}{m+1}\int u^{m+1}\,dx
+3\int\Big(\epsilon\frac{\partial u}{\partial x}\Big)^2\,dx=0
\end{gather*}
are necessary and sufficient
for the function \thetag{3.5} to be a weak modulo $O_{\mathcal{D}'}(\epsilon^2)$
asymptotic of the solution of problem  \thetag{3.1}, \thetag{3.2}.
\end{theorem}

The key point in the proof of these statements
is the derivation of model equations for $S_i$, $\phi_i$
and their analysis. To state the result, we introduce some
auxiliary notation and functions. We write
$$
s_i=\mathop{\rm sign} A_i,\quad s=s_1 s_2,\quad \kappa=\frac2{m-1},
$$
and the convolutions $\lambda_l=\lambda_l(\sigma)$, $l=0,1,2$,
are calculated by formulas of the form \thetag{2.16}.
By $\theta^*_0\in(0,1)$ we denote the number determined
by the relation
$$
1-\theta^*_0=2(\theta^*_0)^\kappa,
$$
and introduce a function $k=k(\sigma)$
determined in the case $s=1$ and $\theta=\theta^*_0$ by the formula
\begin{equation}
k=\frac{\lambda_0(1-\theta)}{1+\theta-2\lambda_0\theta}\Big|_{\theta=\theta^*_0},
\tag {3.6a}
\end{equation}
and in all other cases by the formula
\begin{equation}
k=\frac{1-\theta\lambda_0+s(\lambda_0-1)\theta^\kappa-\sqrt{I}}
{1+\theta-2\theta\lambda_0},
\tag {3.6}
\end{equation}
where
\begin{equation}
I=\big(1-s\theta^\kappa-\lambda_0(\theta^\kappa+\theta)\big)^2
+2\lambda_0\theta^\kappa\big(1-s\theta-\theta^\kappa(1+s)\big).
\tag {3.7}
\end{equation}
Further, we introduce the convolutions
$\lambda_{n,j}=\lambda_{n,j}(\sigma)$:
\begin{equation}
\lambda_{n,j}=\frac1{a_n}\int\omega^{n-j}(\theta\eta+\sigma)
\omega^j(\eta)\,d\eta
\tag {3.8}
\end{equation}
and the functions $L_{i,n}=L_{i,n}(\sigma)$
\begin{gather}
\begin{aligned}
L_{1,n}=&\theta^{n-1}\big[(k+s\theta^{\kappa-1})^n
-(s\theta^{\kappa-1})^n\big]+(1-k)^n-1 \\
&+\sum^{n-1}_{j=1}c^j_n\theta^{n-j}
(k+s\theta^{\kappa-1})^{n-j}(1-k)^j\lambda_{n,j},
\end{aligned} \tag {3.9}\\
L_{2,n}=(1+\theta^3) k^2-2k(1-s\theta^{\kappa+2})
+2\theta^2(k+s\theta^{\kappa-1})(1-k)\lambda_2,
\tag {3.10}
\end{gather}
where $c^j_n$ are binomial coefficients.
Now for all $m\geq2$ and any sign of $s$
except the case
\begin{equation}
s=-1,\quad m=3,
\tag {3.11}
\end{equation}
we defined a function $\sigma_0=\sigma_0(\tau)$ as a solution
of the problem
\begin{equation}
\frac{d\sigma_0}{d\tau}
=\frac{\mathcal{F}_m(\tau+\sigma_0)-\mathcal{J}'_m(\tau+\sigma_0)}
{r_m+\mathcal{J}'_m(\tau+\sigma_0)},\quad
\sigma_0\to0\quad\text{as }\tau\to-\infty.
\tag {3.12}
\end{equation}
Here $\mathcal{J}'_m(\sigma)=\partial\mathcal{J}_m(\sigma)/\partial\sigma$,
\begin{gather}
\mathcal{J}_m(\sigma) =\sigma k
\Big[\frac{1-2s\theta^\kappa-\theta^{2\kappa-1}}{1+s\theta^{\kappa-1}}-\theta k\Big]
+2\theta(\theta k+s\theta^{\kappa})(1-k)\lambda_1,
\tag {3.13}\\
r_m=\theta^{\kappa-1}\frac{s-\theta^\kappa}{1+s\theta^{\kappa-1}},
\tag {3.14}\\
\mathcal{F}_m(\sigma)=(1-\theta^2)^{-1}
\Big\{\frac{2m}{m+1}\cdot\frac{a_{m+1}}{a_2a_m}
(L_{1,m+1}-L_{2,m})+L_{2,m}\Big\}
-(1-r_m)L_{1,m}. \tag {3.15}
\end{gather}
In the exceptional case \thetag{3.11}
we find the function $\sigma_0=\sigma_0(\tau)$ as the solution
of the problem
\begin{equation}
\frac{d\sigma_0}{d\tau}
=\frac{L^-_{1,3}(\tau+\sigma_0)-\mathcal{J}'_3(\tau+\sigma_0)}
{1+\mathcal{J}'_3(\tau+\sigma_0)},\quad
\sigma_0\to0\quad\text{as }\tau\to-\infty,
\tag {3.16}
\end{equation}
where $\mathcal{J}'_3(\sigma)=\partial\mathcal{J}_3(\sigma)/\partial\sigma$,
$\mathcal{J}_3(\sigma)=-\sigma k(\sigma)$,
and $L^-_{1,3}=L_{1,3}\big|_{s=-1}$.

Using the functions $\sigma_0(\tau)$
and setting $\sigma=\tau+\sigma_0(\tau)$,
we define $\psi_1=\psi_1(\tau)$ in all the cases
except the case \thetag{3.11}
by the formula
\begin{equation}
\begin{aligned}
\psi_1(\tau)
=&(1+s\theta^{\kappa-1})^{-1}
\Big\{(\sigma_0+\tau)k(\sigma_0+\tau)-\sigma_0 &
&+\int^{\tau}_{-\infty} L_{1,m}
\big(\tau'+\sigma_0(\tau')\big)\,d\tau'\Big\}.
\end{aligned}\tag {3.17}
\end{equation}
In the exceptional case \thetag{3.11}
we define
\begin{equation}
\psi_1(\tau) =(1+\theta)^{-1}
\Big\{\theta\sigma(k^2-2k)-\sigma_0+2\theta^2(1-k)^2\lambda_1
+\int^{\tau}_{-\infty}
\big(\mathcal{F}_3+(1-r_3)L^-_{1,3}\big)\,d\tau'\Big\}.
\tag {3.18}
\end{equation}

\begin{theorem} \label{thm3.3}
There exists a number $0<\theta^*<1$ depending on $m$ and $s$
such that for $0<c_1\leq\theta\leq \theta^*$
problems \thetag{3.12} and \thetag{3.15}
have unique solutions.
For such $\theta$ the weak asymptotic of the solution
of problem \thetag{3.1}, \thetag{3.2}
has the form \thetag{3.5}, where
\begin{equation}
\begin{gathered}
S_1=\theta A_2 k(\tau+\sigma_0),\quad
S_2=- A_2 k(\tau+\sigma_0),\\
\varphi_{i0}=V_it+x^0_i,\quad
V_i=a_m A^{m-1}_i,\quad
\tau=\beta_1\frac{\varphi_{20}(t)-\varphi_{10}(t)}{\epsilon},\\
\psi_2(\tau)=\big(\psi_1(\tau)+\sigma_0(\tau)\big)/\theta,\quad
\varphi_{i1}=\psi_i/\beta_i,\quad i=1,2.
\end{gathered}
\tag {3.19}
\end{equation}
\end{theorem}

\begin{proof}{Proof of Theorem \ref{thm3.3}}
We represent the weak asymptotic of the solution
in the form \thetag{3.5} and apply
formulas from Sec. 1.
By substituting these expansions into \thetag{3.4}
and equating the coefficients of
$\delta'(x-x^*)$, $\epsilon\delta'(x-\phi_i)$,  and
$\epsilon\delta'(x-x^*)$  with zero,
we obtain the following system of equations
similar to system \thetag{2.30}-- \thetag{2.34}:
\begin{gather}
\frac{d}{d\tau}\sum^2_{i=1}K^{(1)}_{i1}=0,
\tag {3.20}\\
\frac{d}{d\tau}\Big\{\sum^2_{i=1}K^{(2)}_{i1}+R_2\Big\}=0,
\tag {3.21}\\
c\frac{d}{d\tau}\sum^2_{i=1}\Big\{K^{(1)}_{i}\varphi_{i1}
+\frac{V_i}{c}\tau K^{(1)}_{i1}\Big\}=f_m,
\tag {3.22}\\
c\frac{d}{d\tau}\Big\{\sum^2_{i=1}\Big(K^{(2)}_{i}\varphi_{i1}
+\frac{V_i}{c}\tau K^{(2)}_{i1}\Big)+R^{(1)}_2\Big\}=F_m,
\tag {3.23}\\
V_i K^{(1)}_{i0}=a_m K^{(m)}_{i0},
\tag {3.24}\\
V_i K^{(2)}_{i0}=\frac{2m}{m+1}\frac{a_{m+1}}{a_2} K^{(m+1)}_{i0}
-3\frac{a'_2}{a_2}\beta^2_i K^{(2)}_{i0},
\tag {3.25}
\end{gather}
where we use the notation
\begin{gather}
f_m=a_m\Big(\sum^2_{i=1}K^{(m)}_{i1}+R_m\Big),
\tag {3.26}\\
F_m=\frac{2m}{m+1}\frac1{a_2}f_{m+1}
-3\frac{a'_2}{a_2}\Big(\sum^2_{i=1}\beta^2_i K^{(2)}_{i1}+Q_2\Big).
\tag {3.27}
\end{gather}
Clearly,
expressions \thetag{3.24}, \thetag{3.25}
imply the same relations between the amplitude $A_i$,
the velocity $V_i$, and the coefficient $\beta_i$
as in the problem of motion of a single solitary wave:
\begin{equation}
V_i=a_m A^{m-1}_i,\quad
\beta^2_i=\gamma^{m-1}A^{m-1}_i,\quad
\gamma^{m-1}=\frac{(2m/(m+1))a_{m+1}-a_2a_m}{3a'_2}.
\tag {3.28}
\end{equation}
Since $\omega(\eta)$ is the solution of (0.5),
by using the properties of the coefficients $a_1$ and $a'_2$,
we can rewrite the coefficient $\gamma$
in the form \thetag{3.3}.

Let us consider \thetag{3.20} and \thetag{3.21}.
Integrating by $\tau$, we arrive at an equation of the
form \thetag{2.36}.
We set
\begin{equation}
A_i=s_i\frac{\beta^\kappa_i}{\gamma},\quad
K^{(1)}_{11}=s_2\frac{\beta^{\kappa-1}_2}{\gamma}k,\quad
\kappa=\frac{2}{m-1}.
\tag {3.29}
\end{equation}
Then, solving this equation for $k=k(\sigma)$,
we obtain formulas \thetag{3.6}, \thetag{3.7}.
Since $\lambda_0\leq1/\sqrt{\theta}$,
we have $1+\theta-2\theta\lambda_0\geq\mathop{\rm const}>0$ for $\theta<1$.
Next, relation \thetag{3.7} readily implies that
for $s=-1$
the inequality $I\geq\mathop{\rm const}>0$ holds uniformly
in $\sigma\in\mathbb{R}^1$ and $\theta\in(0,1)$.
Let us consider the case $s=1$.
Obviously,
the second term in \thetag{3.7} is positive only for
$\theta<\theta^*_0$.
For $\theta=\theta^*_0$ the first term
takes the form $(1+\theta)^2(1-\lambda_0)^2/4$.
Since $\lambda_0(0)$ is a monotonically decreasing function of $\theta$
and $\lambda_0(0)\big|_{\theta=1}=1$,
for all $\theta\geq\theta^*_0$ there exists
a value $\sigma=\sigma^*(m,\theta)$ such that
the first term is zero at this point.
Thus we obtain the condition $\theta<\theta^*_0$
for the estimate
$$
I\geq\mathop{\rm const}>0
$$
to be uniform with respect to $\tau$.

Finally, the function $k$ defined in \thetag{3.6a}, \thetag{3.6}
is real and smooth only if
$$
\theta\leq \theta^*_1,\quad
\theta^*_1=\theta^*_0\,\,\text{ for }s=1
\quad\text{and}\quad
\theta^*_1=1-\mu,\quad\mu>0,\quad\text{ for }s=-1.
$$

Let us consider Eqs. \thetag{3.22} and \thetag{3.23}.
Using formulas \thetag{3.9}, \thetag{3.10}, and \thetag{3.29},
we transform these equations to the form
\begin{gather}
\frac{d}{d\tau}
\big\{-\sigma k+s\theta^{\kappa-1}\psi_1+\theta\psi_2\big\}=L_{1,m},
\tag {3.30}\\
\begin{aligned}
&(1-\theta^2)\frac{d}{d\tau}
\big\{-\sigma k(2s\theta^\kappa+\theta k)+\theta^{2\kappa-1}\psi_1+\theta\psi_2
+2\theta(s\theta^\kappa+\theta k)(1-k)\lambda_1\big\} \\
&=\frac{2m}{m+1}\frac{a_{m+1}}{a_2a_m}(L_{1,m+1}-L_{2,m})+L_{2,m}.
\end{aligned}\tag {3.31}
\end{gather}
Here we have used the notation
\begin{equation}
\psi_i=\beta_i\varphi_{i1},\quad i=1,2.
\tag {3.32}
\end{equation}
We supplement relations \thetag{3.30} and  \thetag{3.31}
by the identity
\begin{equation}
\sigma=\tau+\sigma_0,\quad
\sigma_0=\theta\psi_2-\psi_1,
\tag {3.33}
\end{equation}
and consider this system in the exceptional case \thetag{3.11}.
Since $\kappa=1$ in this case, \thetag{3.30} takes the form
\begin{equation}
\{1-(\sigma k)'\}\frac{d\sigma}{d\tau}=1+L^-_{1,3}(\sigma),
\tag {3.34}
\end{equation}
where the prime denotes the derivative with respect to $\sigma$
and the superscript $-$ means that we set $s=-1$ in the formula
for $L_{1,3}(\sigma)$.
Now by using the explicit form of the function $k$, we can
verify that
$$
1-\frac{\sqrt{1+\theta}}{1-\sqrt{\theta}}\leq k\leq 0,
\quad
1-\frac{c\theta }{(1-\sqrt{\theta})^2}\leq 1-(\sigma k)'
\leq \frac{\sqrt{1+\theta}}{1-\sqrt{\theta}},
$$
where the constant $c$ independent of $\theta$
is the maximum value of the function $-\sigma\lambda'_0(\sigma)$.

The above formulas readily imply that there exists $\theta^*_2>0$
such that for $0<\theta\leq \theta^*_2$
equation \thetag{3.34}
can be solved for all $\tau\in \mathbb{R}^1$.
Obviously, problem \thetag{3.16} follows from
relations \thetag{3.33} and \thetag{3.34}
and the condition $\varphi_{i1}\to0$ as $\tau\to-\infty$.

Let us study the general case.
We set $\theta\psi_2=\psi_1+\sigma-\tau$ and,
by using \thetag{3.30},
define $d\psi_1/d\tau$ as a function of $\sigma$ and $\tau$.
Then \thetag{3.31} can be transformed to the form
\begin{equation}
\big(r_m+\mathcal{J}'_m(\sigma)\big)\frac{d\sigma}{d\tau}
=\mathcal{F}_{m}(\sigma)+r_m,
\tag {3.35}
\end{equation}
where the notation \thetag{3.13}, \thetag{3.15}
are used.
Let us prove that \thetag{3.35} is globally solvable
for, at least, sufficiently small $\theta$.
For this purpose, we first note that for $\theta\ll1$,
\begin{equation}
\begin{gathered}
\lambda_0=\lambda^0_0+O(\theta^2),\quad \lambda^0_0=\frac1{a_2}\omega(\sigma),\\
\lambda_1=\theta\lambda^0_1+O(\theta^3),\quad \lambda^0_1=\frac1{a_2}\omega'(\sigma)
\int\eta^2\omega(\eta)\,d\eta,
\end{gathered}
\tag {3.36}
\end{equation}
Therefore,
$$
k=s_1\theta^\kappa\lambda^0_0+O(\theta^{\kappa+1}),
$$
and we calculate
\begin{gather*}
\mathcal{J}_3=\frac12\sigma\lambda^0_0\theta+O(\theta^2),\quad
\mathcal{J}_4=\sigma\lambda^0_0\theta+O(\theta^{4/3}),\\
\mathcal{J}_5=-2s_1\sigma\lambda^0_0\theta^{3/2}+O(\theta^2),\quad
\mathcal{J}_m=-s_2\sigma\lambda^0_0\theta^{2\kappa}+O(\theta)\quad\text{for }m\geq6.
\end{gather*}
At the same time, we have
$$
r_3=\frac12+O(\theta),\quad
r_4=1+O(\theta^{1/3}),\quad
r_m=1+O(\theta^\kappa)\quad\text{for } m\geq5.
$$
Hence there exists $\theta^*_m$ such that
$$
r_m+\mathcal{J}'_m(\sigma)\geq\mathop{\rm const}>0\quad\text{for}\quad \theta<\theta^*_m.
$$
It remains to note that $\mathcal{J}'_m(\sigma)$ and $\mathcal{F}_m(\sigma)$ are
bounded by a constant and tend to zero
both as $\sigma\to\pm\infty$ and $\theta\to0$.
This implies that \thetag{3.35} is solvable
for $\theta\leq\theta^*=\min\{\theta^*_1,\theta^*_m\}$,
and hence problem \thetag{3.12} is solvable.

The last step in deriving formulas \thetag{3.19}
is the calculation of the function $\psi_1$.
In the exceptional case \thetag{3.11}
formula \thetag{3.18} is obtained by integrating \thetag{3.33}
for $m=3$ and $s=-1$.
In the general case
we obtain \thetag{3.17} by integrating \thetag{3.32}.
The proof of Theorem \ref{thm3.3} is complete.
\end{proof}

The proof of Theorem \ref{thm3.2} and
of the asymptotic non-uniqueness of the weak asymptotic
is carried out in the same way as the proofs in Sec. 2.


\begin{thebibliography}{99}

%1
\bibitem{AblS}
M. J. Ablowitz, M. J.,  \& Segur, H.,
{\it Solitons and the Inverse Scattering Transform},
SIAM (Philadelphia, 1981).


%5
\bibitem{BonaSa}
Bona, J. L, \& Sant, J. C.,
Dispersive blow-up of solutions of generalized
Korteweg-de Vries equations,
J. Differential Equations {\bf 103}, 3--57 (1993).

%8
\bibitem{Dan}
Danilov, V. G.,
Generalized solutions describing singularity interaction,
International Journal of Math. and Mathematical Sci. {\bf 30},
1--14 (2002).

%6
\bibitem{DaSh1}
Danilov, V. G., \& Shelkovich, V. M.,
Propagation and interaction of shock waves of quasilinear
equation,  Nonlinear Studies {\bf 8}, 135--170 (2001).

%7
\bibitem{DaSh2}
Danilov, V. G., \& Shelkovich, V. M.,
Propagation and interaction of nonlinear waves to
quasilinear equations.
In: {\it Proc. Eighth Inter. Conference on Hyperbolic
Problems},  Birkhauser (2001).

%9
\bibitem{DaSh3}
Danilov, V. G., \& Shelkovich, V. M.,
Propagation of infinitely narrow $\delta$-solitons,
\newline
http://arXiv.org/abs/math-ph/0012002

%13
\bibitem{DaOmRa} Danilov, V. G.,  Omel$'$yanov, G. A., \& Radkevich, E. V.,
Hugoniot-type conditions and weak solutions to the
phase-field system,
European J. Appl. Math. \textbf{10}, 55--77 (1999).

%4
\bibitem{IlKal}
Il'in, A. M., \& Kalyakin, L. A.,
Perturbation of finite-soliton solutions to the
Korteweg-de Vries equation,
Dokl. Ross. Akad. Nauk {\bf 336}, 595--598 (1994);
English transl. in Phys. Dokl. {\bf 40}, 603--606 (1995).

%11
\bibitem{LaxLev1}
Lax, P. D., \& Levermore, C. D.,
The small dispersion limit of the KdV equation, I,
Comm. Pure and Appl. Math., 253--291 (1983).

%12
\bibitem{LaxLev2}
Lax, P. D., \& Levermore, C. D.,
The small dispersion limit of the KdV equation, II,
Comm. Pure and Appl. Math., 577--593, (1983).


%14
\bibitem{LinSc}
Linares, F., \& Scialom, M.,
On the smoothing properties of solutions to the modified
Korteweg-de Vries equations,
J. Differential Equations {\bf 106}, 141--154 (1993).

%2
\bibitem{MasOm1}
Maslov, V. P., \& Omel'yanov, G. A.,
Asymptotic soliton-form solutions of equations with small dispersion,
Uspekhi Mat. Nauk {\bf 36}, 63--126 (1981);
English transl. in Russian Math. Surveys {\bf 36}, 73--149 (1981).

%3
\bibitem{MasOm2}
Maslov, V. P., \& Omel'yanov, G. A.,
{\it  Geometric Asymptotics for Nonlinear PDE}, I,
AMS, MMONO 202 (Providence, RI, 2001).

%15
\bibitem{McSc}
McLaughlin, D. W., \& Scott, A. C.,
A multisoliton perturbation theory,
in: {\it Solitons in Action},
Academic Press, New York, 1978.
pp. 201--256.

%16
\bibitem{MolVak}
Molotkov, I. A., \& Vakulenko, S. A.,
{\it Concentrated Nonlinear Waves\/}  (in Russian),
Leningrad Univ (Le\-nin\-g\-rad, 1988).

%10
\bibitem{Omel}
Omel'yanov, G.A.,
Propagation of singularities: multidimensional effects
(KP and phase field type equations).
In: {\it Proc. I. G. Petrovskii Seminar}, MGU (Moscow, 2002).

\end{thebibliography}

\end{document}