Transmutation operators and a new representation for solutions of perturbed Bessel equations

In the present work we consider a second order singular differential equation

where $\ell$ is a real number, $\ell\geq-\frac{1}{2}$, $q$ is a complex-valued function satisfying the condition

and $\omega$ is a (complex) spectral parameter. Equations of the form (1.1) appear naturally in many real-world applications after a separation of variables and therefore have received considerable attention (see, e.g., [4], [7], [8], [9], [14], [18], [29, Sect. 3.7], [32], [33], [39]).

In [21] a new representation for solutions of one-dimensional Schrödinger equation (the case $\ell=0$ in (1.1)) in the form of Neumann series of Bessel functions (NSFB) [3, 40] was obtained. The representation possesses such remarkable properties as exponentially fast convergence for smooth potentials and uniform error bound for partial sums for all $\omega\in\mathbb{R}$. The idea behind this representation is the existence of the transmutation operator connecting a solution of the equation with the solution of the simpler equation having $q\equiv 0$. A transmutation operator can be realized in the form of a Volterra integral operator. Expanding its integral kernel into a Fourier-Legendre series allowed us to obtain the NSBF representation for the solution.

In [24] we proposed an NSBF representation for the regular solution $\tilde{y}(\omega,x)$ of (1.1) normalized by its asymptotics $\tilde{y}(\omega,x)\sim x^{\ell+1}$ at zero. However, the above mentioned advantages of the regular case were lost for the representation from [24]. The reasons are the following. First, for any fixed $x>0$ the function $\tilde{y}(\omega,x)$ decays as $\frac{1}{\omega^{\ell+1}}$ as $\omega\to\infty$. For that reason any error bound uniform with respect to $\omega\in\mathbb{R}$ is useful only in some neighborhood of $\omega=0$, and for large values of the parameter $\ell$ this neighborhood becomes rather small. Second, the exponential convergence of the partial sums was lost for $\ell\not\in\mathbb{N}$, the convergence was bounded by $\frac{1}{N^{2\ell+3}}$ (here $N$ is the truncation parameter) independently of the smoothness of the potential. The reason was in the underlying Mehler’s integral representation for the solution $\tilde{y}$:

here $R=2\tilde{R}$ and is an even function, and $j_{\ell}$ denotes the spherical Bessel function of the first kind, see [1, Section 10.1]. The integral kernel $\tilde{R}$ admits the following representation

where $c(x)=2^{\ell-3/2}\Gamma\left(\ell+\frac{3}{2}\right)\int_{0}^{x}q(s)\,ds$ and $\tilde{R}_{1}$, extended by 0 onto the whole line, possesses at least $[\ell+2]$ derivatives as a function of $t$. As a result, for non-integer values of $\ell$ the decay rate of the Fourier-Legendre coefficients of $R$ is polynomial and is determined by the first term in (1.2), independent on the smoothness of $q$.

To overcome the first problem, one can use the transmutation operator for the perturbed Bessel equation. Recall that a transmutation operator intertwining (1.1) with the unperturbed Bessel equation

has the form [38, 34, 10, 37, 16]

In such case, if one takes $y(\omega,x)=\omega xj_{\ell}(\omega x)$ as a regular solution of (1.3) (note that this solution remains bounded and does not decay as $\omega\to\infty$, see [1, (9.2.1)]), finds an approximation $K_{N}$ for $K$ satisfying $\|K(x,\cdot)-K_{N}(x,\cdot)\|_{L_{2}(0,x)}\leq\varepsilon_{N}(x)$ and defines an approximate solution $u_{N}$ by using this $K_{N}$ in (1.4), one easily obtains (applying the Cauchy-Schwartz inequality) that

uniformly for $\omega\in\mathbb{R}$. And since $u(\omega,x)$ does not decay as $\omega\to\infty$, such approximation is useful for both small and large values of $\omega$.

However, the problem of constructing an approximation $K_{N}$ for which at the same time the integrals in (1.4) can be easily evaluated and rapid decay of the error bound $\varepsilon_{N}$ as $N\to\infty$ can be proved, is not an easy task, see, e.g., [25] for an attempt. The integral kernel $K$ is a solution of a singular Goursat problem which can be transformed into one of the several equivalent integral equations, see [38], [10], [16]. We do not expect that any of these integral equations can provide an efficient approximation method for the integral kernel. The main reason is the following: we conjecture that the integral kernel has the form

where $\widetilde{K}$ is a nice function (in the sense that it is at least $C^{\infty}$ in the second variable for sufficiently smooth potentials). The base for such conjecture is the formula (4.6) from [9]

proved there under the assumption that $q$ possesses a holomorphic extension onto the disk of radius $2xe\sqrt{1+2|\ell|}$. We believe that the factor $\frac{t^{\ell+1}}{x^{\ell}}$ plays a crucial role in the construction of a good approximation to $K$, however it gets strongly obscured in the known integral equations for $K$.

In [23] we proposed to use a series representation for the integral kernel $R$ and an Erdelyi-Kober operator to obtain a series representation for the integral kernel $K$, and for integer values of $\ell$ the partial sums provided a good approximation for $K$. For non-integer values of $\ell$, the approximation obtained was not practical to substitute it into (1.4).

In the present paper we show that the Legendre polynomials are not the best choice for representing the kernel $R$, and one needs to use the Jacobi polynomials instead. The motivation for this claim is the formula (1.2) and the following representation for $R$

which follows from formula (4.4) from [9] and is proved there under the assumption that $q$ possesses a holomorphic extension onto the disk of radius $2xe\sqrt{1+2|\ell|}$.

We show that the factor $\bigl{(}1-(t/x)^{2}\bigr{)}^{\ell+1}$ should be used as a weight for the Fourier-Jacobi expansion, and $\tilde{R}$ can be represented as

where $P_{n}^{(\alpha,\beta)}$ stands for the Jacobi polynomials, leading to the following representation for the integral kernel $K$

where we denoted

As a result, the following representation for the regular solution of (1.1) is obtained

We show the uniform with respect to $\omega\in\mathbb{R}$ convergence rate estimates for (1.8), prove faster than polynomial convergence for $C^{\infty}$ potentials and present efficient for numerical implementation recurrent formulas to calculate the coefficients $\beta_{n}$ in which only integration is used, and no differentiation is required.

For the derivative of the regular solution $u(\omega,x)$ we obtain the following representation

possessing all the remarkable properties of (1.8): uniform with respect to $\omega\in\mathbb{R}$ error bounds for truncated series, faster than polynomial convergence for $C^{\infty}$ potentials and efficient for numerical implementation formulas to calculate the coefficients $\gamma_{n}$.

The representation (1.6) for the integral kernel $K$ together with the decay rate estimates for the coefficients $\beta_{n}$ can be of the independent interest for studying the transmutation operator for perturbed Bessel equation and for solving the inverse spectral problem like it was done in the regular case in [19], [11], [20].

We would like to mention that despite of lots of technical details in proofs, application of the representations (1.8) and (1.9) for numerical solution of boundary value or spectral problems for equation (1.1) is very simple. All that one needs is a solution $u_{0}$ of the equation

and its derivative $u_{0}^{\prime}$ (both can be obtained numerically). Then one calculates the coefficients $\{\beta_{n}\}_{n=0}^{N_{1}}$ and $\{\gamma_{n}\}_{n=0}^{N_{2}}$ using (7.7)–(7.9) and (7.11). To estimate optimal values of $N_{1}$ and $N_{2}$ one utilizes the equalities (8.1) and (8.3). Computation of both $u(\omega,x)$ and $u^{\prime}(\omega,x)$ for each fixed $\omega$ reduces to the computation of the values of several Bessel functions. As a result, for example, several hundreds eigenvalues of a Sturm-Liouville problem can be obtained in less than a second.

The paper is organized as follows. In Section 2 we present some preliminary information about transmutation operators and asymptotic estimates of the solutions. In Section 3 we study the smoothness of the integral kernel $\tilde{R}$ (Propositions 3.3 and 3.4) and convergence rate of its Fourier-Jacobi expansion (Theorem 3.5 and Proposition 3.6). In Section 4 we show that the integral kernel $K$ possesses Fourier-Jacobi representation (1.6) (Theorem 4.1), extend this representation to a wider class of potentials and prove its convergence rate (Theorem 4.3). In Section 5 we obtain the representation (1.8) (Theorem 5.2). In Section 6 we present necessary facts about the derivative of the regular solution, its relation with the transmutation operator and prove the representation (1.9) (Theorems 6.2 and 6.3). In Section 7 we obtain the recurrent formulas for the coefficients $\beta_{n}$ and $\gamma_{n}$. In Section 8 we present numerical results for solution of two spectral problems and study the decay of the coefficients $\beta_{n}$. In Appendix A we prove that the transmutation integral kernel depends continuously on the potential (Corollary A.3). Finally, in Appendix B we prove that the regular solution of equation (1.1) is non-vanishing on the whole $(0,b]$ for any sufficiently large negative $\lambda=\omega^{2}$ (Proposition B.1).

Denote

This function is a regular solution of equation (1.3) satisfying the asymptotic condition

By $u(\omega,x)$ we denote the regular solution of (1.1) satisfying the same asymptotic condition at 0.

Assume that $q\in C[0,b]$. Then [38], [34], [10], [16] there exists a unique continuous kernel $K(x,t)$ such that for all $\omega\in\mathbb{C}$

The integral kernel $K$ satisfies

The operator $T$ is called the transmutation operator.

The existence of the transmutation operator (2.3) can be established for wider class of potentials. For the purposes of this work we need that the integral kernel $K$ as a function of $t$ belongs to $L_{2}(0,x)$ for each fixed $x$. Let us introduce the notation [17]

Then the following condition on $q$ is sufficient for the integral kernel $K$ to belong to $L_{2}(0,x)$:

For the proof we refer the reader to [34, §2], the case of non-integer values of $\ell$ is also covered taking [13] into account.

The difference between the solutions $u(\omega,x)$ and $b_{\ell}(\omega x)$ satisfies the inequality [17, (2.18)]

from which it follows for the potential $q$ satisfying condition (2.6) that

More can be said about the asymptotic behavior of the solution as $\omega\to\infty$ if the potential $q$ possesses several derivatives on the whole segment $[0,b]$. Indeed, similarly to the proof of Proposition 4.5 from [24] (see also [15, Theorem 4]) one can see that if $q\in W_{1}^{2p-1}[0,b]$, $p\in\mathbb{N}$, i.e., $q$ possesses $2p-1$ derivatives, the last one belonging to $L_{1}(0,b)$, then

where

and $\mathcal{R}_{p}$ as a function of $\omega$ is an even entire function of exponential type $x$.

Application of the Paley-Wiener theorem leads to the following Mehler-type integral representation for the solution $u(\omega,x)$, see [24]

where

By $W_{2}^{\alpha}(\mathbb{R})$, $\alpha\geq 0$ we denote the fractional-order Sobolev space, also called Bessel potential space [2, Chap. 7] consisting of functions satisfying $f\in L_{2}(\mathbb{R})$ and $(1+|\xi|^{2})^{\alpha/2}\mathcal{F}[f](\xi)\in L_{2}(\mathbb{R})$, where $\mathcal{F}$ is the Fourier transform operator. Then, extending the integral kernel $\tilde{R}$ as a function of $t$ by 0 outside of $[-x,x]$ and denoting the resulting function again by $\tilde{R}$ we have the following [24, Proposition 4.1]. $\tilde{R}$ is a continuous, even, compactly supported on $[-x,x]$ function and such that $\tilde{R}\in W_{2}^{\ell+3/2-\mu-\varepsilon}(\mathbb{R})$ for any sufficiently small $\varepsilon>0$.

Additionally, if $q\in W_{1}^{2p-1}[0,b]$, then [24, (4.24)]

where $\tilde{R}_{p}$ (extended by 0 outside of $[-x,x]$) is a continuous, even function satisfying $\tilde{R}_{p}\in W_{2}^{\ell+p+3/2-\varepsilon}(\mathbb{R})$ for any sufficiently small $\varepsilon>0$.

The integral kernels $\tilde{R}$ and $K$ are related by the following relations [23, (2.5) and (2.6)]

and

here $n$ is an arbitrary integer satisfying $n>\ell+1$.

In this section we show that the integral kernel $\tilde{R}$ can be expanded into a Fourier-Jacobi series (1.5), study error bounds for the remainder of the truncated series and the decay rate of the coefficients $\beta_{n}$.

We recall some definitions from the approximation theory used in this section. By $[[x]]$ we denote the largest integer smaller than $x$, and let $\{\{x\}\}:=x-[[x]]$. Then $\{\{x\}\}\in(0,1]$.

Following [31] we say that a function $f$ is Lipschitz of an order $\alpha>0$ on $I$ (may be a segment or the whole line) if

1. 1.

   there exist the derivatives of $f$ of all orders up to the order $[[\alpha]]$;

2. 2.

   $f^{(m)}\in L_{\infty}(I)$ for all $m\leq[[\alpha]]$;

3. 3.

   $f^{([[\alpha]])}$ satisfies the Lipschitz condition of order $\{\{\alpha\}\}$, i.e., there exists a positive constant $A$ such that

   $$|f^{([[\alpha]])}(x)-f^{([[\alpha]])}(y)|\leq A|x-y|^{\{\{\alpha\}\}}\qquad\forall x,y\in I.$$

We will denote the class of such functions by $\operatorname{Lip}(\alpha,I)$.

Following [28] we say that a function $f$ belongs to the Zygmund class $\operatorname{Zyg}(\alpha,I)$ on some $I$ for some $\alpha>0$ if

1. 1.

   $f^{(m)}\in C(I)$ for all $m\leq[[\alpha]]$;

2. 2.

   $f^{([[\alpha]])}$ satisfies the Zygmund condition for some constant $C$, i.e.,

   $$|f^{([[\alpha]])}(x+h)-2f^{([[\alpha]])}(x)+f^{([[\alpha]])}(x-h)|\leq Ch^{\{\{\alpha\}\}}\qquad\forall x,h:\ \{x-h,x,x+h\}\subset I.$$

Note that $\operatorname{Lip}(\alpha,I)\subset\operatorname{Zyg}(\alpha,I)$.

The following proposition immediately follows from [28, Theorems 1 and 3] and shows the relationship between the decay rate of the Fourier transform and the smoothness of a function.

Following [26] we introduce the following notations and definitions. Let $\mathcal{P}_{n}$ be the set of algebraic polynomials of degree not greater than $n$. Let $W_{\alpha}(x)=(1-x^{2})^{\alpha/2}$, $x\in[-1,1]$, $\alpha>-1/2$. We define as best weighted polynomial approximation of a function $f$ such that $fW_{\alpha}\in L_{2}(-1,1)$ the quantity

Denote by $S(\alpha;f,x)$ the orthonormal expansion of $f$ with respect to the system of normalized Jacobi polynomials $\{\tilde{P}_{n}^{(\alpha,\alpha)}(x)\}_{n=0}^{\infty}$, that is

where

Then the best weighted polynomial approximation of the function $f$ coincides with the norm of the remainder of its Fourier-Jacobi expansion:

Let $S_{k}^{(\alpha)}$, $k=1,2,\ldots$ be the set of functions $f$ satisfying

1. 1.

   $f$ is a $k$-times iterated integral of $f^{(k)}$ in $(-1,1)$;

2. 2.

   $f^{(l)}W_{\alpha+l}\in L_{2}(-1,1)$, $l=0,1,\ldots,k$.

By $S_{0}^{(\alpha)}$ we define the set of functions $f$ satisfying $W_{\alpha}f\in L_{2}(-1,1)$.

Consider the following modulus of continuity:

where $f^{\ast}(\theta)$ is defined by $f^{\ast}(\theta):=f(\cos\theta)$, $0\leq\theta\leq\pi$.

Then the following result holds

In the rest of this section we apply Proposition 3.1 and Theorem 3.2 to estimate the decay rate of the coefficients and of the remainder of the Fourier-Jacobi expansion of the integral kernel $\tilde{R}$.

Let $x>0$ be fixed. Consider the functions

The expansion (1.5) reduces to the Fourier-Jacobi expansion of the function $h$.