Spectral parameter power series representation for solutions of linear system of two first order differential equations

TL;DR

This paper introduces a spectral parameter power series (SPPS) method for solving linear systems of two first-order differential equations, providing a recursive integral-based series representation for solutions with broad applicability.

Contribution

The paper develops a new SPPS representation for solutions of first-order linear systems with arbitrary matrix coefficients, including a transformation for more general systems and numerical implementation guidance.

Findings

SPPS series converges and exists for the systems considered.

Transformation reduces general systems to the form suitable for SPPS.

Numerical examples demonstrate the effectiveness of the method.

Abstract

A representation in the form of spectral parameter power series (SPPS) is given for a general solution of a one dimension Dirac system containing arbitrary matrix coefficient at the spectral parameter, \[ B \frac{dY}{dx} + P(x)Y = \lambda R(x)Y,\] where $Y=(y_1,y_2)^T$ is the unknown vector-function, $\lambda$ is the spectral parameter, $B = \begin{pmatrix}0 & 1 \\ -1 & 0\end{pmatrix}$, and $P$ is a symmetric $2\times 2$ matrix, $R$ is an arbitrary $2\times 2$ matrix whose entries are integrable complex-valued functions. The coefficient functions in these series are obtained by recursively iterating a simple integration process, beginning with a non-vanishing solution for one particular $\lambda = \lambda_0$. The existence of such solution is shown. For a general linear system of two first order differential equations \[ P(x)\frac{dY}{dx}+Q(x)Y = \lambda R(x)Y,\ x\in [a,b], \] where $P$, $Q$, $R$ are $2\times 2$ matrices whose entries are integrable complex-valued functions, $P$ being invertible for every $x$, a transformation reducing it to a type considered above is shown. The general scheme of application of the SPPS representation to the solution of initial value and spectral problems as well as numerical illustrations are provided.

Tables

Table 1. Table 1: Eigenvalues of the spectral problem ( 4.5 ), ( 4.6 ) and absolute errors of the approximate eigenvalues obtained using the representation ( 2.16 ) truncated to N = 100 𝑁 100 N=100 terms, using additionally the spectral shift technique and of those reported in [ 2 ] .

$n$	${\lambda }_n$	Abs. error,	Abs. error,	Abs. error,
		directly from (2.16)	using spectral shift	reported in [2]
-100	-313.9101939150852		$5.3\cdot {10}^{-6}$
-50	-156.8314900718780		$9.6\cdot {10}^{-8}$
-20	-62.58649828127890		$3.8\cdot {10}^{-10}$
-10	-31.17522014114365	$1.8\cdot {10}^{-3}$	$5.1\cdot {10}^{-12}$
-7	-21.75442521496494	$1.7\cdot {10}^{-7}$	$5.1\cdot {10}^{-13}$
-5	-15.47654249528427	$6.4\cdot {10}^{-10}$	$6.0\cdot {10}^{-14}$
-2	-6.079080595285440	$3.5\cdot {10}^{-13}$	$1.4\cdot {10}^{-14}$	$3.5\cdot {10}^{-13}$
-1	-2.977189710951455	$5.1\cdot {10}^{-14}$	$5.1\cdot {10}^{-14}$	$2.6\cdot {10}^{-13}$
0	1	$1.9\cdot {10}^{-14}$	$1.9\cdot {10}^{-14}$	$3.1\cdot {10}^{-13}$
1	3.478833069692201	$1.4\cdot {10}^{-14}$	$1.4\cdot {10}^{-14}$	$1.2\cdot {10}^{-12}$
2	6.578592238156064	$3.1\cdot {10}^{-13}$	$1.2\cdot {10}^{-14}$
5	15.97642352997195	$1.1\cdot {10}^{-9}$	$5.1\cdot {10}^{-14}$
7	22.25436259528469	$4.2\cdot {10}^{-7}$	$5.1\cdot {10}^{-12}$
10	31.67518895778715	$1.5\cdot {10}^{-3}$	$5.0\cdot {10}^{-12}$
20	63.08649039551696		$3.7\cdot {10}^{-10}$
50	157.3314888061299		$9.5\cdot {10}^{-8}$
100	314.4101935985044		$5.2\cdot {10}^{-6}$

Table 2. Table 2: Eigenvalues of the spectral problem ( 4.9 ) and absolute errors of the approximate eigenvalues obtained converting the problem to a Dirac system and directly using the SPPS representation from [ 17 ] with two different spectral shift strategies.

		Abs. error,	Abs. error,	Abs. error,
$n$	${\lambda }_n$	using the Dirac	using spectral shift to	using a constant step
		system approach	the nearest eigenvalue	size spectral shift
0	1.519865821099347	$6.4\cdot {10}^{-13}$	$2.0\cdot {10}^{-14}$	$4.2\cdot {10}^{-14}$
1	4.943309822144690	$4.0\cdot {10}^{-13}$	$6.8\cdot {10}^{-14}$	$6.3\cdot {10}^{-13}$
2	10.28466264508758	$8.7\cdot {10}^{-13}$	$4.4\cdot {10}^{-13}$	$2.6\cdot {10}^{-12}$
3	17.55995774641423	$1.3\cdot {10}^{-12}$	$3.8\cdot {10}^{-13}$	$3.5\cdot {10}^{-12}$
5	37.96442586193434	$4.2\cdot {10}^{-13}$	$3.1\cdot {10}^{-13}$	$4.7\cdot {10}^{-12}$
10	123.4977068009282	$3.9\cdot {10}^{-12}$	$1.4\cdot {10}^{-12}$	$4.5\cdot {10}^{-12}$
25	678.9217784771679	$1.0\cdot {10}^{-11}$	$1.0\cdot {10}^{-12}$	$6.0\cdot {10}^{-12}$
50	2604.036332024594	$1.4\cdot {10}^{-11}$	$1.4\cdot {10}^{-11}$	$3.9\cdot {10}^{-11}$
75	5779.062267233881	$2.6\cdot {10}^{-11}$	$1.1\cdot {10}^{-11}$	$5.9\cdot {10}^{-8}$
100	10204.07191390758	$7.1\cdot {10}^{-11}$	$1.1\cdot {10}^{-7}$	$3.8\cdot {10}^{-5}$
150	22804.07903279700	$1.3\cdot {10}^{-10}$	$1.1\cdot {10}^{-2}$	$2.4\cdot {10}^{-3}$

Equations

$$B\frac{dY}{dx}+P(x)Y=\lambda R(x)Y,$$

$$P(x)\frac{dY}{dx}+Q(x)Y=\lambda R(x)Y,x\in [a,b],$$

$$\begin{cases}v^{\prime}+p_{1}(x)u+q(x)v=\lambda\bigl{(}r_{11}(x)u+r_{12}(x)v\bigr{)},\\ -u^{\prime}+q(x)u+p_{2}(x)v=\lambda\bigl{(}r_{21}(x)u+r_{22}(x)v\bigr{)},\end{cases}$$

$$B\frac{dY}{dx}+P(x)Y=\lambda R(x)Y,Y(x)=\left(\begin{array}{c}u(x)\\ v(x)\end{array}\right),$$

$$B=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right),P(x)=\left(\begin{array}{cc}{p}_1(x)& q(x)\\ q(x)& p_2(x)\end{array}\right),R(x)=\left(\begin{array}{cc}{r}_{11}(x)& r_{12}(x)\\ r_{21}(x)& r_{22}(x)\end{array}\right),$$

$$B\frac{dY}{dx}+P(x)Y=0$$

$$X^{(0)}(x)$$

$$Y^{(0)}(x)$$

$$Z^{(n)}(x)$$

$$X^{(n+1)}(x)$$

$$Y^{(n+1)}(x)$$

$${\tilde{X}}^{(0)}(x)$$

$${\tilde{Y}}^{(0)}(x)$$

$$\bigl{(}p(x)u^{\prime}\bigr{)}^{\prime}+q(x)u=\omega^{2}r(x)u,$$

$${\mathcal{X}}^{(0)}(x)$$

$${\mathcal{X}}^{(n)}(x)$$

$${\tilde{\mathcal{X}}}^{(n)}(x)$$

$$u=c_1{u}_0\sum _{k=0}^{\infty }\frac{{\omega }^{2k}{\tilde{\mathcal{X}}}^{(2k)}}{(2k)!}+c_2{u}_0\sum _{k=0}^{\infty }\frac{{\omega }^{2k}{\mathcal{X}}^{(2k+1)}}{(2k+1)!}.$$

$$\omega \frac{v}{p}=u_0\frac{d}{dx}\left(\frac{u}{u_0}\right)=u^{\prime }-\frac{u_0^{\prime }}{u_0}u.$$

$$\{\begin{array}{l}{v}^{\prime }+\frac{u_0^{\prime }}{u_0}v=\omega ru,\\ -u^{\prime }+\frac{u_0^{\prime }}{u_0}u=-\omega \frac{1}{p}v,\end{array}$$

$$X^{(2n)}=Y^{(2n+1)}\equiv 0,X^{(2n-1)}={\mathcal{X}}^{(2n-1)},Y^{(2n)}={\mathcal{X}}^{(2n)},$$

$${\tilde{X}}^{(2n+1)}={\tilde{Y}}^{(2n)}\equiv 0,{\tilde{X}}^{(2n)}={\tilde{\mathcal{X}}}^{(2n)},{\tilde{Y}}^{(2n+1)}={\tilde{\mathcal{X}}}^{(2n+1)}.$$

$$Y=c_1{Y}_1+c_2{Y}_2=c_1\left(\begin{array}{c}{u}_1\\ v_1\end{array}\right)+c_2\left(\begin{array}{c}{u}_2\\ v_2\end{array}\right),$$

$$\left(\begin{array}{c}{u}_1\\ v_1\end{array}\right)=\sum _{n=0}^{\infty }\frac{{\lambda }^n}{n!}\left(\begin{array}{c}f{\tilde{X}}^{(n)}\\ g{\tilde{Y}}^{(n)}\end{array}\right),\left(\begin{array}{c}{u}_2\\ v_2\end{array}\right)=\sum _{n=0}^{\infty }\frac{{\lambda }^n}{n!}\left(\begin{array}{c}f{X}^{(n)}\\ g{Y}^{(n)}\end{array}\right).$$

$$Y_1(x_0)=\left(\begin{array}{c}f(x_0)\\ 0\end{array}\right),Y_2(x_0)=\left(\begin{array}{c}0\\ g(x_0)\end{array}\right).$$

$$\{\begin{array}{l}{v}^{\prime }+p_{1}u+qv=h_1,\\ -u^{\prime }+qu+p_{2}v=h_2,\end{array}$$

$$u(x_0)=v(x_0)=0,$$

$$u(x)=f(x)\int_{x_{0}}^{x}\biggl{(}-\frac{h_{2}(t)}{f(t)}+\frac{p_{2}(t)}{f^{2}(t)}\int_{x_{0}}^{t}\bigl{(}f(s)h_{1}(s)+g(s)h_{2}(s)\bigr{)}\,ds\biggr{)}dt$$

$$v(x)=g(x)\int_{x_{0}}^{x}\biggl{(}\frac{h_{1}(t)}{g(t)}+\frac{p_{1}(t)}{g^{2}(t)}\int_{x_{0}}^{t}\bigl{(}f(s)h_{1}(s)+g(s)h_{2}(s)\bigr{)}\,ds\biggr{)}dt.$$

$${\left(\frac{u^{\prime }}{p_2}\right)}^{\prime }+\left(p_1-{\left(\frac{q}{p_2}\right)}^{\prime }-\frac{q^2}{p_2}\right)u$$

$${\left(\frac{v^{\prime }}{p_1}\right)}^{\prime }+\left(p_2+{\left(\frac{q}{p_1}\right)}^{\prime }-\frac{q^2}{p_1}\right)v$$

$$\begin{gathered} u(x)=f(x){\int }_{x_0}^x\frac{p_2(t)}{f^2(t)}{\int }_{x_0}^{t}f(s)\left(h_1(s)-\frac{q(s)}{p_2(s)}{h}_2(s)-{\left(\frac{h_2(s)}{p_2(s)}\right)}^{\prime }\right)dsdt \\ +c_{11}f(x)+c_{12}f(x){\int }_{x_0}^x\frac{p_2(s)}{f^2(s)}ds \end{gathered}$$

$$\begin{gathered} v(x)=g(x){\int }_{x_0}^x\frac{p_1(t)}{g^2(t)}{\int }_{x_0}^{t}g(s)\left(h_2(s)-\frac{q(s)}{p_1(s)}{h}_1(s)+{\left(\frac{h_1(s)}{p_1(s)}\right)}^{\prime }\right)dsdt \\ +c_{21}g(x)+c_{22}g(x){\int }_{x_0}^x\frac{p_1(s)}{g^2(s)}ds. \end{gathered}$$

$$u(x_0)=v(x_0)=0,u^{\prime }(x_0)=-h_2(x_0),v^{\prime }(x_0)=h_1(x_0).$$

$$c_{12}=-\frac{h_2(x_0)f(x_0)}{p_2(x_0)}\text{and}{c}_{22}=\frac{h_1(x_0)g(x_0)}{p_1(x_0)}.$$

Full text

[TABLE]

where $Y=(y_{1},y_{2})^{T}$ is the unknown vector-function, $\lambda$ is the spectral parameter, $B=\begin{pmatrix}0&1\\ -1&0\end{pmatrix}$, and $P$ is a symmetric $2\times 2$ matrix, $R$ is an arbitrary $2\times 2$ matrix whose entries are integrable complex-valued functions. The coefficient functions in these series are obtained by recursively iterating a simple integration process, beginning with a non-vanishing solution for one particular $\lambda=\lambda_{0}$. The existence of such solution is shown.

For a general linear system of two first order differential equations

[TABLE]

where $P$, $Q$, $R$ are $2\times 2$ matrices whose entries are integrable complex-valued functions, $P$ being invertible for every $x$, a transformation reducing it to a system ($\ast$) is shown.

The general scheme of application of the SPPS representation to the solution of initial value and spectral problems as well as numerical illustrations are provided.

1 Introduction

The spectral parameter power series (SPPS) representation for solutions of second-order linear differential equations [14], [17] has proven to be an efficient tool for solving (analytically and numerically) and studying a variety of problems, see the review [13] and recent papers [4], [5], [6], [11], [23], [26]. The SPPS method starts with a non-vanishing solution of the equation for one fixed value of the spectral parameter and by performing a series of recursive integrations produces coefficients of the Taylor series of the solution with respect to the spectral parameter. The procedure can be easily and efficiently implemented numerically allowing one to solve a variety of spectral problems with remarkable accuracy.

Later the SPPS representation was extended to solutions of singular second order differential equations [10], of equations with polynomial dependence on the spectral parameter [21] and recently of linear differential equations of arbitrary order [18]. However the SPPS representation for the solutions of linear systems of differential equations was constructed only for Zakharov-Shabat system and a particular case of one dimensional Dirac system [21], in both cases by transforming the system into a certain Sturm-Liouville equation. Even though the general one dimensional Dirac system can be transformed into a Sturm-Liouville equation with potential polynomially dependent on the spectral parameter, see [3], and the result from [21] may be applied, we are not aware of such transformation for more complicated right-hand sides of the system. We opted for a different approach which allowed us to deal with arbitrary linear systems of two first order differential equations.

It is worth adding that the SPPS representation of solutions of one dimension Schrödinger equations allowed us to establish the mapping theorem for transmutation operators [9] which finally has led us to the development of two new methods, analytic approximation of transmutation operators (AATO) [20] and Neumann series of Bessel functions (NSBF) representation for the solutions [16] allowing one, in particular, to obtain hundreds of approximate eigenvalues with non deteriorating accuracy. The representation proposed in this paper opens possibility to extend the AATO and NSBF methods onto one-dimensional Dirac systems.

The paper is organized as follows. In Section 2 we consider a one-dimensional Dirac system whose right-hand side may contain arbitrary matrix-function coefficient at the spectral parameter. In Subsection 2.1 we introduce the formal powers starting from one non-vanishing particular solution corresponding to zero value of the spectral parameter. In Subsection 2.2 we show how the general solution of the system can be written in the terms of these formal powers (Theorem 2.2). In Subsection 2.3 we show the existence of a non-vanishing particular solution. In Subsection 2.4 we describe the spectral shift procedure. In Subsection 2.5 we extend the results onto discontinuos coefficients. In Section 3 we show that by simple transformation the general linear system of two first order differential equations can be reduced to the form covered in Section 2. In Section 4 we propose the general scheme for application of the SPPS representation to the solution of initial value and spectral problems, present numerical results for a particular Dirac system, show how a Sturm-Liouville spectral problem can be reformulated as a spectral problem for a Dirac system and discuss possible advantages of such problem reformulation.

2 The spectral parameter power series representation for solutions of generalized Dirac systems

Consider the following system

[TABLE]

or in the matrix form,

[TABLE]

where

[TABLE]

$p_{i},q,r_{ij}\in C[a,b]$, $i,j\in\{1,2\}$ are complex-valued functions of the real variable $x$, and $\lambda$ is an arbitrary complex constant. In the case when $R(x)=\begin{pmatrix}1&0\\ 0&1\\ \end{pmatrix}$ the system (2.2) is known as one-dimensional Dirac system [22], and in the case when $p_{1}=\bar{p}_{2}$, $q=0$ and $\begin{pmatrix}r_{11}&r_{12}\\ r_{21}&r_{22}\\ \end{pmatrix}=\begin{pmatrix}0&i\\ i&0\\ \end{pmatrix}$ the system (2.1) is known as Zakharov-Shabat system [1, 27].

2.1 A system of generalized formal powers

Suppose that the homogeneous system

[TABLE]

possesses a solution $Y=(f,g)^{T}$ such that both functions $f$ and $g$ are non-vanishing on $[a,b]$. From now on we will call such solution as non-vanishing solution of the homogeneous system (2.3). Let $x_{0}$ be a point from the segment $[a,b]$. Consider the following systems of functions defined by the recursive relations

[TABLE]

Similarly we use as the initial functions

[TABLE]

and define functions $\widetilde{Z}^{(n)}$, $\widetilde{X}^{(n+1)}$ and $\widetilde{Y}^{(n+1)}$, $n\geq 0$ using the same formulas (2.6)–(2.8) changing correspondingly all the functions $X^{(n)}$, $Y^{(n)}$ and $Z^{(n)}$ by $\widetilde{X}^{(n)}$, $\widetilde{Y}^{(n)}$ and $\widetilde{Z}^{(n)}$.

Example 2.1*.*

Following [17], let us consider a Sturm-Liouville equation

[TABLE]

where $p\in C^{1}[a,b]$, $q,r\in C[a,b]$ are complex-valued functions such that $p$ does not vanish on $[a,b]$. Suppose that a function $u_{0}$ is a non-vanishing solution corresponding to $\omega=0$. The following systems of functions were introduced in [17]

[TABLE]

and it was proved that the general solution of (2.11) has the form

[TABLE]

Let $\omega\neq 0$. Consider a function $v$ defined by

[TABLE]

Then equation (2.11) is equivalent to the system

[TABLE]

a particular case of (2.1) having $q=u_{0}^{\prime}/u_{0}$ and $p_{1}=p_{2}=0$. A non-vanishing particular solution of (2.14) corresponding to $\omega=0$ can be taken in the form $(u,v)^{T}=(u_{0},1/u_{0})^{T}$ and one can easily verify that the formal powers defined by (2.4)–(2.10) satisfy for all $n\geq 0$

[TABLE]

2.2 The SPPS representation

The following theorem gives the general solution of the system (2.2).

Theorem 2.2**.**

Suppose that the homogeneous system (2.3) possesses a solution $Y_{0}=(f,g)^{T}$ such that both functions $f$ and $g$ are non-vanishing on $[a,b]$. Then a general solution of the system (2.2) has the form

[TABLE]

where $c_{1}$ and $c_{2}$ are arbitrary complex constants and

[TABLE]

Here the formal powers $X^{(n)}$, $Y^{(n)}$, $\widetilde{X}^{(n)}$, $\widetilde{Y}^{(n)}$ are constructed starting with the particular solution $Y_{0}$ according to (2.4)–(2.10). Both series converge uniformly on $[a,b]$. The solutions $Y_{1}$ and $Y_{2}$ satisfy the following initial conditions:

[TABLE]

The proof of Theorem 2.2 requires several lemmas. In the first lemma we consider the nonhomogeneous system obtained from the left-hand side of (2.1) and construct a right-inverse operator for this system.

Lemma 2.3**.**

Under the conditions of Theorem 2.2 the solution of the nonhomogeneous system

[TABLE]

with the initial conditions

[TABLE]

where $h_{1,2}\in C[a,b]$ and $x_{0}\in[a,b]$, can be written in the form

[TABLE]

and

[TABLE]

Proof.

First we show how one can obtain the formulas (2.20) and (2.21) in the case when $h_{1,2}$, $q$, $p_{1,2}\in C^{1}[a,b]$ and the functions $p_{1}$ and $p_{2}$ are non-vanishing on $[a,b]$.

By differentiation and simple algebraic transformations one can verify that if the functions $u$ and $v$ satisfy (2.18), then these functions satisfy the following nonhomogeneous Sturm-Liouville equations

[TABLE]

Note that since the pair of functions $f$ and $g$ is a solution of the homogeneous system (2.3), the function $f$ is a solution of homogeneous part of equation (2.22), i.e., \left(\frac{f^{\prime}}{p_{2}}\right)^{\prime}+\Bigl{(}p_{1}-\bigl{(}\frac{q}{p_{2}}\bigr{)}^{\prime}-\frac{q^{2}}{p_{2}}\Bigr{)}f=0, and $g$ is a solution of the equation \left(\frac{g^{\prime}}{p_{1}}\right)^{\prime}+\Bigl{(}p_{2}+\bigl{(}\frac{q}{p_{1}}\bigr{)}^{\prime}-\frac{q^{2}}{p_{1}}\Bigr{)}g=0. We recall that if for the operator $L=\frac{d}{dx}p\frac{d}{dx}+q$ a non-vanishing function $u_{0}$ satisfying $Lu_{0}=0$ is known, then the operator $L$ possesses the Polya’s factorization $L=\frac{1}{u_{0}}\partial pu_{0}^{2}\partial\frac{1}{u_{0}}$, see, e.g., [25]. Using the Polya’s factorization and the Abel’s formula, general solutions of equations (2.22) and (2.23) can be written in the form

[TABLE]

and

[TABLE]

A pair $u$, $v$ of solutions of equations (2.22), (2.23) is a solution of the Cauchy problem (2.19) for the nonhomogeneous Dirac system (2.18) if and only if the functions $u$, $v$ satisfy

[TABLE]

Hence we obtain from (2.24)–(2.26) that $c_{11}=c_{21}=0$,

[TABLE]

Now we integrate by parts the expression f\cdot\Bigl{(}\frac{h_{2}}{p_{2}}\Bigr{)}^{\prime} in (2.24) and use the equality $f^{\prime}=qf+p_{2}g$ to obtain the formula (2.20):

[TABLE]

Similarly one obtains the second formula (2.21).

The validity of the formulas (2.20) and (2.21) without any additional requirements on the functions $q$, $p_{1,2}$ and $h_{1,2}$ can be checked directly. Indeed, using $g^{\prime}=-p_{1}f-qg$ we obtain

[TABLE]

It follows from the formulas (2.20) and (2.21) that

[TABLE]

Since the functions $f$ and $g$ satisfy (2.3), by simple algebraic manipulations we find that

[TABLE]

and integrating by parts in (2.29) we obtain

[TABLE]

Hence the expression in the brackets in (2.28) is equal to zero and we have verified the first equation in (2.18). Second equation can be verified similarly. The initial conditions (2.19) follow directly from the formulas (2.20) and (2.21). ∎

The next lemma in an analogue of the Abel’s formula. We construct the general solution of the homogeneous system (2.3) starting from a known particular solution.

Lemma 2.4**.**

Under the conditions of Theorem 2.2 the general solution of the homogeneous system (2.3) has the form

[TABLE]

where $\kappa:=f(x_{0})g(x_{0})$ and $c_{1}$ and $c_{2}$ are arbitrary complex constants.

Proof.

Similarly to the proof of Lemma 2.3 under the additional assumptions on $p_{1,2}$ and $q$ the formula (2.31) can be obtained from the expressions (2.24) and (2.25) considering $h_{1}=h_{2}=0$ and two initial conditions $\begin{pmatrix}u(x_{0})\\ v(x_{0})\\ \end{pmatrix}=\begin{pmatrix}f(x_{0})\\ 0\\ \end{pmatrix}$ (giving us the expression at $c_{1}$) and $\begin{pmatrix}u(x_{0})\\ v(x_{0})\\ \end{pmatrix}=\begin{pmatrix}0\\ g(x_{0})\\ \end{pmatrix}$ (giving us the expression at $c_{2}$).

Without any additional assumption verification of the fact that (2.31) is a solution of (2.3) can be done by the direct substitution. The solution given by the formula (2.31) is general since the expressions at $c_{1}$ and at $c_{2}$ are linearly independent which can be seen from their values at $x=x_{0}$. ∎

Lemma 2.5**.**

Under the conditions of Theorem 2.2 let us define

[TABLE]

where $\|\cdot\|$ denotes max-norm on $[a,b]$. Then the following estimates hold for the functions $X^{(n)}$, $Y^{(n)}$, $\widetilde{X}^{(n)}$, $\widetilde{Y}^{(n)}$, $Z^{(n)}$, $\widetilde{Z}^{(n)}$, $n\geq 0$.

[TABLE]

Proof.

The proof is straightforward by induction. For $n=0$, (2.34) follows directly from the definition of the constant $c$. Suppose (2.34) is true for some $n$. Then assuming $x\geq x_{0}$ we obtain from (2.6) and (2.35)

[TABLE]

and from (2.7) and (2.34) we obtain

[TABLE]

The proof for the case $x\leq x_{0}$ and for the functions $Y^{(n)}$, $\widetilde{X}^{(n)}$, $\widetilde{Y}^{(n)}$, $\widetilde{Z}^{(n)}$ is similar. ∎

Corollary 2.6**.**

Under the conditions of Lemma 2.5 the following estimates hold for any $n\geq 0$

[TABLE]

Proof.

The result follows from Lemma 2.5. Indeed, we have, for example, for the function $X^{(n)}$

[TABLE]

∎

Now we present the proof of Theorem 2.2.

Proof of Theorem 2.2.

It follows from Corollary 2.6 that both series in (2.16) converge uniformly on $[a,b]$ as well as the series of termwise derivatives, hence it is possible to apply the Dirac operator $B\frac{d}{dx}+P$ termwise to the series. Consider the second series in (2.16). As it follows from the definitions (2.4)–(2.8) and Lemmas 2.3 and 2.4, the functions $X^{(n)}$, $Y^{(n)}$, $n\geq 0$ satisfy

[TABLE]

where for $n=0$ the symbols $X^{(-1)}$ and $Y^{(-1)}$ appear only to unify the notation. Therefore,

[TABLE]

That is, the second series in (2.16) is a solution of the system (2.2). The proof for the first series in (2.16) is the same.

The linear independence of the solutions $(u_{1},v_{1})^{T}$ and $(u_{2},v_{2})^{T}$ can be obtained considering their values at $x=x_{0}$. Indeed, as it follows from the definitions (2.4)–(2.10),

[TABLE]

∎

Remark 2.7*.*

Taking into account the results of Example 2.1 one can see that Theorem 2.2 contains the SPPS representation from [17] as a particular case.

2.3 Construction of a non-vanishing particular solution

In this subsection we show how the general solution of the homogeneous system (2.3) can be obtained from Theorem 2.2 and discuss how one can select a solution $(u,v)^{T}$ of (2.3) such that both functions $u$ and $v$ are non-vanishing on $[a,b]$.

The system (2.3) can be rewritten either as

[TABLE]

or as

[TABLE]

Both systems can be considered as particular cases of a system of the type (2.1) taking $\lambda=-1$ and corresponding matrices $P$ and $R$. The homogeneous system associated to (2.39) possesses a non-vanishing particular solution $(u,v)^{T}=(1,1)^{T}$, while the homogeneous system associated to (2.40) possesses a non-vanishing particular solution $(u,v)^{T}=(\exp(\int q(s)\,ds),\exp(-\int q(s)\,ds))^{T}$. Hence, the general solution of the homogeneous system (2.3) can be obtained from Theorem 2.2 applied either to the system (2.39) or to the system (2.40).

In the case when all the coefficients $p_{1,2}$ and $q$ are real-valued functions, it is possible to construct a non-vanishing solution (complex-valued) explicitly. The following well-known proposition can be used.

Proposition 2.8**.**

Let the coefficients $p_{1,2}$ and $q$ of the system (2.3) are real-valued functions and $(u_{1},v_{1})^{T}$ and $(u_{2},v_{2})^{T}$ are two linearly independent real-valued solutions of the system (2.3). Then the linear combination

[TABLE]

is a non-vanishing solution of the system (2.3), i.e., both functions $u$ and $v$ do not have zeros on $[a,b]$.

Proof.

Suppose that at some point $x_{0}$ we have $u(x_{0})=0$. Since both functions $u_{1}$ and $u_{2}$ are real valued, $u_{1}(x_{0})=u_{2}(x_{0})=0$ which means that at the point $x_{0}$ $(u_{1}(x_{0}),v_{1}(x_{0}))^{T}=(0,v_{1}(x_{0}))^{T}=c\cdot(0,v_{2}(x_{0}))^{T}=c\cdot(u_{2}(x_{0}),v_{2}(x_{0}))^{T}$ for some constant $c$, a contradiction with the linear independency of the solutions. ∎

In the general situation when one or several of the coefficients of the system (2.3) may possess complex values, we are not aware of any explicit method of constructing a non-vanishing solution. However the situation is not that bad, at least a non-vanishing solution always exists and there are plenty of them.

Proposition 2.9**.**

Let $(u_{1},v_{1})^{T}$ and $(u_{2},v_{2})^{T}$ be two linearly independent solutions of (2.3). Then there exists a linear combination

[TABLE]

such that both functions $u$ and $v$ are non-vanishing on $[a,b]$.

Proof.

The proof is based on the Sard’s theorem and follows the proofs of Proposition 2.2 and Corollary 2.3 from [8], see also [17, Remark 5]. Consider the complex projective space $\mathbb{CP}^{1}$, i.e., the quotient of $\mathbb{C}^{2}\setminus\{0\}$ by the action of $\mathbb{C}^{*}$. Proposition 2.2 from [8] states the following. Let $I\subset\mathbb{R}$ be an interval. A differentiable map $f:I\to\mathbb{CP}^{1}$ cannot be surjective. And the proof is that the Sard’s theorem implies that the image $f(I)$ has measure zero.

In Corollary 2.3 [8] for two differentiable functions $y_{1}$ and $y_{2}$ which do not vanish simultaneously, authors consider a differentiable map

[TABLE]

If a linear combination $c_{1}y_{1}+c_{2}y_{2}$ vanishes at some point $x_{0}\in I$, then the determinant $\begin{vmatrix}y_{1}(x_{0})&-c_{2}\\ y_{2}(x_{0})&c_{1}\end{vmatrix}$ is equal to zero, which implies that $(y_{1}(x_{0}),y_{2}(x_{0}))$ is proportional to $(-c_{2},c_{1})$, hence $[-c_{2}:c_{1}]$ belongs to the image $f(I)$. Since the image $f(I)$ has a measure zero, the set of complex constants $[c_{1}:c_{2}]\in\mathbb{CP}^{1}$ for which the linear combination $c_{1}y_{1}+c_{2}y_{2}$ vanishes at some point, has measure zero.

Now we consider $y_{1}=u_{1}$, $y_{2}=u_{2}$ and the corresponding map $f_{1}$ given by (2.41). As was mentioned in the proof of Proposition 2.8, the functions $u_{1}$ and $u_{2}$ do not vanish simultaneously, hence the image $f_{1}([a,b])$ has measure zero. The same reasoning applies to $y_{1}=v_{1}$, $y_{2}=v_{2}$ and the corresponding map $f_{2}$. The union $f_{1}([a,b])\cup f_{2}([a,b])$ also has measure zero, hence for almost all points $[c_{1}:c_{2}]\in\mathbb{CP}^{1}$ both linear combinations $c_{1}u_{1}+c_{2}u_{2}$ and $c_{1}v_{1}+c_{2}v_{2}$ do not vanish at any point of the segment $[a,b]$. ∎

Remark 2.10*.*

As can be seen from the proof of Proposition 2.9, one can obtain a non-vanishing solution of (2.3) by taking two linearly independent solutions $(u_{1},v_{1})^{T}$ and $(u_{2},v_{2})^{T}$ of (2.3), choosing randomly a point on $\mathbb{CP}^{1}$, i.e., some complex numbers $c_{1}$ and $c_{2}$, and verifying if the linear combination $c_{1}(u_{1},v_{1})^{T}+c_{2}(u_{2},v_{2})^{T}$ vanishes on $[a,b]$. If not, we are done. If yes, repeating the process. Since the set of “bad” coefficients has measure zero, the non-vanishing solution in most cases will be obtained on the first try.

2.4 Spectral shift

Similarly to the Taylor series, an approximation of the solution given by a truncation of the series (2.16) is more accurate near the origin, while the accuracy deteriorates as the absolute value of the parameter $\lambda$ increases. In [10, 12, 17, 21] the spectral shift technique was introduced and successfully applied to improve the accuracy of the approximation for the large $\lambda$.

Under the assumption that the matrix $R$ is symmetric, i.e., $r_{12}\equiv r_{21}$ or equivalently

[TABLE]

application of the spectral shift is straightforward. Let a non-vanishing particular solution $(u,v)^{T}$ of (2.2) is known for some $\lambda=\lambda_{0}$. We can rewrite the system as

[TABLE]

which is again a system of the type (2.1) due to the assumption (2.42), and $(u,v)^{T}$ is a non-vanishing solution of the system (2.43) corresponding to $\lambda-\lambda_{0}=0$. Hence one can construct the formal powers and obtain the general solution of (2.43) (which also is a general solution of (2.1)) as the series with respect to the spectral parameter $\Lambda:=\lambda-\lambda_{0}$ using Theorem 2.2.

For the general case, when the condition (2.42) is not satisfied, we introduce the new unknown vector function $U$ by

[TABLE]

The system (2.43) for the new unknown $U$ takes the form

[TABLE]

a system which is again of the form (2.1). If $(u,v)^{T}$ is a non-vanishing solution for (2.1) corresponding to $\lambda=\lambda_{0}$, one can take $\frac{1}{w(x)}(u,v)^{T}$ as a non-vanishing solution for (2.44) corresponding to $\lambda-\lambda_{0}=0$ and construct the SPPS representation with respect to $\Lambda=\lambda-\lambda_{0}$.

2.5 Discontinuous coefficients

Following [28, Chapter 1] we can consider the system (2.1) with the coefficients $p_{i}$, $q$, $r_{ij}\in L^{1}(a,b)$. In such case a vector function $(u,v)^{T}$ is called a solution of (2.1) if both functions $u,v$ are absolutely continuous on $[a,b]$ and satisfy the system almost everywhere. With the slight modification to the proof the SPPS representation is valid for this case.

Suppose that functions $f,g\in\operatorname{AC}[a,b]$ are such that the following assumption holds.

[TABLE]

Note that it is sufficient for $f$ and $g$ to be non-vanishing on $[a,b]$. Then one can define the systems of functions $X^{(n)}$, $Y^{(n)}$, $Z^{(n)}$, $\widetilde{X}^{(n)}$, $\widetilde{Y}^{(n)}$ and $\widetilde{Z}^{(n)}$ by (2.4)–(2.10).

Lemma 2.11**.**

Let $f$ and $g$ be absolutely continuous functions satisfying (2.45). Define (to simplify the formulas we assume that $x>x_{0}$)

[TABLE]

Then the functions $X^{(n)}$, $Y^{(n)}$, $Z^{(n)}$, $\widetilde{X}^{(n)}$, $\widetilde{Y}^{(n)}$ and $\widetilde{Z}^{(n)}$, $n\geq 0$, are absolutely continuous and satisfy the following estimates

[TABLE]

The proof is by induction, similarly to the proof of Proposition 5 from [7]. We left the details to the reader.

Recall that the space $AC[a,b]$ of absolutely continuous functions coincides with the Sobolev space $W_{1}^{1}[a,b]$, hence the series $\sum_{n=0}^{\infty}u_{n}(x)$, where $u_{n}\in AC[a,b]$, converges to an absolutely continuous function if the series converges at some point $x_{0}\in[a,b]$ and the series of the derivatives $\sum_{n=0}^{\infty}u_{n}^{\prime}$ converges in $L^{1}(a,b)$ norm.

Now using estimates from Lemma 2.11 one can easily verify that the series $\sum_{n=0}^{\infty}\lambda^{n}X^{(n)}$, $\sum_{n=0}^{\infty}\lambda^{n}Y^{(n)}$ and $\sum_{n=0}^{\infty}\lambda^{n}Z^{(n)}$ are absolutely continuous functions and that Theorem 2.2 holds with the following change: $f$ and $g$ are required to be absolutely continuous functions satisfying (2.45).

3 General linear system

Consider a general linear system of two first order differential equations

[TABLE]

where $Y=(y_{1},y_{2})^{T}$ is the unknown vector-function and $\mathcal{P}$, $\mathcal{Q}$, $\mathcal{R}$ are $2\times 2$ matrices whose entries are continuous complex-valued functions. Assume additionally that $\det\mathcal{P}\neq 0$ for all $x\in[a,b]$.

Multiplying (3.1) by $B\mathcal{P}^{-1}$ we arrive to the system

[TABLE]

where $Q=B\mathcal{P}^{-1}\mathcal{Q}$ and $R=B\mathcal{P}^{-1}\mathcal{R}$. In general, the system (3.2) is not of the type (2.2) since the condition $\operatorname{tr}BQ(x)\equiv 0$ need not be satisfied. So we may proceed similarly to Subsection 2.4 and introduce new unknown vector function $U$ by

[TABLE]

for which the system (3.2) takes the form

[TABLE]

and since $B^{2}=-I$ one easily checks that

[TABLE]

hence the system (3.3) is of the type (2.2) and Theorem 2.2 can be applied to it.

4 Numerical illustration

4.1 General scheme and implementation details

The general scheme of application of the SPPS representation to the approximate solution of initial value and spectral problems for the system (2.1) is similar to that for the SPPS representation for the Sturm-Liouville equation, see [17], [19], [21]. Consider an initial value problem

[TABLE]

and a spectral problem given by the following boundary conditions

[TABLE]

where $y_{1}$, $y_{2}$, $a_{1}$, $a_{2}$, $b_{1}$, $b_{2}$ are some complex constants satisfying $|a_{1}|+|a_{2}|\neq 0$ and $|b_{1}|+|b_{2}|\neq 0$. We would like to stress out that more complicated problems like having mixed or spectral parameter dependent boundary conditions can be treated as well.

Construct a non-vanishing particular solution of (2.3) as described in Subsection 2.3. According to Remark 2.10 one can take two linearly independent solutions $(u_{1},v_{1})^{T}$ and $(u_{2},v_{2})^{T}$ satisfying $u_{1}(x_{0})\neq 0$, take some finite set of complex constants $c_{1},\ldots,c_{K}$ (they can be taken, e.g., by randomly choosing modulus $\rho_{k}$ and phase $\varphi_{k}$ in the polar representation $c_{k}=\rho_{k}e^{i\varphi_{k}}$) and select as the non-vanishing solution the one having the least value of the expression

[TABLE] 2. 2.

Calculate the formal powers $X^{(n)}$, $Y^{(n)}$, $\widetilde{X}^{(n)}$ and $\widetilde{Y}^{(n)}$, $n=0,\ldots,N$ according to (2.4)–(2.10). The number $N$ may be estimated either using the bounds from Lemma 2.5 or simply by verifying that

[TABLE]

is sufficiently small. E.g., it is equal to zero in machine-precision arithmetic. Here $\|\cdot\|$ denotes max-norm on $[a,b]$. 3. 3.

For the solution of the initial value problem (4.1) one calculates approximate solutions $\widetilde{Y}_{1}$ and $\widetilde{Y}_{2}$ using truncated sums from (2.16).

In the particular case when $x_{0}=a$ the solution of (4.1) due to (2.17) is given by

[TABLE]

In the general case, one finds the constants $c_{1}$ and $c_{2}$ in (2.15) by solving the linear system of equations

[TABLE] 4. 4.

For the solution of the spectral problem (4.2) in the particular case $x_{0}=a$ note that a solution satisfying the first boundary condition in (4.2) is given (up to a multiplicative constant) by

[TABLE]

This solution satisfies the second boundary condition in (4.2) iff. the following characteristic equation is satisfied:

[TABLE]

The function $\Xi$, known as characteristic function of the spectral problem, is analytic. By truncating the series representations (2.16) for the functions $u_{1}$, $u_{2}$, $v_{1}$ and $v_{2}$ one obtains approximate characteristic function which results to be a polynomial. Its (complex) roots closest to the origin approximate the exact eigenvalues, while more distant roots result to be spurious. See [10, Section 7.2] for further discussion on how these spurious roots can be discarded.

For the general case one considers the general solution given by (2.15) and substitutes it into boundary conditions (4.2). The existence of non-trivial solution is equivalent to the following condition:

[TABLE]

Truncating the series representations for $u_{1}$, $u_{2}$, $v_{1}$ and $v_{2}$ in (4.4) on obtains a polynomial approximating the characteristic function of the problem. 5. 5.

If more than few closest to zero eigenvalues are needed, one may apply several spectral shift procedures described in Subsection 2.4. Since for the one-dimensional Dirac system the distance between consequent eigenvalues remains bounded for all eigenvalues, see [22, Chap. 7, §2], the following simple recipe showed to deliver excellent results. We use the SPPS representation to find the eigenvalues $\lambda_{0}$ and $\lambda_{\pm 1}$ and general solutions corresponding to $\lambda_{\pm 1}$. Than we look for complex coefficients which give us non-vanishing solutions corresponding to $\lambda_{\pm 1}$ using the same criteria as in (4.3). Now we use $\lambda_{1}$ and $\lambda_{-1}$ as centers for the spectral shift procedure and use the corresponding non-vanishing solutions to calculate the formal powers. Resulting SPPS representations give us $\lambda_{2}$ and $\lambda_{-2}$ and corresponding general solutions. And so on, having found $\lambda_{n}$ and $\lambda_{-n}$ and corresponding non-vanishing solutions we use them as the new centers for the spectral shift procedure (the total step being $\lambda_{n}-\lambda_{n-1}$ and $\lambda_{-n}-\lambda_{-(n-1)}$, bounded quantity as $n\to\infty$) until required number of eigenvalues be find.

We would like to emphasize that all steps of the proposed scheme can be realized numerically, there is no need for the integrals in (2.4)–(2.10) to be available in the closed form. We refer the reader to [10], [19] and [20] for additional details and only mention that in the following examples all the functions involved were represented by their values on the uniform mesh and Newton-Cotes 6 point rule was used for indefinite integration. All computations were performed in double machine precision in Matlab 2017.

4.2 Example: spectral problem for a Dirac system

Consider the following spectral problem (Example 3.4 from [2])

[TABLE]

with the boundary conditions

[TABLE]

The characteristic function for this problem can be written in the terms of Airy functions, see [2], allowing one to compute exact eigenvalues with any desired precision using, e.g., Wolfram Mathematica.

We used 2001 points mesh to represent all the functions involved in this example and computed the formal powers for $n\leq 100$. In Table 1 we present the absolute errors of the approximate eigenvalues obtained either directly from the SPPS representation (2.16) or using the spectral shift procedure. In the latter case the eigenvalues $\lambda_{n}$ for $|n|\leq 100$ were computed. Even direct application of the SPPS representation delivers more accurate values than those reported in [2] requiring only 0.1sec of computation time, while the spectral shift technique allows one to obtain two hundreds eigenvalues with a good accuracy in about 20 seconds.

4.3 Example: application to Sturm-Liouville spectral problems

In Example 2.1 we showed how a Sturm-Liouville equation can be transformed into a one dimensional Dirac system. A spectral problem for equation (2.11) can be transformed into a spectral problem for the system (2.14) as well. Indeed, consider a boundary condition of the form

[TABLE]

Due to (2.13) we have that

[TABLE]

hence (4.7) is equivalent to the following (spectral parameter dependent) boundary condition

[TABLE]

The general scheme presented in Subsection 4.1 can be applied to spectral parameter dependent boundary conditions of the form (4.8) with minimal modifications. Note that rewriting the equation $(p(x)u^{\prime})^{\prime}+q(x)u=0$ as an equivalent system

[TABLE]

one can apply the results of Subsection 2.3 to construct a non-vanishing particular solution $u_{0}$ of (2.11) as well.

One possible advantage of reformulating a Sturm-Liouville spectral problem as an equivalent Dirac system and applying the proposed SPPS representation instead of the SPPS representation from [17] consists in the following. The eigenvalues $\lambda_{n}=\omega_{n}^{2}$ of the Sturm-Liouville problem grow as $\frac{\pi}{b-a}n^{2}$ as $n\to\infty$. However the spectral shift technique works best if the change of the spectral parameter remains bounded on each step. As a result, obtaining $n$ eigenvalues requires $O(n^{2})$ steps using the SPPS representation from [17]. In contrary, Dirac-based approach requires only $O(n)$ steps greatly reducing computation time if one needs large number of eigenvalues. We would like to point out that the method proposed in [16] is better suited for computing large sets of eigenvalues, nevertheless the SPPS representation is simpler and still being used for numerous applications, see, e.g., [4], [5], [6], [23], [26].

Consider the following spectral problem (the first Paine problem, [24], see also [20, Example 7.4]) to illustrate this advantage numerically.

[TABLE]

We computed approximate eigenvalues for this problem using the Dirac system reformulation as well as directly the SPPS representation from [17]. For the latter we applied two different spectral shift strategies, to the nearest new eigenvalue and constant step size. The first one requires less steps but is known to fail eventually due to increasing gaps between consecutive eigenvalues, the second requires more steps but allowed us to obtain more accurate results previously, see [10], [19].

We used 50001 points mesh to represent all the functions involved in this example and computed the formal powers for $n\leq 100$. The spectral shift for the constant step was taken equal to $5$. Note that such large number of mesh points was taken in order to avoid integration errors and to be able to illustrate SPPS-related behavior, compare to [21, Example 6.1]. In Table 2 we present the absolute errors of the approximate eigenvalues. As one can observe, the performance of all three approaches was comparable for the first 50 eigenvalues. For larger index eigenvalues the accuracy of the Dirac system based approach remains essentially the same while the accuracy of the eigenvalues obtained using directly the SPPS representation from [17] started to deteriorate (surprisingly more rapidly when the uniform step size was used).

Acknowledgements

The authors acknowledge the support from CONACYT, Mexico via the project 222478. N. Gutiérrez Jiménez would like to express his gratitude to the Mathematical department of Cinvestav where he completed the PhD program (the presented paper contains part of the obtained results) and to CONACYT, Mexico for the financial support making it possible.