\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 230, pp. 1--16.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}
\begin{document}
\title[\hfilneg EJDE-2016/230\hfil Integrability of the limit model]
{Integrability of the limit model for a vacuum
diode and a solution to the singular boundary-value problem}
\author[A. A. Kosov, E. I. Semenov, A. V. Sinitsyn \hfil EJDE-2016/230\hfilneg]
{Alexander A. Kosov, Edward I. Semenov, Alexander V. Sinitsyn}
\address{Alexander A. Kosov \newline
Matrosov Institute for System Dynamics and Control Theory,
Siberian Branch of Russian Academy of Sciences,
ISDCT SB RAS Lermontov str., 134, Post Box 292
664033, Irkutsk, Russia}
\email{aakosov@yandex.ru}
\address{Edward I. Semenov \newline
Matrosov Institute for System Dynamics and Control Theory,
Siberian Branch of Russian Academy of Sciences,
ISDCT SB RAS Lermontov str., 134, Post Box 292
664033, Irkutsk, Russia}
\email{edwseiz@gmail.com}
\address{Alexander V. Sinitsyn
Universidad Nacional de Colombia,
Departamento de Matematicas
Bogota, Colombia}
\email{avsinitsyn@yahoo.com}
\thanks{Submitted August 4, 2015. Published August 23, 2016.}
\subjclass[2010]{34B16, 34B60, 35Q83}
\keywords{First integrals; singular boundary value problem; vacuum diode}
\begin{abstract}
This article concerns singular boundary-value problems for a vacuum diode model.
We prove the integrability of a system of nonlinear differential
equations and construct a complete system of the first integrals; thus
developing a method for solving singular boundary value problems.
Also we study the asymptotic behavior of the solution in a neighborhood
of the singular point.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\allowdisplaybreaks
\section{Introduction}
Modeling a plasma as a flow of charged particles interacting in a vacuum
usually needs the application of the Vlasov-Maxwell
or Vlasov-Poisson equations \cite{bat93,cd2000,b2007,lxz2014}.
When solving these nonlinear systems of partial differential equations (PDEs)
with initial and boundary conditions, it is necessary to find the solutions
and ascertain their properties in a number of conditions
(positiveness, monotonicity, singularity, etc.).
We study a simpler model described by a system of ordinary differential equations
with boundary conditions that retains the principal physical properties
of the initial model, and is a more efficient way to overcome possible mathematical
difficulties. This approach allows us
to obtain a limit model for plane vacuum diode magnetic insulation.
This model is given by a system of two singular second-order ODEs \cite{ben98}.
Related results to this limit model and the boundary-value problem are found
in \cite{ben98,sin2001} where analytical and numerical methods are combined.
In this article, we concentrate on exact analytical methods
to integrate the corresponding nonlinear systems and transform
the singular boundary value problem into a system of nonlinear equations.
This article is structured as follows.
Section 2 presents the description of the model, the formulation of the
corresponding singular boundary-value problem, and the definition of a solution
taht differs from the classical solution.
Section 3 presents the mathematical model in Hamiltonian form.
Using Liouville theorem, we prove the system integrability and construct
a complete system of four first integrals. Here the Hamiltonian form
of equations allows us to apply the classic integration methods for
nonlinear systems developed in analytical mechanics \cite{w1927}.
Section 4 describes the method of constructing a solution of the singular
boundary value problem, which is based on the use of the complete system of
first integrals.
Section 5 presents another method of integrating nonlinear systems,
which is based on the replacement
of the variables represented in a special form proposed by the authors
in \cite{ss2010}.
This method reduces the problem of integration to the same quadratures
described in Section 3. The quadratures are
represented by some combinations of elementary and elliptic functions.
Section 6 describes a class of exact solutions for the nonlinear system in
explicit form. The cases, when the explicit solutions imply exact solutions of
non-singular boundary-value problems, and the cases, when a sequence of
explicit solutions approximate a solution of the initial singular
boundary-value problem (in some sense), are given.
Section 7 shows an asymptotic representation of the solution for the boundary
value problem in the vicinity of the singular point.
In such representations, in particular,
the electric potential approximation up to $O(x^{4/3})$, agree with
the estimates of the upper and lower solutions (obtained by a different method)
to the boundary-value problem formulated by the authors \cite{sin2001}.
Section 8 gives examples illustrating the numerical results obtained.
\section{Model and problem statement}
The limit model of a plane vacuum diode was proposed in \cite{ben98}.
The model consists of the two second-order nonlinear ordinary differential
equations:
\begin{equation} \label{e1}
\frac{d^2\varphi}{dx^2}=j\frac{(1+\varphi)}{\sqrt{(1+\varphi)^2-a^2-1}},
\quad
\frac{d^2a}{dx^2}=j\frac{a}{\sqrt{(1+\varphi)^2-a^2-1}},
\end{equation}
where the independent variable $x\in[0,1]$ denotes a relative distance from the
cathode ($x=1$ corresponds to the anode). The function $\varphi(x)$ describes
the distribution of the electric potential in the process of moving from the
cathode to anode; $a(x)$ is the potential of the magnetic field; the
model's parameter $j$ denotes the density of current through the diode.
System \eqref{e1} describes the electric and magnetic fields inside the diode,
and its solution shall satisfy the following boundary conditions:
\begin{gather}
\varphi(0)=0, \quad a(0)=0, \quad \varphi'(0)=\frac{d\varphi}{dx}(0)=0, \label{e2}\\
\varphi(1)=\varphi_1, \quad a(1)=a_1. \label{e3}
\end{gather}
The boundary-value problem \eqref{e1}--\eqref{e3} is singular: after substituting
conditions \eqref{e2} into equations \eqref{e1} for $x=0$, the denominator vanishes.
Therefore, the classical definition of the solution as a pair of functions
$(\varphi(x), a(x))$ satisfying \eqref{e2} and \eqref{e3} and converting
\eqref{e1} into an identity on an interval $x\in[0,1]$ (the derivatives
at the ends of the interval are considered as unilateral)
cannot be applied to this problem. So, it is
necessary to define the concept of a solution for \eqref{e1}--\eqref{e3}.
The parameter $j$ is free, see \cite{ben98}, and must be
found together with the solution of a boundary value problem
\eqref{e1}--\eqref{e3}.
Let $\Omega=\{(\varphi, a): (1+\varphi)^2-a^2-1>0\}$.
On each compact subset $\Omega$ the
right sides of \eqref{e1} have bounded partial derivatives. Therefore,
the conditions of existence and uniqueness of solutions for the initial-value
problem \eqref{e1} are satisfied. Furthermore, due to obvious symmetry we may confine
our consideration to investigation of only the solutions with positive values
of $1+\varphi(x)$, $a(x)$, i.e., consider the problem only in domain
$\Omega_{+}=\Omega\cap\{(\varphi,a): 1+\varphi>0, \; a>0\}$.
Let the conditions of the theorem on existence and uniqueness be satisfied for the
right end, i.e. $\Theta_1=(1+\varphi_1)^2-1-a_1^2>0$.
\begin{definition} \label{def1}\rm
A pair $(\varphi(x), a(x))$ of
twice differentiable functions in $]0, 1]$ assuming values on $\Omega_{+}$
is a solution of \eqref{e1}--\eqref{e3} when
\begin{itemize}
\item[(1)] $\varphi(1)=\varphi_1$, $a(1)=a_1$;
\item[(2)] on each interval $x\in[\epsilon, 1]$, $0<\epsilon<1$ -- after
the substitution -- $(\varphi, a)$ satisfies \eqref{e1};
\item[(3)] there are limits $\lim_{\epsilon\to +0}\varphi(\epsilon)=0$,
$\lim_{\epsilon\to +0}a(\epsilon)=0$, $\lim_{\epsilon\to +0}
\varphi'(\epsilon)=0$.
\end{itemize}
At $x=0$ this function is redefined by the first two relations in its
property (3).
\end{definition}
Properties of the above definition:
\begin{itemize}
\item[(a)] This definition does not necessitate substitution of the boundary
conditions at the interval's left end into the system, what allows us to avoid the
division by zero;
\item[(b)] It does not impose any restrictions on the behavior of the first
derivative $a'(x)$ and the second derivatives at the interval's left end;
\item[(c)] It may obviously be upgraded also for the case, when
$\Theta_1=(1+\varphi_1)^2-1-a_1^2=0$ and when there is a singularity
at the right end of the interval.
Furthermore, this definition allows for $\lim_{\epsilon\to +0}a'(\epsilon)$
not to exist or to be infinity.
\end{itemize}
The principal objective of the work initiated by the authors is to develop a
method for constructing a solution of the singular boundary-value problem
\eqref{e1}--\eqref{e3} in terms of Definition \ref{def1}. To this end it is necessary to show
that system \eqref{e1} is integrable in quadratures and to construct a complete
system of the first integrals. Furthermore, formulas in explicit form,
which approximate the solution of the boundary-value problem at $x=0$, are obtained.
\section{Representation of the problem in Hamiltonian form and its integrability}
Let
\begin{gather*}
t=x,\quad q_1=\varphi(x),\quad q_2=a(x),\quad p_1=-\varphi'(x), \\
p_2=a'(x),\quad q=col(q_1,q_2)\in R^2,\quad p=col(p_1,p_2)\in R^2, \\
H(q,p)=\frac{1}{2}(-p_1^2+p_2^2)+j\sqrt{(1+q_1)^2-q_2^2-1}.
\end{gather*}
Hence \eqref{e1} is equivalent to
\begin{equation} \label{e4}
\dot{q}=\frac{\partial H}{\partial p}, \quad
\dot{p}=-\frac{\partial H}{\partial q}.
\end{equation}
In terms of the new variables the boundary conditions may be rewritten as
\begin{gather}
q_1(0)=0, \quad q_2(0)=0, \quad p_1(0)=0, \label{e5} \\
q_1(1)=Q_1\equiv \varphi_1, \quad q_2(1)=Q_2\equiv a_1. \label{e6}
\end{gather}
The new boundary-value problem \eqref{e4}--\eqref{e6} is equivalent to the initial
problem \eqref{e1}--\eqref{e3}.
It differs from the initial system only by the Hamiltonian
form of the system of differential equations. The Hamiltonian form
of \eqref{e4}--\eqref{e6}
allows one to apply the integration technique developed for solving problems of
analytical mechanics \cite{w1927}.
The Hamiltonian system \eqref{e4} has the energy integral
\begin{equation}
J_1\equiv H(q,p)=\frac{1}{2}(-p_1^2+p_2^2)+j
\sqrt{(1+q_1)^2-q_2^2-1}=c_1=:{\rm const}. \label{e7}
\end{equation}
It can readily be seen that another first integral of system \eqref{e4} has the form
\begin{equation}
J_2=(1+q_1)p_2+q_2p_1=c_2=:{\rm const}. \label{e8}
\end{equation}
The first integrals $J_1$, $J_2$ do not depend explicitly on time.
The rank of the Jacobi matrix for the first integrals $J_1$, $J_2$ in
$\mathcal{D}=\{(q_1, q_2, p_1, p_2): (q_1, q_2)\in\Omega_{+}$,
$(p_1, p_2)\in\mathbb{R}^2\}$
becomes less than $2$ only for the set
\begin{align*}
\mathcal{M}=\Big\{&
(q_1, q_2, p_1, p_2): (q_1, q_2)\in\Omega_{+},\;
p_1=\pm\frac{\sqrt{j}q_2}{((1+q_1)^2-q_2^2-1)^{1/4}},\\
&p_2=\pm\frac{\sqrt{j}(1+q_1)}{((1+q_1)^2-q_2^2-1)^{1/4}}
\Big\}.
\end{align*}
This set $\mathcal{M}$ is two-dimensional, and it does not separate the
four-dimensional set $\mathcal{D}$ into sub-domains.
The first integrals $J_1$, $J_2$ are functionally independent in
$\mathcal{D}\setminus \mathcal{M}$. So, due to the Liouville theorem \cite{w1927},
the diode model \eqref{e4} is integrable in this domain.
To construct the two necessary first
integrals, let us express the momentum from \eqref{e7} and \eqref{e8} in terms of the
coordinates
\begin{gather}
p_1=-\frac{c_2q_2}{w^2}\mp\frac{\sqrt{c_2^2-2w^2(c_1-j\sqrt{w^2-1})}}{w^2}(1+q_1),
\label{e9} \\
p_2=\frac{c_2(1+q_1)}{w^2}\pm\frac{\sqrt{c_2^2-2w^2(c_1-j\sqrt{w^2-1})}}{w^2}q_2,
\label{e10}
\end{gather}
where $w^2=(1+q_1)^2-q_2^2$. These formulas may be written in the form
$p_1=\frac{\partial\Phi}{\partial q_1}$, $p_2=\frac{\partial\Phi}{\partial q_2}$,
where function $\Phi(q_1,q_2,c_1,c_2)$ is given by the formula
\begin{equation}
\begin{aligned}
&\Phi(q_1,q_2,c_1,c_2)\\
&=\frac{c_2}{2}\ln\big| \frac{1+q_1+q_2}{1+q_1-q_2}\big|
\pm
\int_{\sqrt{(1+q_1)^2-q_2^2}}^{\sqrt{(1+Q_1)^2-Q_2^2}}
\frac{\sqrt{c_2^2-2s^2(c_1-j\sqrt{s^2-1})}}{s}ds.
\end{aligned} \label{e11}
\end{equation}
The integral in \eqref{e11} may be reduced to elementary and elliptic functions.
By the Liouville theorem, the two other necessary first integrals $J_3$, $J_4$
are expressed via function $\Phi(q_1,q_2,c_1,q_2)$:
\begin{gather}
\begin{aligned}
J_3&\equiv\frac{\partial\Phi}{\partial c_2}\\
&=\frac{1}{2}\ln\big| \frac{1+q_1+q_2}{1+q_1-q_2}\big|
\pm \int_{\sqrt{(1+q_1)^2-q_2^2}}^{\sqrt{(1+Q_1)^2-Q_2^2}}
\frac{c_2ds}{s\sqrt{c_2^2-2s^2(c_1-j\sqrt{s^2-1})}}=c_3,
\end{aligned} \label{e12} \\
J_4\equiv\frac{\partial\Phi}{\partial c_1}
=\mp\int_{\sqrt{(1+q_1)^2-q_2^2}}^{\sqrt{(1+Q_1)^2-Q_2^2}}
\frac{sds}{\sqrt{c_2^2-2s^2(c_1-j\sqrt{s^2-1})}}=t+c_4, \label{e13}
\end{gather}
where $c_3$ and $c_4$ are constants.
Integrals in \eqref{e12} and \eqref{e13} are used only when the radicals
in the denominators are nonzero. Since the identity $(1+q_1)p_1+q_2p_2\equiv 0$
holds on the whole set $\mathcal{M}$, and equality
$$
(1+q_1)p_1+q_2p_2=\mp\sqrt{c_2^2-2w^2(c_1-j\sqrt{w^2-1})}
$$
follows from \eqref{e9}, \eqref{e10}, the nonzero radicals in the denominators
of \eqref{e12}, \eqref{e13} provide for independence of the first integrals
$J_1$ and $J_2$.
The system of four first integrals \eqref{e7}, \eqref{e8}, \eqref{e12}
and \eqref{e13} obtained may be used in domain $\Omega_{+}$ for the purpose
of solving initial and boundary-value problems, while including problem
\eqref{e4}--\eqref{e6}.
\section{Solving the singular boundary-value problem}
Taking $t=1$ in \eqref{e13} and applying conditions \eqref{e6} to the
right end of the interval, one obtains $c_4=-1$. Similarly, from \eqref{e12} we have
$ c_3=\frac{1}{2}\ln\big| \frac{1+Q_1+Q_2}{1+Q_1-Q_2}\big|$.
At $t=0$, from conditions \eqref{e2} and \eqref{e5}
for the right end of the interval and from
integrals \eqref{e7}, \eqref{e8} one obtains $p_1(0)=0$, $p_2(0)=c_2$, $c_2^2=2c_1$.
Thus, we introduce the following functions
\begin{gather*}
F(u, v, Q_1, Q_2)=\int_1^{\sqrt{(1+Q_1)^2-Q_2^2}}
\frac{sds}{\sqrt{u^2(1-s^2)+2vs^2\sqrt{s^2-1}}}, \\
G(u, v, Q_1, Q_2)=\int_1^{\sqrt{(1+Q_1)^2-Q_2^2}}
\frac{uds}{s\sqrt{u^2(1-s^2)+2vs^2\sqrt{s^2-1}}}.
\end{gather*}
These functions may be represented as combinations of elementary and elliptic
functions of their arguments. So, basing on \eqref{e12} and \eqref{e13},
by the relations between the arbitrary constants derived from the boundary conditions,
we have the following statement.
\begin{theorem} \label{thm1}
If $c_2=u_{*}$ and $j=v_{*}$ represent a solution of the system of two nonlinear
equations
\begin{equation}
F(u, v, Q_1, Q_2)=1, \quad
G(u, v, Q_1, Q_2)=\frac{1}{2}\ln\big|\frac{1+Q_1+Q_2}{1+Q_1-Q_2}\big|, \label{e14}
\end{equation}
then the solution of the boundary-value problem \eqref{e4}--\eqref{e6}
exists and represents a solution of the initial-value problem for \eqref{e4}
with the initial conditions at the right end of the interval assigned by
\begin{gather}
q_1(1)=Q_1, \quad
p_1(1)=-\frac{u_{*}Q_2}{Z^2}
-\frac{\sqrt{u_{*}^2(1-Z^2)+2v_{*}Z^2\sqrt{Z^2-1}}}{Z^2}(1+Q_1), \label{e15}\\
q_2(1)=Q_2,\quad
p_2(1)=\frac{u_{*}(1+Q_1)}{Z^2}
+\frac{\sqrt{u_{*}^2(1-Z^2)+2v_{*}Z^2\sqrt{Z^2-1}}}{Z^2}Q_2, \label{e16}
\end{gather}
from \eqref{e9} and \eqref{e10}, where $Z^2=(1+Q_1)^2-Q_2^2$.
\end{theorem}
If \eqref{e14} is inconsistent, this does not mean that the boundary-value
problem \eqref{e4}--\eqref{e6} does not possess any solution.
There are situations when equalities \eqref{e14} are violated.
This takes place when the signs change in the
expressions \eqref{e9} and \eqref{e10} describing the momentum expressed via the
coordinates in the solution of boundary value problem.
\section{Integration by replacement of the variables}
System \eqref{e1} has solutions of the form \cite{ss2010}:
\begin{equation}
\begin{gathered}
\varphi(x)=\frac{1}{2\gamma}z(x)[e^{\omega(x)}+\gamma^2 e^{-\omega(x)}]-1, \\
a(x)=\frac{1}{2\gamma}z(x)[e^{\omega(x)}-\gamma^2 e^{-\omega(x)}].
\end{gathered}\label{e17}
\end{equation}
Here, $z(x)$, $\omega(x)$ are functions to be defined; $\gamma\ne 0$ is an
arbitrary constant.
To this end we substitute the ansatz \eqref{e17} into \eqref{e1} and, after simple
transformations, we obtain the system
\begin{equation}
\begin{gathered}
A\big[\frac{1}{2\gamma}e^{\omega(x)}+\frac{\gamma}{2}e^{-\omega(x)}\big]
+B\big[\frac{1}{2\gamma}e^{\omega(x)}-\frac{\gamma}{2}e^{-\omega(x)}\big]=0,
\\
A\big[\frac{1}{2\gamma}e^{\omega(x)}-\frac{\gamma}{2}e^{-\omega(x)}\big]
+B\big[\frac{1}{2\gamma}e^{\omega(x)}+\frac{\gamma}{2}e^{-\omega(x)}\big]=0,
\end{gathered}\label{e18}
\end{equation}
where
$$
A:=z''(x)+z(x)\omega^{\prime^2}(x)-j\frac{z(x)}{\sqrt{z^2(x)-1}},\quad
B:=2z'(x)\omega'(x)+z(x)\omega^{\prime\prime}(x).
$$
Hence, $A=0$, $B=0$, i.e.
\begin{gather}
z''(x)+z(x)\omega^{\prime^2}(x)-j\frac{z(x)}{\sqrt{z^2(x)-1}}=0, \label{e19} \\
2z'(x)\omega'(x)+z(x)\omega^{\prime\prime}(x)=0. \label{e20}
\end{gather}
Integrating \eqref{e20}, we obtain
\begin{equation}
\omega'(x)=\frac{C_1}{z^2(x)}, \label{e21}
\end{equation}
where $C_1>0$ is an arbitrary constant. Substituting \eqref{e21} in \eqref{e19},
we obtain a single nonlinear second-order ODE for the function
$z(x)$
\begin{equation}
z''(x)+\frac{C_1^2}{z^3(x)}-j\frac{z(x)}{\sqrt{z^2(x)-1}}=0,\label{e22}
\end{equation}
with the initial-boundary conditions
\begin{gather}
\frac{1}{2\gamma}z(0)\big[e^{\omega(0)}+\gamma^2 e^{-\omega(0)}\big]=1,\quad
\frac{1}{2\gamma}z(0)\big[e^{\omega(0)}-\gamma^2 e^{-\omega(0)}\big]=0,
\label{e23} \\
\frac{1}{2\gamma}z'(0)\big[e^{\omega(0)}+\gamma^2 e^{-\omega(0)}\big]
+\frac{1}{2\gamma}z(0)\big[e^{\omega(0)}-\gamma^2 e^{-\omega(0)}\big]\omega'(0)=0,
\label{e24} \\
\frac{1}{2\gamma}z(1)\big[e^{\omega(1)}+\gamma^2 e^{-\omega(1)}\big]
=\varphi_1+1, \quad
\frac{1}{2\gamma}z(1)\big[e^{\omega(1)}-\gamma^2 e^{-\omega(1)}\big]=a_1. \label{e25}
\end{gather}
Conditions \eqref{e23}--\eqref{e25} have been borrowed from \eqref{e2} and
\eqref{e3} for the ansatz \eqref{e17}.
Now by \eqref{e21} we have
$\omega'(0)=\frac{C_1}{z^2(0)}$.
Multiplication of equation \eqref{e22} by the function $z'(x)\ne 0$ and
subsequent integration allow one to obtain
\begin{equation}
z^{\prime^2}(x)-\frac{C_1^2}{z^2(x)}-2j\sqrt{z^2(x)-1}=C_2. \label{e26}
\end{equation}
Here $C_2$ is an integration constant, which is chosen on account of the
initial boundary conditions \eqref{e23}, \eqref{e24}.
From conditions \eqref{e23}, one obtains
$z(0)=1$. Hence $\omega(0)=\ln\gamma$, $\gamma>0$.
Therefore, from \eqref{e26} we have
${z'}^2 (0)-C_1^2=C_2$, and from \eqref{e24} we have
$z'(0)=0$, hence $C_2=-C_1^2$. So, function $z(x)$, which satisfies the
conditions $z(0)=1$, $z'(0)=0$, may represent a solution of equation
$$
z'(x)=\pm\frac{\sqrt{C_1^2(1-z^2(x))+2jz^2(x)\sqrt{z^2(x)-1}}}{z(x)},
$$
from which it follows that
\begin{equation}
x=\pm\int_1^z \frac{sds}{\sqrt{C_1^2(1-s^2)+2js^2\sqrt{s^2-1}}}. \label{e27}
\end{equation}
So, we obtain $z$ as an implicit function of $x$.
Now, using \eqref{e21}, \eqref{e27} and
condition $\omega(x)=\ln\gamma$ at x = 0, it is possible to find $\omega$ as a
function of $z$:
$$
\omega(z)=\ln\gamma\pm\int_1^z
\frac{C_1ds}{s\sqrt{C_1^2(1-s^2)+2js^2\sqrt{s^2-1}}}.
$$
Now it is time to substitute conditions \eqref{e3} at the right end into
this equality and into \eqref{e27}. Taking into account the inequality
$\Theta(1)>0$ to
choose signs before integrals, we return to the system of equations \eqref{e14}, what
allows to define the values $C_1$ and $j$.
Therefore, the integration by the technique of replacement of variables lead
us to the solvability conditions for the boundary value problem similar
to those obtained with the application of the Liouville theorem.
It should be noted that integration \eqref{e1} by means of replacements
of the variables similar \eqref{e17}, was already carried out earlier by
a number of authors.
In \cite{dsz2013} replacement with a hyperbolic sine and a cosine is used
for separation of variables and the solution for \eqref{e1} singular
initial value problem on the left end of $[0, 1]$
\begin{equation}
\varphi(0)=0,\quad a(0)=0,\quad \varphi'(0)\equiv\frac{d\varphi}{dx}(0)=0,
\quad a'(0)\equiv\frac{d a}{dx}(0)=\beta\in\mathbb{R}. \label{e28}
\end{equation}
J. Batt stated the original approach to the solution of a singular boundary
value problem in several lectures at Irkutsk in August, 2009
(The authors do not know, whether these results were published or not).
In a preprint \cite{var2012} replacement \eqref{e17} is used with $\gamma=1$ and
also the above-stated singular initial value problem on the left end of a piece
$[0, 1]$ is considered.
This article differs from \cite{var2012} and others by two main stands:
We give strict definition of the solution of a singular BVP, and
we use an IVP only on the right end of the interval, at $x=1$
where conditions for existence and uniqueness are satisfied.
We will show the fundamental difference of our approach from that
of \cite{var2012} on a simple example.
\begin{example} \label{examp1} \rm
Let us consider a singular BVP
\begin{gather*}
\frac{d^2\varphi(x)}{dx^2}=\frac{2}{a(x)},\quad
\frac{d^2 a(x)}{dx^2}=-\frac{1}{\varphi(x)}, \\
\varphi(0)=0,\quad a(0)=0,\quad \varphi'(0)\equiv\frac{d\varphi}{dx}(0)=0, \\
\varphi(1)=\varphi_1=3,\quad a(1)=a_1=3/2.
\end{gather*}
The exact solution of this BVP in terms of Definition \ref{def1} is given by the
functions
\begin{equation}
\varphi(x)=3x^{4/3},\quad a(x)=\frac{3}{2}x^{2/3}. \label{e29}
\end{equation}
Here we have $a'(x)=x^{-1/3}$, therefore $\lim_{x\to +0}a'(x)=+\infty$.
So there is no number $\beta\in\mathbb{R}$ such that \eqref{e29} would be the
solution of corresponding~\cite{var2012} the IVP \eqref{e28}.
\end{example}
This example shows that transition from a BVP with singularity on the left
end of a piece to an IVP on the left end of a piece isn't correct and can
lead to loss of the solution of the BVP.
We will note also that singular BVP can have not less than two
solutions \cite{dsz2013} in certain cases.
\section{Constant solution $z(x)$}
In this section we consider the solutions of the initial system \eqref{e1},
which are constant in terms of variable $z$. Substituting
$z=z(x)\equiv {\rm const}$ into \eqref{e21} and \eqref{e22},
we obtain the following equations:
\begin{equation}
\omega'(x)=\frac{C_1}{z^2},\quad
\frac{C_1^2}{z}-j\frac{z}{\sqrt{z-1}}=0,\quad
z\equiv {\rm const}. \label{eD1}
\end{equation}
Let us show that the assumption of equality of $z(x)\equiv {\rm const}$
is not compatible with the boundary conditions at the interval's left end
for any value of parameter $\gamma$; there is no compatibility already with
the first two equalities of \eqref{e2}. Adding the first two
equalities from \eqref{e24},
one obtains $\gamma^{-1}z(0) e^{\omega(0)}=1$. When subtracting
one of these equalities from the other, one obtains $\gamma z(0) e^{-\omega(0)}=1$.
The result of multiplication of the two equalities obtained yields $z^2(0)=1$.
The latter is the initial condition, and function $z(x)$, which is constant,
shall satisfy this condition. Hence $z(x)\equiv 1$, what makes the denominator
in \eqref{e22} vanish. Therefore, constant solutions of the form
$z(x)\equiv {\rm const}$ are not compatible with the boundary conditions \eqref{e2}
at the interval's left end, when $x=0$.
Now introduce the condition $x_1=1$ and verify the condition of compatibility
of constant solutions $z(x)\equiv {\rm const}$ with the boundary conditions
\eqref{e3} at the right end of the interval, when $x=x_1$.
Adding equalities \eqref{e25}, one obtains
\begin{equation}
\gamma^{-1}z(x_1) e^{\omega(x_1)}=1+\varphi_1+a_1. \label{eD2}
\end{equation}
The result of subtraction of one of these equalities from the other yields
$\gamma z(x_1) e^{-\omega(x_1)}=1+\varphi_1-a_1$.
Similarly, the product of the two equalities obtained gives
$z^2(x_1)=(1+\varphi_1)^2-a_1^2=\Theta(x_1)+1$, and this is a positive number,
which shall satisfy the condition $z^2(x_1)\ne 1$.
Let $\Theta(x)=(1+\varphi(x))^2-a^2(x)-1$.
Since $(\varphi_1, a_1)\in \Omega_+$, it is obvious that $\Theta(x_1)>0$,
and condition $z^2(x_1)\ne 1$ holds.
It is also obvious from \eqref{e28} that
$C_1^2\Theta^{1/2}(x_1)-j(\Theta(x_1)+1)^2=0$.
Therefore, on account of parity and oddness of, respectively, \eqref{e21}
and \eqref{e22} with
respect to argument $z$, an arbitrary constant $C_1$ must be fixed, i.e.
$C_1=\sqrt{j}\Theta^{-1/4}(x_1)(\Theta(x_1)+1)$, and
$z(x_1)=\sqrt{\Theta(x_1)+1}$. Now, from \eqref{e21}, it is possible to obtain
$\omega'=\sqrt{j}\Theta^{-1/4}(x_1)$, whence
$\omega(x)=\sqrt{j}\Theta^{-1/4}(x_1)(x-x_1)+\omega_1$.
The value of the arbitrary constant $\omega_1$ is found from equality \eqref{e29}:
$\omega_1=\ln\frac{\gamma(1+\varphi_1+a_1)}{\sqrt{\Theta(x_1)+1}}$.
Now, substituting the obtained values of $z(x)$ and $\omega(x)$
in \eqref{e17},
we can formulate the solution of system \eqref{e1} in terms of the initial variables
\begin{equation}
\begin{aligned}
\varphi(x)=(1+\varphi_1)\cosh\lambda(x)+a_1\sinh\lambda(x)-1, \\
a(x)=(1+\varphi_1)\sinh\lambda(x)+a_1\cosh\lambda(x),
\end{aligned} \label{eD3}
\end{equation}
where $\lambda(x)=\sqrt{j}\Theta^{-1/4}(x_1)(x-x_1)$. Note that \eqref{eD1} gives an exact
solution of \eqref{e1}. It satisfies conditions \eqref{e3} at the right end of the interval,
where $x=x_1=1$, and does not satisfy conditions \eqref{e2} at the interval's
left end, where $x=0$. Furthermore, equality $x_1=1$ is not obligatory, so,
solution \eqref{eD1} may be used also in assigning the boundary conditions at any
point different from~$1$.
The parameter $j$ remains arbitrary. So, we have to choose this parameter
to satisfy not only \eqref{e3} but also
the following additional boundary conditions:
\begin{equation}
\varphi(x_0)=\varphi_0,\quad a(x_0)=a_0,\quad 00$, and conditions \eqref{e2} correspond to $\Theta(0)=0$,
variables \eqref{eD1} do not represent a solution of the boundary-value
problem \eqref{e1}--\eqref{e3}.
We intend to demonstrate that under some additional conditions imposed on
the solution obtained not at the endpoints of interval $[0, 1]$ but rather on
a countable set of its interior points one can obtain a solution of the
boundary-value problem (in the sense of Definition \ref{def1})
with the aid of \eqref{eD1}.
However, system \eqref{e1} describes the functions $\varphi(x)$ and $a(x)$
not on the whole interval $(0,1)$ but only on a sequence of subintervals.
Let there be given three numerical sequences $x_k$, $\alpha_{2k}$,
$\beta_{2k}$ $(k=1,2,\ldots)$ such that
\begin{gather*}
1=x_1>x_2>\dots >x_k>x_{k+1}>\dots >0,\quad \lim_{k\to +\infty}x_k=0, \\
0<\alpha_{2k}<\alpha_{*}<1,\quad 0<\beta_{2k}<\beta_{*}<1.
\end{gather*}
Consider system \eqref{e1} with condition \eqref{e3} on the interval
$[x_2, x_1]$ for the case when $x=x_1$
(we try to obtain a solution of this problem in the form \eqref{eD1}).
Compute the values $\varphi(x_2)$, $a(x_2)$ by formulas \eqref{eD1}. With
regard to the next interval $[x_3,x_2]$ let us assume that $\varphi(x)$, $a(x)$
satisfy the perturbed system \eqref{e1} with some unknown additive perturbations in
the right-hand side. This may be caused, for example, by the presence of some
contamination inside the diode or by an external disturbed field.
Let functions $\varphi(x)$, $a(x)$ remain unknown within the interval
$[x_3,x_2]$, but we know the relation between the values at the endpoints
\begin{equation}
\varphi(x_3)=\alpha_2\varphi(x_2),\quad a(x_3)=\beta_2 a(x_2). \label{eD7}
\end{equation}
In accordance with the previous computations of $\varphi(x_2)$ and $a(x_2)$,
it is possible to use \eqref{eD5} to find values of $\varphi(x_3)$, $a(x_3)$,
and to employ \eqref{eD1} in order to construct functions $\varphi(x)$,
$a(x)$ on the next interval $[x_4,x_3]$ (implying the replacement of
$x_1$ with $x_3$). Continuing this process sequentially,
it is possible to obtain exact solutions of \eqref{e1} in
the form \eqref{eD1} for the intervals $[x_{2k},x_{2k-1}]$
$(k=1,2,\ldots)$, while functions $\varphi(x)$, $a(x)$ remain unknown
on the intervals $(x_{2k-1},x_{2k-2})$ $(k=2, 3,\ldots)$.
Since functions $\cosh(\cdot)$, $\sinh(\cdot)$ are bounded on the restricted
interval $x\in [0, 1]$, conditions $0<\alpha_{*}<1$ and $0<\beta_{*}<1$
presume that functions $\varphi(x)$, $a(x)$ tend to zero
on intervals $[x_{2k}, x_{2k-1}]$ $(k=1,2,\ldots)$, when $k\to +\infty$.
This property is characterized by derivatives $\varphi'(x)$, $a'(x)$.
If additional perturbations in the right-hand sides of \eqref{e1}, which are effective
on the intervals $(x_{2k-1}, x_{2k-2})$ $(k=2, 3,\ldots)$, are bounded uniformly
in all cases $(k=2, 3\ldots)$, then functions $\varphi(x)$, $\varphi'(x)$
and $a(x)$ vanish, when $k\to +\infty$ also on the intervals
$(x_{2k-1},x_{2k-2})$. Consequently, functions $\varphi(x)$, $a(x)$ constructed
in the process of successive application of formulas \eqref{eD1} give a solution of
the boundary-value problem \eqref{e2} in the sense of Definition \ref{def1}. The diagrams
representing the functions $\varphi(x)$, $a(x)$ in the case, when
$j=\frac{1}{2}$, $\varphi_1=a_1=1$, $\alpha_{2k}=\beta_{2k}=\frac{1}{2}$,
$x_2=\frac{1}{2}$, $x_3=\frac{1}{4}$,
$x_4=\frac{1}{8}$, $x_5=\frac{1}{16},\ldots$, are shown below in
Figures \ref{fig1} and \ref{fig2}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig1} % {kss_fig1.eps}}
\end{center}
\caption{Function $\varphi(x)$. Solid line corresponds to formula \eqref{eD1}.
Dashed straight line connects the points where the exact behavior of the
function is unknown}
\label{fig1}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig2} %{kss_pic2.eps}}
\end{center}
\caption{Function $a(x)$. Solid line shows the solution obtained by formula
\eqref{eD1}. Dashed shows where the exact form of the function is unknown}
\label{fig2}
\end{figure}
\section{Asymptotic behavior}
Consider the function
$\Theta(x)=(1+\varphi(x))^2-a(x)^2-1$. Differentiating this
function several times and using \eqref{e1}, we obtain
the following identities
\begin{gather}
\Theta(x)\equiv (1+\varphi(x))^2-a(x)^2-1, \nonumber \\
\Theta'(x)\equiv 2(1+\varphi(x))\varphi'(x)-2a(x)a'(x), \nonumber\\
\Theta''(x)\equiv 2(\varphi'(x))^2-2(a'(x))^2+2j
\big(\Theta^{1/2}(x)+\Theta^{-1/2}(x)\big), \label{e37} \\
\Theta'''(x)\equiv j\big(3\Theta^{-1/2}(x)-\Theta^{-3/2}(x)\big)\Theta'(x).
\label{e38}
\end{gather}
Integrating the differential equation \eqref{e37}, we reduce its order by a unit
\begin{equation}
\Theta''(x)+c=2j(3\Theta^{1/2}(x)+\Theta^{-1/2}(x)). \label{e39}
\end{equation}
Here $c$ is an arbitrary constant. If we multiply equation \eqref{e39} by
$\Theta'(x)\ne 0$ and integrate the expression, we obtain
$$
\Theta'(x)=\pm\sqrt{8j}\sqrt{\Theta^{3/2}(x)
-\frac{c}{4j}\Theta(x)+\Theta^{1/2}(x)+\tilde{C}}.
$$
Here $\tilde{C}$ is an arbitrary constant. It is possible to transform this
equation by the replacement $H(x)=\Theta^{1/2}(x)$
$$
H(x)H'(x)=\pm\sqrt{2j}\sqrt{H^3(x)-\frac{c}{4j}H^2(x)+H(x)+\tilde{C}}.
$$
On account of the initial condition $\Theta(0)=0$ we have $H(0)=0$ and
$\tilde{C}=0$.
And finally we obtain
\begin{equation}
\int_{0}^{H}\frac{\sqrt{s}ds}{\sqrt{s^2-\frac{c}{4j}s+1}}=\pm\sqrt{2j}x. \label{e40}
\end{equation}
This formula determines $\Theta(x)$ as a function of argument $x$, where the
integral in the left hand side of \eqref{e38} is reduced to elliptic functions. And
in the case, when $c^2-16j^2=0$, the integral in \eqref{e38} is computed in terms of
elementary functions. From \eqref{e38} it is hardly ever possible to understand the
nature and the behavior of function $H(x)$ (i.e. $\Theta(x)$). So, it is necessary
to conduct asymptotic analysis. The arbitrary constant $c$ in \eqref{e37} may be
taken according to the boundary conditions \eqref{e2} and \eqref{e3}.
To this end, we take an arbitrary small positive number $\epsilon>0$ and
consider the initial value problem for \eqref{e1} under the following conditions:
$$
\varphi(\epsilon)=\epsilon,\quad \varphi'(\epsilon)=\epsilon,\quad
a(\epsilon)=\epsilon,\quad a'(\epsilon)=c_2
$$
where $c_2=p_2(0)$ is the value of the constant in the first integral $J_2$
in the solution of boundary-value problem \eqref{e1}--\eqref{e3}.
Let us denote the solution
of this initial-value problem as $(\bar{\varphi}(x), \bar{a}(x))$, and function
$\Theta$ taken along the solution is denoted as $\overline{\Theta}(x)$. Note,
for any sufficiently small $\epsilon>0$ the initial point lies in domain
$\Omega_{+}$, hence function $\overline{\Theta}(x)$ is defined at least on some
interval $[\epsilon, \epsilon+\delta)$ and satisfies identities \eqref{eD6},
\eqref{e37}.
Here $\delta>0$ is a positive number. From \eqref{eD6} we have
$$
\overline{\Theta}''(\epsilon)=
2\epsilon^2-2c_2^2+2j(\sqrt{2\epsilon}+\frac{1}{\sqrt{2\epsilon}}).
$$
Substituting this value into \eqref{e37} one can find
$c=2c_2^2-2\epsilon^2+4j\sqrt{2\epsilon}$.
It is important that this representation of constant $c$ contains only positive
powers of parameter $\epsilon$.
Now, passing to the limit with $\epsilon\to +0$, we can state that
function $\Theta(x)$ satisfies equation \eqref{e37} on the solution of the
boundary-value problem \eqref{e1}--\eqref{e3}, when $c=2c_2^2$ and $\Theta(0)=0$.
So, we shall look for an approximate solution of \eqref{e37} in the form
$\Theta(x)=kx^{\alpha}$, where the positive parameters, i.e. $\alpha$
and $k$, are to be defined. By equating (after the substation
into \eqref{e37}) the main terms in the left-hand and right-hand sides, and next by
separating the variables, we obtain $\alpha=\frac{4}{3}$ and
$k=\big(\frac{9j}{2}\big)^{2/3}$. Therefore, when values of
$x>0$ are small, function $\Theta(x)$ taken along the solution of the
boundary-value problem may be approximated by the formula
\begin{equation}
\tilde{\Theta}(x)=\big(\frac{9j}{2}\big)^{2/3}x^{4/3}. \label{e41}
\end{equation}
Substituting \eqref{e39} into \eqref{e1}, it is possible to decompose
this system into two independent scalar equations, for which the following approximate solutions
\begin{gather}
\tilde{\varphi}(x)=\frac{1}{2}\big(\frac{9j}{2}\big)^{2/3}x^{4/3},
\label{e42}\\
\tilde{a}(x)=c_2x\Big(1+\frac{1}{14}\big(\frac{9j}{2}\big)^{2/3}x^{4/3}\Big).
\label{e43}
\end{gather}
may be obtained similarly, on account of the conditions at the left end of the
interval, where $x=0$.
As obvious from \eqref{e40} and \eqref{e41}, the curve, which maps the solution
of the boundary-value problem onto the plane $(\varphi, a)$, may be approximated
in the vicinity of the origin by the function
\begin{equation}
a=c_22^{3/4}\big(\frac{9j}{2}\big)^{-1/2}\varphi^{3/4}
\big(1+\frac{1}{7}\varphi\big). \label{e44}
\end{equation}
Given by equation \eqref{e42}, the curve passes from inside of domain
$\Omega_{+}$ to its boundary at the origin. It has a vertical tangent at this
point, which touches also the boundary of domain $\Omega_{+}$.
It is incorrect to consider the initial-value problem \eqref{e1} with the
conditions \eqref{e2} and $a'(0)=c_2$ because the conditions of the theorem
on existence are violated under these initial conditions, and the substitution
of these initial conditions into the equation implies division by zero.
Meanwhile, using \eqref{e40}, \eqref{e41}, it is possible
to correctly state the initial-value problem for an arbitrary small positive
value of the independent variable $x=x_{0}=\epsilon>0$.
\section{Examples}
Consider problem \eqref{e1}--\eqref{e3} with the following conditions
for the right end of the interval: $q_1(1)=q_2(1)=Q_1=Q_2=1$. When solving the
corresponding system of equations \eqref{e14} with the aid of the iterative method,
one finds its approximate solution $c_2=u_{*}=0.8798287042$,
$j=v_{*}=0.5337203307$. This solution corresponds to the following values of
momentum $p_1(1)=-1.444231410$, $p_2(1)=1.162030057$ computed at
the right end of the interval by formulas \eqref{eD5} and \eqref{eD6}.
The Figures \ref{fig3} and \ref{fig4} represent diagrams demonstrating
components of the boundary-value problem solution (bold red line)
obtained by numerical right to left integration from
$t=1$ to $t=0$, and the solution asymptotics in the vicinity of
the interval's left end (thin blue line) computed by formulas
\eqref{e39}--\eqref{e42} at $t=1$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig3} %{kss_pic3.eps}
\end{center}
\caption{Solution of the boundary-value problem for the first
coordinate $q_1(t)$ (bold red line) and its asymptotic computed by formula
\eqref{e40} (blue).}
\label{fig3}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig4} %{kss_pic4.eps}}
\end{center}
\caption{Solution of the boundary value problem on the second
coordinate $q_2(t)$ (red) and it asymptotic by formula \eqref{e41}
(blue). The asymptotic behavior almost coincided with solution,
the curves merge}
\label{fig4}
\end{figure}
\begin{figure}[ht]
\includegraphics[width=0.6\textwidth]{fig5} % kss_pic8.eps
\caption{Solution of boundary value problem on the first
coordinate $q_1(t)$ (red) and its asymptotic by formula \eqref{e40}
(blue) for $q_1(1)=10, q_2(1)=1$}
\label{fig5}
\end{figure}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig6} % kss_pic9.eps}}
\end{center}
\caption{Solution of the boundary value problem on the second
coordinate $q_2(t)$ (red) and its asymptotic by formula \eqref{e41}
(blue) for $q_1(1)=10, q_2(1)=1$}
\label{fig6}
\end{figure}
Consider now the boundary-value problem \eqref{e1}--\eqref{e3} with the values,
which substantially differ at the right end of the interval, i.e.
$q_1(1)=Q_1=10$, $q_2(1)=Q_2=1$.
Solving system \eqref{e14} by the iteration method, one can
find its approximate solution $c_2=u^{*}=0.5404932672$, $j=v_{*}=12.14221503$.
It corresponds to the following values of momentum at the right end of the interval:
$p_1(1)=-16.33935466$, $p_2(1)=1.534531629$. Consider only
the following two diagrams of coordinate variations for the given variant of
initial data.
In this variant, the curves do not merge completely as we see in
figure \ref{fig4}. Likewise for all other diagrams, we obviously see
high precision of the solution for the boundary-value problem in the vicinity
of the left end of the interval, represented by formulas \eqref{e39}--\eqref{e42}.
\subsection*{Acknowledgments}
This research was supported by the Russian Foundation for Basic Research
(project no.~15-08-06680).
\begin{thebibliography}{99}
\bibitem {bat93} J. Batt, K. Fabian;
\emph{Stationary solutions of the relativistic
Vlasov-Maxwell system of plasma physics}, Chinese Annals of Mathematics,
Ser.~B, 14(3) (1993), pp.~253--1278.
\bibitem {ben98} N. Ben Abdallah, P. Degond, F. Mehats;
\emph{Mathematical models of magnetic insulation}, Physics of Plasma,
5 (1998), pp.~1522-1534.
\bibitem {cd2000} L. Caffarelli, J. Dolbeault, P. A. Markowich, C. Schmeiser;
\emph{On Maxwellian equilibria of insulated semiconductors}, Interfaces
Free Bound, 2 (2000), pp.~331--339.
\bibitem {b2007} P. Crispel, P. Degond, M. Vignal;
\emph{A plasma expansion model based on the full Euler-Poisson system},
Mathematical Models and Methods in Applied Sciences, 17 (2007), pp.~1129--1158.
\bibitem {dsz2013} S. Djebali, O. Saifi, S. Zahar;
\emph{Singular boundary-value problems with variable
coefficients on the positive half-line},
Electronic Journal of Differential Equations, Vol. (2013), N~73, pp.~1--18.
\bibitem{lxz2014} Yuanjie Lei, Linjie Xiong, Huijiang Zhao;
\emph{One-species Vlasov-Poisson-Landau system near Maxwellians in the whole space},
Kinetic and Related Models (KRM). 2014. V.~7. N~3. pp.~551--590.
\bibitem {ss2010} E. I. Semenov, A. V. Sinitsyn;
\emph{Mathematical model of magnetic insulation vacuum diode and its exact solutions},
News of Irkutsk State University. Series Mathematics, 3 (2010), N~1, pp.~78--91.
\bibitem {sin2001} A. V. Sinitsyn;
\emph{Positive solutions of nonlinear singular
boundary-value problem of magnetic insulation}, Mathematical
Models and Computer Simulation. 13 (2001), N~5, pp.~37--52.
\bibitem {var2012} V. P. Varin;
\emph{An analysis of a vacuum diode model}.
Keldysh Institute preprints, 8 (2012), 19 p.
\bibitem {w1927} E.T. Whittaker;
\emph{A treatise on the analytical dynamics of Particles and
Rigid Bodies: With an Introduction to the Problem of Three Bodies},
Cambridge University Press, 1927.
\end{thebibliography}
\end{document}