On the High-Order Quasi Exactly Solvable Differential Equations

Abstract

In this paper, we present a new method for solving a class of high-order quasi exactly solvable ordinary differential equations. With this method, the computed solution is expressed as a linear combination of the canonical polynomials associated with the given differential operator. An iterative algorithm summarizing the procedure is presented and its efficiency is demonstrated through considering two applied problems.

Share and Cite:

Alzanki, T. , Shaaban, M. and El-Daou, M. (2019) On the High-Order Quasi Exactly Solvable Differential Equations. American Journal of Computational Mathematics, 9, 234-250. doi: 10.4236/ajcm.2019.94018.

1. Introduction

Let us consider the ordinary differential operator D of order $\nu \ge 2$,

$\left(Du\right)\left(x\right):=\underset{i=0}{\overset{\nu }{\sum }}\text{ }\text{ }{A}_{i}\left(x\right)\frac{{\text{d}}^{i}u}{\text{d}{x}^{i}}$ (1)

where $\left\{{A}_{i}\left(x\right);i=0,1,\cdots ,\nu \right\}$ is a set of $\nu +1$ polynomials, $\text{deg}\left[{A}_{i}\left(x\right)\right]=i+p$ for all $i=0,1,\cdots ,\nu$, where p is a prescribed nonnegative integer called the height of D, with

${A}_{i}\left(x\right):=\underset{j=0}{\overset{i+p}{\sum }}\text{ }{A}_{ij}{x}^{j},\text{\hspace{0.17em}}\text{\hspace{0.17em}}i=0,1,\cdots ,\nu .$ (2)

This paper is concerned with the solution of the following problem: Given $n\ge 0$ and $p\ge 0$, construct two polynomials, $y\left(x\right)=\underset{j=0}{\overset{n}{\sum }}{y}_{j}{x}^{j}$ of degree n and $\lambda \left(x\right)=\underset{j=0}{\overset{p}{\sum }}{\lambda }_{j}{x}^{j}$ of degree p such that the pair $\left\{y\left(x\right),\lambda \left(x\right)\right\}$ satisfies the following ordinary differential equation exactly

$\left(Dy\right)\left(x\right)=\lambda \left(x\right)y\left(x\right).$ (3)

When $\nu =2$, Equation (3) is called quasi exactly solvable (QES). This class of QES problems has applications in various fields of engineering, chemistry and quantum mechanics. Many different techniques to solve QES equations are reported in the literature: Among these are the Functional Ansatz Method, Constraint polynomial approach, asymptotic iteration method and Lie algebraic method (see  -  ). The case $\nu =2$ was also discussed very recently in , where the authors developed a new approach based on a special set of polynomials associated with the differential operator D called canonical polynomials. The main objective of this paper is to extend that canonical polynomials approach to solve equations of the form (3) with arbitrary order $\nu \ge 2$. More precisely, we present a procedure to construct a pair of polynomials based on the canonical polynomial associated with D. While the existing method for solving QES requires the solution of a nonlinear algebraic system with dimensions depending on the desired degree of y, the canonical polynomial approach presented in  requires a nonlinear algebraic system of dimensions depending on p only. This advantage is due to the fact that the sequence of canonical polynomials enjoys the permanence characteristic .

The canonical polynomials (to be explained shortly) appeared for the first time in  wherein Lanczos developed an efficient method, called the Tau method, to approximate the exact solution of differential equations in terms of a finite number of canonical polynomials. Later on, the concept of the canonical polynomial was generalized in  to develop a recursive approach of the Tau method that can apply to more complex differential equations. And it was due to the computational efficiency of the canonical polynomials that makes the Tau method more competitive compared to other existing approximation methods (more details can be found in  -  ).

Section 2 will concentrate on the construction of the canonical polynomials associated with the $\nu$ th differential operator (1) and on their computation. In Section 3 we present an algorithm that allows to obtain the pair of polynomials $\left\{y\left(x\right),\lambda \left(x\right)\right\}$ in an effective way. Two examples confirming our results are discussed in Section 4.

2. The Canonical Polynomials

Let D be the differential operator defined in (1). In this section we recall the main features of the canonical polynomials associated with D (see  ), and we give an algorithm for computing them. First rewrite (3):

$\left(Dy\right)\left(x\right)-\lambda \left(x\right)y\left(x\right)=\left(\underset{i=1}{\overset{\nu }{\sum }}\text{ }{A}_{i}\left(x\right)\frac{{\text{d}}^{i}y}{\text{d}{x}^{i}}\right)+\left({A}_{0}\left(x\right)-\lambda \left(x\right)\right)y\left(x\right)=0$

That is

$\left({D}^{*}y\right)\left(x\right):=\underset{i=0}{\overset{\nu }{\sum }}\text{ }{A}_{i}^{*}\left(x\right)\frac{{\text{d}}^{i}y}{\text{d}{x}^{i}}=0$

where ${A}_{0}^{*}\left(x\right):={A}_{0}\left(x\right)-\lambda \left(x\right)$ and ${A}_{i}^{*}\left(x\right)={A}_{i}\left(x\right),i=1,2,\cdots ,\nu$. So, for the sake of simplicity, we shall hide the asterisk "*" and carry out the analysis for

$\left(Dy\right)\left(x\right):=\underset{i=0}{\overset{\nu }{\sum }}\text{ }{A}_{i}\left(x\right)\frac{{\text{d}}^{i}y}{\text{d}{x}^{i}}=0$ (4)

keeping in mind that ${A}_{0}\left(x\right)$ involves the unknown coefficients of $\lambda \left(x\right)$.

Definition 1. For any integer $k\ge 0$, ${Q}_{k}^{*}\left(x\right)$ is called a kth canonical function of D if $D{Q}_{k}^{*}={x}^{k}$.

The following notation will enable us to formulate the next theorem:

${A}_{rj}^{\left(k\right)}:={A}_{rj}{\delta }_{kr}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\text{where}\text{\hspace{0.17em}}\text{ }{\delta }_{kr}=\left(\begin{array}{l}\frac{k!}{\left(k-r\right)!}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{ }k\ge r\\ 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{ }0\le k (5)

${\alpha }_{j}^{\left(k\right)}:=\underset{i=0}{\overset{\nu }{\sum }}\text{ }\text{ }{A}_{i,i+j}^{\left(k\right)}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}{\beta }_{j}^{\left(k\right)}:=\underset{i=j}{\overset{\nu }{\sum }}\text{ }\text{ }{A}_{i,i-j}^{\left(k\right)}$ (6)

Theorem 1. Under the above assumptions and notation, the canonical functions associated with the differential operator (4) are formally generated by the recursion:

${Q}_{k+p}^{*}=\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left[{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}^{*}\right],\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=0,1,2,\cdots$ (7)

provided ${\alpha }_{p}^{\left(k\right)}\ne 0$.

In particular, if $p=0$ and ${\alpha }_{0}^{\left(k\right)}\ne 0$, then

${Q}_{k}^{*}=\frac{1}{{\alpha }_{0}^{\left(k\right)}}\left[{x}^{k}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}^{*}\right],\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=0,1,2,\cdots$

Proof. For $k=0,1,2,\cdots$

$\begin{array}{c}D{x}^{k}=\underset{r=0}{\overset{\nu }{\sum }}\left[\underset{j=0}{\overset{r+p}{\sum }}\text{ }\text{ }{A}_{rj}{x}^{j}\right]k\left(k-1\right)\left(k-2\right)\cdots \left(k-r+1\right){x}^{k-r}\\ =\underset{r=0}{\overset{\nu }{\sum }}\left[\underset{j=0}{\overset{r+p}{\sum }}\text{ }\text{ }{A}_{rj}{x}^{j}\right]\frac{k!}{\left(k-r\right)!}{x}^{k-r}=\underset{r=0}{\overset{\nu }{\sum }}\underset{j=0}{\overset{r+p}{\sum }}\left({A}_{rj}{\delta }_{kr}\right){x}^{j+k-r}\\ ={x}^{k}\underset{r=0}{\overset{\nu }{\sum }}\underset{j=0}{\overset{r+p}{\sum }}\text{ }\text{ }{A}_{rj}^{\left(k\right)}{x}^{j-r}={x}^{k}\left\{\underset{j=0}{\overset{p}{\sum }}\left(\underset{i=0}{\overset{\nu }{\sum }}\text{ }\text{ }{A}_{i,i+j}^{\left(k\right)}\right){x}^{j}+\underset{j=1}{\overset{\nu }{\sum }}\left(\underset{i=j}{\overset{\nu }{\sum }}\text{ }\text{ }{A}_{i,i-j}^{\left(k\right)}\right){x}^{-j}\right\}\\ =\underset{j=0}{\overset{p}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{x}^{k+j}+\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{x}^{k-j}={\alpha }_{p}^{\left(k\right)}{x}^{k+p}+\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{x}^{k+j}+\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{x}^{k-j}\\ ={\alpha }_{p}^{\left(k\right)}{x}^{k+p}+\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}D{Q}_{k+j}^{*}+\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}D{Q}_{k-j}^{*}.\end{array}$

Since D is linear, the latter yields:

$D\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}^{*}\right\}={\alpha }_{p}^{\left(k\right)}{x}^{k+p}.$ (8)

If ${\alpha }_{p}^{\left(k\right)}=0$, then ${x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}{\beta }_{j}^{\left(k\right)}{Q}_{k-j}^{*}$ is an exact solution.

If ${\alpha }_{p}^{\left(k\right)}\ne 0$, then we obtain the desired formula for $\left\{{Q}_{k}^{*}\right\}$ :

${Q}_{k+p}^{*}=\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}^{*}\right\}.$

In particular, if $p=0$ then

$\begin{array}{c}D{x}^{k}=\underset{r=0}{\overset{\nu }{\sum }}\left[\underset{j=0}{\overset{r}{\sum }}\text{ }\text{ }{A}_{rj}{x}^{j}\right]\frac{k!}{\left(k-r\right)!}{x}^{k-r}\\ ={\alpha }_{0}^{\left(k\right)}{x}^{k}+\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{x}^{k-j}\\ ={\alpha }_{0}^{\left(k\right)}{x}^{k}+\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}D{Q}_{k-j}^{*}.\end{array}$

If ${\alpha }_{0}^{\left(k\right)}=0$, then ${x}^{k}-\underset{j=1}{\overset{\nu }{\sum }}{\beta }_{j}^{\left(k\right)}{Q}_{k-j}^{*}$ is an exact solution.

If ${\alpha }_{0}^{\left(k\right)}\ne 0$, then

${Q}_{k}^{*}=\frac{1}{{\alpha }_{0}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}^{*}\right\}.$

This completes the proof.

For illustration, when $k=0$, Equation (7) gives:

${Q}_{p}^{*}=\frac{1}{{\alpha }_{p}^{\left(0\right)}}\left\{1-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(0\right)}{Q}_{j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(0\right)}{Q}_{0-j}^{*}\right\}=\frac{1}{{\alpha }_{p}^{\left(0\right)}}\left\{1-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(0\right)}{Q}_{j}^{*}\right\}.$

When $k=1$,

$\begin{array}{c}{Q}_{p+1}^{*}=\frac{1}{{\alpha }_{p}^{\left(1\right)}}\left\{x-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(1\right)}{Q}_{j+1}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(1\right)}{Q}_{1-j}^{*}\right\}\\ =\frac{1}{{\alpha }_{p}^{\left(1\right)}}\left\{x-{\alpha }_{p-1}^{\left(1\right)}{Q}_{p}^{*}-\underset{j=0}{\overset{p-2}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(1\right)}{Q}_{j+1}^{*}-{\beta }_{1}^{\left(1\right)}{Q}_{0}^{*}\right\}.\end{array}$

Proceeding this way, we find that for $k=0,1,2,3,\cdots$

${Q}_{k+p}^{*}\in \text{span}\left\{1,x,\cdots ,{x}^{k}\right\}\cup \text{span}\left\{{Q}_{0}^{*},{Q}_{1}^{*},\cdots ,{Q}_{p-1}^{*}\right\}.$

We are able now to formulate one of the main results of this paper:

Theorem 2. For all $k=0,1,2,\cdots$, each ${Q}_{k+p}^{*}$ can be written in the form

${Q}_{k+p}^{*}\left(x\right)={Q}_{k+p}\left(x\right)+{R}_{k+p}$ (9)

where ${Q}_{k+p}\left(x\right)$ is a polynomial of degree k, called a canonical polynomial associated with D, and generated by the self starting recursive formula:

$\begin{array}{l}{Q}_{1}={Q}_{2}=\cdots ={Q}_{p-1}=0\\ \text{and}\text{\hspace{0.17em}}\text{ }{Q}_{k+p}=\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}\right\},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(\text{if}\text{\hspace{0.17em}}{\alpha }_{p}^{\left(k\right)}\ne 0\right)\end{array}$ (10)

and where ${R}_{k+p}$ is a linear combination of the undefined canonical polynomials $\left\{{Q}_{0}^{*},{Q}_{1}^{*},\cdots ,{Q}_{p-1}^{*}\right\}$, called residual, and written as ${R}_{k+p}=\underset{r=0}{\overset{p-1}{\sum }}{\rho }_{k+p}^{\left(r\right)}{Q}_{r}^{*}$,

where $\left\{{\rho }_{k+p}^{\left(r\right)};k\ge 0\right\}$ are sequences of constants given by the self starting recursion

$\begin{array}{l}{\rho }_{j}^{\left(i\right)}=\left\{\begin{array}{l}1\text{ }\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{ }i=j\\ 0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{otherwise}\text{ }\end{array}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\text{with}\text{\hspace{0.17em}}\text{ }0\le i,j\le p-1\\ \text{and}\text{\hspace{0.17em}}\text{ }{\rho }_{k+p}^{\left(r\right)}=\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{\rho }_{k+j}^{\left(r\right)}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{\rho }_{k-j}^{\left(r\right)}\right\},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(\text{if}\text{\hspace{0.17em}}{\alpha }_{p}^{\left(k\right)}\ne 0\right)\end{array}$ (11)

Proof. This follows by an induction argument once (9) is inserted in (10) and the terms are rearranged:

$\begin{array}{c}{Q}_{k+p}^{*}=\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}\left({Q}_{k+j}+{R}_{k+j}\right)-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}\left({Q}_{k-j}+{R}_{k-j}\right)\right\}\\ =\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{R}_{k+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{R}_{k-j}\right\}\\ =\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{R}_{k+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{R}_{k-j}\right\}\\ =\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}\right\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{R}_{k+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{R}_{k-j}\right\}\end{array}$

$\begin{array}{l}=\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}\right\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}\left(\underset{r=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\rho }_{k+j}^{\left(r\right)}{Q}_{r}^{*}\right)-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}\left(\underset{r=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\rho }_{k-j}^{\left(r\right)}{Q}_{r}^{*}\right)\right\}\\ =\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{Q}_{k+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}\right\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\underset{r=0}{\overset{p-1}{\sum }}\frac{1}{{\alpha }_{p}^{\left(k\right)}}\left\{-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(k\right)}{\rho }_{k+j}^{\left(r\right)}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{\rho }_{k-j}^{\left(r\right)}\right\}{Q}_{r}^{*}\\ \equiv {Q}_{k+p}+{R}_{k+p}\end{array}$

yielding Equations (10) and (11) as required.

3. Construction of Solution

This section is concerned with the construction of the two polynomials $\left\{y\left(x\right),{A}_{0}\left(x\right)\right\}$ that satisfy Equation (4).

Theorem 3. The above notation and assumptions hold. Suppose that the

coefficients $\left\{{A}_{0j},j=0,1,\cdots ,p\right\}$ of ${A}_{0}\left(x\right)=\underset{j=0}{\overset{p}{\sum }}\text{ }\text{ }{A}_{0j}{x}^{j}$ satisfy the following system

consisting of $p+1$ algebraic equations:

${A}_{0,p}=-n!\underset{i=1}{\overset{\nu }{\sum }}\frac{{A}_{i,i+p}}{\left(n-i\right)!}$ (12)

$\underset{j=0}{\overset{p-1}{\sum }}\left[\underset{i=0}{\overset{\nu }{\sum }}\frac{{A}_{i,i+j}}{\left(n-i\right)!}\right]{\rho }_{n+j}^{\left(r\right)}+\underset{j=1}{\overset{\nu }{\sum }}\left[\underset{i=j}{\overset{\nu }{\sum }}\frac{{A}_{i,i-j}}{\left(n-i\right)!}\right]{\rho }_{n-j}^{\left(r\right)}=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}r=0,1,\cdots ,p-1$ (13)

where $\left\{{\rho }_{k}^{\left(r\right)}\right\}$ are given by (11). Then

${Y}_{n}\left(x\right):={x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{Q}_{n+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}$ (14)

is an exact polynomial solution of Equation (4) where $\left\{{\alpha }_{j}^{\left(n\right)},{\beta }_{j}^{\left(n\right)}\right\}$ are parameters determined in terms of ${\left\{{A}_{ij},j=0,1,\cdots ,p\right\}}_{i=0}^{\nu }$ as defined in (5)-(6), and $\left\{{Q}_{k}\right\}$ is a sequence of canonical polynomials associated with D and recursively generated by (10)

Proof. Let $n\ge 1$. Setting $k=n$ in Equation (8) we get

$D\left\{{x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{Q}_{n+j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}^{*}\right\}={\alpha }_{p}^{\left(n\right)}{x}^{n+p}.$ (15)

If condition (12) holds, then (5) implies that ${\alpha }_{p}^{\left(n\right)}=0$ and consequently the right hand side of (15) vanishes:

$D\left\{{x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{Q}_{n+j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}^{*}\right\}=0.$

In other words,

${Y}_{n}^{*}\left(x\right):={x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{Q}_{n+j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}^{*}$

becomes an exact solution, but not necessarily an exact polynomial due to the appearance of the undefined canonical polynomials $\left\{{Q}_{0}^{*},{Q}_{1}^{*},\cdots ,{Q}_{p-1}^{*}\right\}$. However, in order to be an exact polynomial, ${Y}_{n}^{*}\left(x\right)$ must be independent of the p

undefined canonical functions ${\left\{{Q}_{i}^{*}\right\}}_{i=0}^{p-1}$. This can be achieved by an appropriate

adjustment of the p coefficients $\left\{{A}_{00},{A}_{01}j,\cdots ,{A}_{0,p-1}\right\}$ of ${A}_{0}\left(x\right)$ as explained next. Using (9) in ${Y}_{n}^{*}\left(x\right)$ we can write:

$\begin{array}{c}{Y}_{n}^{*}\left(x\right)={x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{Q}_{n+j}^{*}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}^{*}\\ ={x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}\left({Q}_{n+j}+{R}_{n+j}\right)-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}\left({Q}_{n-j}+{R}_{n-j}\right)\\ =\left({x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{Q}_{n+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{R}_{n+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{R}_{n-j}\right)\end{array}$

$=\left({x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{Q}_{n+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}\right)$ (16)

$\text{ }+\underset{r=0}{\overset{p-1}{\sum }}\left(-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{\rho }_{n+j}^{\left(r\right)}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{\rho }_{n-j}^{\left(r\right)}\right){Q}_{r}^{*}$ (17)

Working out the coefficients of ${\left\{{Q}_{r}^{*}\right\}}_{r=0}^{p-1}$ in (17) we find that due to Equation (13) we get:

$\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{\rho }_{n+j}^{\left(r\right)}+\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{\rho }_{n-j}^{\left(r\right)}=0.$

Thus ${Y}_{n}^{*}\left(x\right)$ reduces to the polynomial (16)

${Y}_{n}\left(x\right)={x}^{n}-\underset{j=0}{\overset{p-1}{\sum }}\text{ }\text{ }{\alpha }_{j}^{\left(n\right)}{Q}_{n+j}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}$.

The following corollary follows immediately from the previous theorem:

Corollary 4. If $p=0$ and

${A}_{00}=-n!\underset{i=1}{\overset{\nu }{\sum }}\frac{{A}_{ii}}{\left(n-i\right)!}$ (18)

then

${Y}_{n}\left(x\right)={x}^{n}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(n\right)}{Q}_{n-j}$

is an exact polynomial solution of Equation (4) where

${Q}_{k}=\frac{1}{{\alpha }_{0}^{\left(k\right)}}\left\{{x}^{k}-\underset{j=1}{\overset{\nu }{\sum }}\text{ }\text{ }{\beta }_{j}^{\left(k\right)}{Q}_{k-j}\right\},\text{\hspace{0.17em}}\text{\hspace{0.17em}}k=0,1,2,\cdots$ (19)

(note that ${\beta }_{j}^{\left(k\right)}=0$ if $j and therefore all ${Q}_{k}$ ’s are defined).

Computational aspects

For computational purposes, one can reduce the height of D from p to zero by differentiating (4) p times. This is due to the following trivial identity:

$\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left(UV\right)=\underset{k=0}{\overset{p}{\sum }}\text{ }\text{ }{C}_{p}^{k}{U}^{\left(k\right)}{V}^{\left(p-k\right)},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{C}_{p}^{k}=\frac{p!}{k!\left(p-k\right)!}$ (20)

Applying this identity to (4) we get

$\begin{array}{l}\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left[\left(Dy\right)\left(x\right)\right]:=\underset{i=0}{\overset{\nu }{\sum }}\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left({A}_{i}\left(x\right)\frac{{\text{d}}^{i}y}{\text{d}{x}^{i}}\right)\\ =\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left[{A}_{\nu }{y}^{\left(\nu \right)}\right]+\underset{k=1}{\overset{\nu -1}{\sum }}\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left({A}_{i}\left(x\right)\frac{{\text{d}}^{i}y}{\text{d}{x}^{i}}\right)+\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left[{A}_{0}y\right]\end{array}$ (21)

From (20),

$\begin{array}{c}\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left[{A}_{0}y\right]=\underset{k=0}{\overset{p}{\sum }}\text{ }\text{ }{C}_{p}^{k}{A}_{0}^{\left(k\right)}{y}^{\left(p-k\right)}=\underset{k=0}{\overset{p-1}{\sum }}\text{ }\text{ }{C}_{p}^{k}{A}_{0}^{\left(k\right)}{y}^{\left(p-k\right)}+{A}_{0}^{\left(p\right)}y\\ =\underset{k=0}{\overset{p-1}{\sum }}\text{ }\text{ }{C}_{p}^{k}{A}_{0}^{\left(k\right)}{y}^{\left(p-k\right)}+\left(p!{A}_{0,p}\right)y\end{array}$

$\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left[{A}_{\nu }{y}^{\left(\nu \right)}\right]=\underset{k=0}{\overset{p}{\sum }}\text{ }\text{ }{C}_{p}^{k}{A}_{\nu }^{\left(k\right)}{y}^{\left(\nu +p-k\right)}={A}_{\nu }{y}^{\left(\nu +p\right)}+\underset{k=1}{\overset{p}{\sum }}\text{ }\text{ }{C}_{p}^{k}{A}_{\nu }^{\left(k\right)}{y}^{\left(\nu +p-k\right)}$

Inserting the later in (18) we get

$\left({D}^{\left(p\right)}y\right)\left(x\right)={A}_{\nu }{y}^{\left(\nu +p\right)}+\left[\underset{k=1}{\overset{\nu -1}{\sum }}\frac{{\text{d}}^{p}}{\text{d}{x}^{p}}\left({A}_{i}\left(x\right){y}^{\left(i\right)}\right)+\underset{k=0}{\overset{p-1}{\sum }}\text{ }\text{ }{C}_{p}^{k}{A}_{0}^{\left(k\right)}{y}^{\left(p-k\right)}\right]+\left(p!{A}_{0,p}\right)y$ (22)

which is a differential operator of order $\nu +p$ with height 0. Therefore we can apply our results to (22) and reconstruct the solution of the original problem by an antiderivative process. This will reduce the computation cost because the residual subspace of the new operator will be 0.

4. Applications

In the section we solve two applied problems by means of Algorithm ((12)-(13)-(14)-(18)-(19)) formulated in Theorem 3:

Example 1. Modified Manning potential with parameters. Let us consider the Schrodinger’s equation

$-\frac{{\text{d}}^{2}\psi }{\text{d}{x}^{2}}+V\left(x\right)\psi \left(x\right)=E\psi$ (23)

where the potential $V\left(x\right)$ is given by

$V\left(x\right)=-{v}_{1}\text{ }{\text{sech}}^{6}x-{v}_{2}\text{ }{\text{sech}}^{4}x-{v}_{3}\text{ }{\text{sech}}^{2}x$ where ${v}_{1}$, ${v}_{2}$ are given constant parameters and ${v}_{3}$ is an unknown parameter and E is the unknown eigenvalue. We wish to compute E and ${v}_{3}$. This potential describes a double-well potential

whenever ${v}_{1}>0$, ${v}_{2}<0$, ${v}_{3}>0$ and $-\frac{{v}_{3}}{2{v}_{2}}<1$ which was discussed in ( ,  ).

Equation (23) can be written as a 2nd order QES in the form (4) with height $p=1$. Setting $z={\mathrm{tanh}}^{2}x$ and

$\psi \left(x\right)=\mathrm{exp}\left(\frac{\sqrt{{v}_{1}}}{2}{\mathrm{tanh}}^{2}x\right){\left(1-{\mathrm{tanh}}^{2}x\right)}^{\frac{\sqrt{-E}}{2}}\varphi \left(x\right)$

allows to write (23) as

${B}_{2}\left(z\right){\varphi }^{″}\left(z\right)+{B}_{1}\left(z\right){\varphi }^{\prime }\left(z\right)+{B}_{0}\left(z\right)\varphi \left(z\right)=0$ (24)

where

$\begin{array}{l}{B}_{2}\left(z\right)=-4z+4{z}^{2}\\ {B}_{1}\left(z\right)=-2+\left(6+4\sqrt{-E}-4\sqrt{{v}_{1}}\right)z+4\sqrt{{v}_{1}}{z}^{2}\\ {B}_{0}\left(z\right)=\left(\sqrt{-E}-E-\sqrt{{v}_{1}}-{v}_{1}-{v}_{2}-{v}_{3}\right)+\left({v}_{1}+{v}_{2}+3\sqrt{{v}_{1}}+2\sqrt{{v}_{1}}\sqrt{-E}\right)z\end{array}$

We can reduce the height of Equation (24) from $p=1$ to $p=0$ by taking its first derivative:

$\frac{\text{d}}{\text{d}z}\left[{B}_{2}\left(z\right){\varphi }^{″}\left(z\right)+{B}_{1}\left(z\right){\varphi }^{\prime }\left(z\right)+{B}_{0}\left(z\right)\varphi \left(z\right)\right]=0$

which implies that

${A}_{3}\left(z\right){\varphi }^{‴}\left(z\right)+{A}_{2}\left(z\right){\varphi }^{″}\left(z\right)+{A}_{1}\left(z\right){\varphi }^{\prime }\left(z\right)+{A}_{0}\left(z\right)\varphi \left(z\right)=0$ (25)

where

$\begin{array}{l}{A}_{3}\left(z\right)=-4z+4{z}^{2}\\ {A}_{2}\left(z\right)=-6+\left(14+4\left(\sqrt{-E}-\sqrt{{v}_{1}}\right)\right)z+4\sqrt{{v}_{1}}{z}^{2}\\ {A}_{1}\left(z\right)=6+\sqrt{-E}-E+4\left(\sqrt{-E}-\sqrt{{v}_{1}}\right)-\sqrt{{v}_{1}}-{v}_{1}-{v}_{2}-{v}_{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(11\sqrt{{v}_{1}}+2\sqrt{-E}\sqrt{{v}_{1}}+{v}_{1}+{v}_{2}\right)z\\ {A}_{0}\left(z\right)={{B}^{\prime }}_{0}\left(z\right)=3\sqrt{{v}_{1}}+2\sqrt{-E}\sqrt{{v}_{1}}+{v}_{1}+{v}_{2}\end{array}$ (26)

We solved Equation (25) by Algorithm ((12)-(13)-(14)) with $n=5$, ${v}_{1}=1$ and ${v}_{2}=-50$.

First we use Equations (10)-(11) to compute the canonical polynomials associated with Equation (25). Here are some of them:

$\begin{array}{l}{Q}_{0}\left(z\right)=-\frac{1}{20},\\ {Q}_{1}\left(z\right)=\frac{{v}_{3}}{640}-\frac{z}{32}-\frac{71}{160},\\ {Q}_{2}\left(z\right)=-\frac{{v}_{3}^{2}}{11520}+\left(\frac{{v}_{3}}{576}-\frac{173}{288}\right)z+\frac{7{v}_{3}}{128}-\frac{{z}^{2}}{36}-\frac{12259}{1440},\\ {Q}_{3}\left(z\right)=\frac{{v}_{3}^{3}}{122880}+\left(-\frac{{v}_{3}^{2}}{6144}+\frac{127{v}_{3}}{1024}-\frac{17947}{768}\right)z-\frac{523{v}_{3}^{2}}{61440}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(\frac{{v}_{3}}{384}-\frac{13}{12}\right){z}^{2}+\frac{22487{v}_{3}}{7680}-\frac{{z}^{3}}{32}-\frac{1271741}{3840},\end{array}$

$\begin{array}{c}{Q}_{4}\left(z\right)=-\frac{{v}_{3}^{4}}{614400}+\frac{77{v}_{3}^{3}}{30720}+\left(\frac{{v}_{3}^{3}}{30720}-\frac{157{v}_{3}^{2}}{3840}+\frac{129889{v}_{3}}{7680}-\frac{1474407}{640}\right)z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\left(-\frac{{v}_{3}^{2}}{1920}+\frac{91{v}_{3}}{192}-\frac{534}{5}\right){z}^{2}-\frac{219017{v}_{3}^{2}}{153600}+\left(\frac{{v}_{3}}{160}-\frac{247}{80}\right){z}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{681721{v}_{3}}{1920}-\frac{{z}^{4}}{20}-\frac{104477841}{3200}\end{array}$

From (12), we have

$n!\underset{i=1}{\overset{3}{\sum }}\frac{{a}_{ii}}{\left(n-i\right)!}=24\left(-1+\sqrt{-E}\right)=0$

which gives $E=-169$.

From (26), ${A}_{0}=-20$ and therefore ${B}_{0}\left(z\right)=230-{v}_{3}-20z$ where ${v}_{3}$ takes the six values:

$\begin{array}{l}229.2464376447439023,\text{\hspace{0.17em}}\text{\hspace{0.17em}}281.6583890371778959,\text{\hspace{0.17em}}\text{\hspace{0.17em}}344.0001059505978400,\\ 415.4312871412210056,\text{\hspace{0.17em}}\text{\hspace{0.17em}}495.5461344115843794,\text{\hspace{0.17em}}\text{\hspace{0.17em}}584.1176458146749767\end{array}$

which are the zeros of the following polynomial

$P\left({v}_{3}\right)=\frac{{v}_{3}^{6}-2350{v}_{3}^{5}+2256468{v}_{3}^{4}-1132669000{v}_{3}^{3}+313363258912{v}_{3}^{2}-45290385058560{v}_{3}+2670960608870400}{122880}$

whose the plot is given in Figure 1. Further, for $n=5$, we obtain obtain $\varphi \left(z\right)$ for the six values of ${v}_{3}$ are

Figure 1. Plot of $P\left({v}_{3}\right)$ whose the zeros are ${v}_{3}$. Here $n=5$.

$\begin{array}{c}{\varphi }_{5,1}\left(z\right)=848965.2873783716955968+319874.1407437577997468z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+44578.5802899894901557{z}^{2}+2890.7267231984532218{z}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+87.6883905888140244{z}^{4}+{z}^{5}\end{array}$

$\begin{array}{c}{\varphi }_{5,2}\left(z\right)=-3233.8652725518631273+83528.1351716517372919z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+21688.9752060886286782{z}^{2}+1968.4480715340372917{z}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+74.5854027407055260{z}^{4}+{z}^{5}\end{array}$

$\begin{array}{c}{\varphi }_{5,3}\left(z\right)=22.2163496075455889-1266.3331045478599408z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+6294.6494541225146804{z}^{2}+1094.9987219712643521{z}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+58.9999735123505400{z}^{4}+{z}^{5}\end{array}$

$\begin{array}{c}{\varphi }_{5,4}\left(z\right)=-0.1710184447933733+15.8560851714625286z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-173.3804428508567195{z}^{2}+392.8109983168830377{z}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+41.1421782146947486{z}^{4}+{z}^{5}\end{array}$

$\begin{array}{c}{\varphi }_{5,5}\left(z\right)=0.0013228922922234-0.1756444672214074z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+3.0942041854688765{z}^{2}-15.3305321180491017{z}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+21.1134663971039051{z}^{4}+{z}^{5}\end{array}$

$\begin{array}{c}{\varphi }_{5,6}\left(z\right)=-9.8748769085397727×{10}^{-6}+0.0017484340817809z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-0.0437115352456713{z}^{2}+0.3460170974111983{z}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-1.0294114536687442{z}^{4}+{z}^{5}\end{array}$

The graphs of the six functions are shown in Figure 2.

Example 2. The Schrodinger’s equation of $O\left(N\right)$ invariant decatic anharmonic oscillator in N-dimensional spherical coordinates is

$-\frac{{\text{d}}^{2}R\left(r\right)}{\text{d}{r}^{2}}-\frac{N-1}{r}\frac{\text{d}R\left(r\right)}{\text{d}r}+\left[\frac{l\left(l+N-2\right)}{{r}^{2}}+2\left(V\left(r\right)-E\right)\right]R\left(r\right)=0$ (27)

where R stands for the radial wave function. Setting $R\left(r\right)={r}^{\frac{1-N}{2}}\psi \left(r\right)$ transforms Equation (27) to

Figure 2. Plot of the wave function $\Psi \left(x\right)$ for the six values of ${v}_{3}$ that are the roots of $p\left({v}_{3}\right)=0$ that appear at the top of each plot. $n=5$.

$-\frac{{\text{d}}^{2}\psi \left(r\right)}{\text{d}{r}^{2}}+\left[\frac{\mu \left(\mu -1\right)}{{r}^{2}}+2\left(V\left(r\right)-E\right)\right]\psi \left(r\right)=0$ (28)

where $\mu =\mathcal{l}+\frac{1}{2}\left(N-1\right)$, $\mathcal{l}$ being a positive integer.

Further, consider the transformation

$\psi \left(z\right)={r}^{\mu }\mathrm{exp}\left(-\alpha \frac{{r}^{2}}{2}-\beta \frac{{r}^{4}}{4}-\gamma \frac{{r}^{6}}{6}\right)\varphi \left(z\right);\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}z={r}^{2}$

where $\gamma =\sqrt{2}$ and $\left\{\alpha ,\beta \right\}$ are parameters that depend on two unknowns ${\lambda }_{3}$ and ${\lambda }_{4}$ :

$\beta =\frac{{\lambda }_{4}}{\gamma }\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{and}\text{\hspace{0.17em}}\text{ }\text{\hspace{0.17em}}\alpha =\frac{1}{\gamma }\left({\lambda }_{3}-\frac{{\lambda }_{4}^{2}}{4}\right)$

${\lambda }_{3}$ and ${\lambda }_{4}$ should computed with the eigenvalue E. This yields a second order ODE of the form (1) with height $p=2$ :

${B}_{2}\left(z\right)\frac{{\text{d}}^{2}\varphi \left(z\right)}{\text{d}{z}^{2}}+{B}_{1}\left(z\right)\frac{\text{d}\varphi }{\text{d}z}+{B}_{0}\left(z\right)\varphi \left(z\right)=0$

where

$\begin{array}{l}{B}_{2}\left(z\right)=4z\\ {B}_{1}\left(z\right)=4\left(\left(N+2\mathcal{l}\right)/2-\gamma {z}^{3}-\beta {z}^{2}-\alpha z\right)\\ {B}_{0}\left(z\right)=\left(2\alpha \beta -\gamma \left(M+4\right)-2{\lambda }_{2}\right){z}^{2}+\left({\alpha }^{2}-\beta \left(M+2\right)-2{\lambda }_{1}\right)z+2E-\alpha M\end{array}$

with

$M=N+2\mathcal{l},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\gamma =\sqrt{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\beta ={\lambda }_{4}/\gamma ,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\alpha =\left({\lambda }_{3}-{\lambda }_{4}^{2}/4\right)/\gamma$

In order to obtain an equation with height $p=0$, we differentiate it twice to get:

${A}_{4}\left(z\right){\varphi }^{\left(4\right)}\left(z\right)+{A}_{3}\left(z\right){\varphi }^{‴}\left(z\right)+{A}_{2}\left(z\right){\varphi }^{″}\left(z\right)+{A}_{1}\left(z\right){\varphi }^{\prime }\left(z\right)+{A}_{0}\left(z\right)\varphi \left(z\right)=0$ (29)

where

$\begin{array}{l}{A}_{4}\left(z\right)={B}_{2}\left(z\right)=4z\\ {A}_{3}\left(z\right)=2{{B}^{\prime }}_{2}\left(z\right)+{B}_{1}\left(z\right)=8+2\left(2\mathcal{l}+N\right)-4\alpha z-4\beta {z}^{2}-4\gamma {z}^{3}\\ {A}_{2}\left(z\right)={{B}^{″}}_{2}2\left(z\right)+2{{B}^{\prime }}_{1}\left(z\right)+{B}_{0}\left(z\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=-8\alpha +2E-\alpha M+\left({\alpha }^{2}-16\beta -2{\lambda }_{1}-\beta \left(2+M\right)\right)z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }+\left(2\alpha \beta -24\gamma -2{\lambda }_{2}-\gamma \left(4+M\right)\right){z}^{2}\end{array}$

$\begin{array}{l}{A}_{1}\left(z\right)={{B}^{″}}_{1}\left(z\right)+2{{B}^{\prime }}_{0}\left(z\right)=2{\alpha }^{2}-8\beta -4{\lambda }_{1}-2\beta \left(2+M\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{ }+\left(-24\gamma +4\left(2\alpha \beta -2{\lambda }_{2}-\gamma \left(4+M\right)\right)\right)z\\ {A}_{0}\left(z\right)={{B}^{″}}_{0}\left(z\right)=2\left(2\alpha \beta -2{\lambda }_{2}-\gamma \left(4+M\right)\right)\end{array}$

The canonical polynomials are obtained by recursion (10)-(11):

$\begin{array}{l}{Q}_{0}\left(z\right)=\frac{1}{4\alpha \beta -18\gamma -4}\\ {Q}_{1}\left(z\right)=\frac{2-{\alpha }^{2}+11\beta +\left(-2+2\alpha \beta -9\gamma \right)z}{6\left(-2+2\alpha \beta -13\gamma \right)\left(-2+2\alpha \beta -9\gamma \right)}\\ {Q}_{2}\left(z\right)=\left(4-26\alpha -4{\alpha }^{2}+{\alpha }^{4}+52\beta +165{\beta }^{2}+4E-4\alpha \beta E-169\alpha \gamma +26E\gamma \right)\frac{1}{{D}_{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{ }+\left(-4+2{\alpha }^{2}-30\beta +4\alpha \beta -2{\alpha }^{3}\beta +30al{\beta }^{2}-18\gamma +9{\alpha }^{2}\gamma -135\beta \gamma \right)\frac{z}{{D}_{2}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{ }+\left(4-8\alpha \beta +4{\alpha }^{2}{\beta }^{2}+44\gamma -44\alpha \beta \gamma +117{\gamma }^{2}\right)\frac{{z}^{2}}{{D}_{2}}\end{array}$ (30)

where ${D}_{2}:=12\left(-2+2\alpha \beta -17\gamma \right)\left(-2+2\alpha \beta -13\gamma \right)\left(-2+2\alpha \beta -9\gamma \right)$

We have applied Algorithm ((12)-(13)-(14)) for different sets of parameters:

1) For $n=2$, ${\lambda }_{1}=1$, ${\lambda }_{2}=1$, $N=3$, $\mathcal{l}=1$, $\gamma =\sqrt{2}$, $\beta ={\lambda }_{4}/\gamma$, $\alpha =\left({\lambda }_{3}-{\lambda }_{4}^{2}/4\right)/\gamma$, $p=0$ ; $\nu =4$, the unknown $E,{\lambda }_{3},{\lambda }_{4}$ are determined by solving the following system:

$-24-204\sqrt{2}+12{\lambda }_{3}{\lambda }_{4}-3{\lambda }_{4}^{3}=0$ (31)

$\begin{array}{l}-15184570-4233890\sqrt{2}+\left(347972E+79152\sqrt{2}E\right){\lambda }_{4}\\ +\left(200600+34760\sqrt{2}\right){\lambda }_{4}^{2}\\ +\left(179700+454070\sqrt{2}-4656E-544\sqrt{2}E\right){\lambda }_{4}^{3}\\ +\left(-680-40\sqrt{2}-2960E-10236\sqrt{2}E\right){\lambda }_{4}^{4}\\ +\left(-520-4100\sqrt{2}+16E-32\sqrt{2}{E}^{2}\right){\lambda }_{4}^{5}\\ +\left(-11625-825\sqrt{2}+104\sqrt{2}E\right){\lambda }_{4}^{6}+330E{\lambda }_{4}^{7}=0\end{array}$ (32)

$\begin{array}{l}106642188+34864280\sqrt{2}-2087832{\lambda }_{4}^{2}-474912\sqrt{2}{\lambda }_{4}^{2}-1749968{\lambda }_{4}^{3}\\ \text{ }-3254098\sqrt{2}{\lambda }_{4}^{3}+13968{\lambda }_{4}^{4}+1632\sqrt{2}{\lambda }_{4}^{4}+15968{\lambda }_{4}^{5}+59960\sqrt{2}{\lambda }_{4}^{5}\\ \text{ }+110898{\lambda }_{4}^{6}+17820\sqrt{2}{\lambda }_{4}^{6}-264\sqrt{2}{\lambda }_{4}^{7}-1388{\lambda }_{4}^{8}-1155\sqrt{2}{\lambda }_{4}^{9}\\ +\left(-6528{\lambda }_{4}^{4}-27936\sqrt{2}{\lambda }_{4}^{4}+192\sqrt{2}{\lambda }_{4}^{6}+1184{\lambda }_{4}^{7}\right)E=0\end{array}$ (33)

Note that Equation (31) and Equation (33) are linear in ${\lambda }_{3}$ and E respectively. So we compute ${\lambda }_{3}$ and E in terms of ${\lambda }_{4}$ and substitute their expressions in Equation (32) which gives the values of ${\lambda }_{4}$. As a result we get for $E>0$ two sets of solutions:

$\begin{array}{c}\underset{_}{\text{Set1}}\\ {\lambda }_{3}=10.4794241513810915978439288620\\ {\lambda }_{4}=3.57541756044743705617752580655\\ E=10.3921270096757567640709427917\end{array}$

$\begin{array}{c}\underset{_}{\text{Set2}}\\ \begin{array}{l}{\lambda }_{3}=10.4794241513810915978439288620\\ {\lambda }_{4}=2.96736011311581477714971460010\\ E=45.4729663748453390359477629584\end{array}\end{array}$

For Set 1 we have

$\begin{array}{l}\alpha =6.20558774817446509454508658219,\\ \beta =2.09824045820667335037013178954,\\ \gamma =\sqrt{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu =2\end{array}$

and the exact solution of Equation (29) when $n=2$ is

$\begin{array}{c}{\varphi }_{2}\left(z\right)=0.148569159577064357603068196653\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-0.890196601343371995098741638225z+{z}^{2}\end{array}$

Then the wave function for $n=2$ is

${\psi }_{2}\left(z\right)={r}^{\mu }\mathrm{exp}\left(-\alpha \frac{{r}^{2}}{2}-\beta \frac{{r}^{4}}{4}-\gamma \frac{{r}^{6}}{6}\right){\varphi }_{2}\left(z\right);\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}z={r}^{2}$

which is plotted in Figure 3.

2) For $n=4$ ; ${\lambda }_{1}=1$ ; ${\lambda }_{2}=1$ ; $N=6$ ; $\mathcal{l}=1$ ; $\gamma =\sqrt{2}$ ; $p=0$ ; $\nu =4$ ; $M=N+2l$ we obtain a system of equations with unknown $E,{\lambda }_{3},{\lambda }_{4}$. This system has three sets of solutions:

$\begin{array}{c}\text{Set1}\\ \begin{array}{l}{\lambda }_{3}=15.2213022374713691334593627854\\ {\lambda }_{4}=3.35080280384336327546445916262\\ E=122.054198917713672699520321911\end{array}\end{array}$

$\begin{array}{c}\text{Set2}\\ \begin{array}{l}{\lambda }_{3}=14.6173922152479058765833382408\\ {\lambda }_{4}=3.74178054365072115240609371379\\ E=64.518363537276464736709058162\end{array}\end{array}$

Figure 3. Plot of wave function ${\psi }_{2}\left(x\right)$ for the Set 1 ( $n=2$ ).

$\begin{array}{c}\begin{array}{c}\text{Set3}\\ {\lambda }_{3}=14.3050849085847171927078631515\\ {\lambda }_{4}=4.24082480839216204741958029044\\ E=20.6723051499346259885190817057\end{array}\end{array}$

The exact polynomial solution of Equation (29) that corresponds to each set of the computed parameters above are:

$\begin{array}{l}\text{Set}1:\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\varphi }_{4}\left(z\right)=0.05123329350714747-0.55678527469919100z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+1.865828731884573800{z}^{2}-2.378365009178223910{z}^{3}+{z}^{4}\end{array}$

$\begin{array}{l}\text{Set}2:\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\varphi }_{4}\left(z\right)=1.83303106908165349-7.57827320512557611z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+5.04729517062652433{z}^{2}+1.59031054924421090{z}^{3}+{z}^{4}\end{array}$

$\begin{array}{l}\text{Set}3:\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\varphi }_{4}\left(z\right)=40.08040813732159932+35.42885586322148665z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+20.94732836908738593{z}^{2}+5.631928970610121495{z}^{3}+{z}^{4}\end{array}$

Figure 4 shows ${\psi }_{4,1}\left(x\right)$.

3) For $n=8$ ; ${\lambda }_{1}=1$ ; ${\lambda }_{2}=1$ ; $N=6$ ; $\mathcal{l}=1$. Here are the results:

$\begin{array}{c}\text{Set1}\\ {\lambda }_{3}=20.56410237324861185\\ {\lambda }_{4}=3.779556076814648561\\ E=284.2200173241970169\end{array}$ $\begin{array}{c}\text{Set2}\\ {\lambda }_{3}=20.56410237324861185\\ {\lambda }_{4}=3.779556076814648561\\ E=284.220017324197016903\end{array}$

$\begin{array}{c}\text{Set3}\\ {\lambda }_{3}=19.97053471979101493\\ {\lambda }_{4}=4.043790238076919932\\ E=200.814293836308274681\end{array}$ $\begin{array}{c}\text{Set4}\\ {\lambda }_{3}=19.49903355415370773\\ {\lambda }_{4}=4.34676919014164935\\ E=128.39046011355799670\end{array}$

$\begin{array}{c}\text{Set5}\\ {\lambda }_{3}=19.19007126134549327\\ {\lambda }_{4}=4.69538124652545845\\ E=68.10813774437088605\end{array}$ $\begin{array}{c}\text{Set6}\\ {\lambda }_{3}=19.09600441735804689\\ {\lambda }_{4}=5.10347417185941959\\ E=21.86369137050343994\end{array}$

The exact solution that corresponds to Set 1:

Figure 4. Plot of wave function ${\psi }_{4,1}\left(x\right)$ for the Set 1 ( $n=4$ ).

Figure 5. Plot of wave function ${\psi }_{8,1}\left(x\right)$ for the Set 1 ( $n=8$ ).

$\begin{array}{c}{\varphi }_{8,1}\left(z\right)=0.0018666121867654-0.0551016926965569z\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+0.5790267235874925{z}^{2}-2.9461631765945323{z}^{3}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+8.1634406700523406{z}^{4}-12.8729691673963369{z}^{5}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+11.4573865169591070{z}^{6}-5.3265009701024445{z}^{7}+{z}^{8},\end{array}$

Figure 5 shows ${\psi }_{8,1}\left(x\right)$.

5. Conclusion

In this paper we have extended the canonical polynomials approach that was developed in  to solve QES differential equations of arbitrary high order $\nu \ge 2$. While the existing methods for solving QESs require the solution of a nonlinear algebraic system whose dimensions depend on the desired degree of $y\left(x\right)$, our new approach requires solving a nonlinear algebraic system of dimensions depending on p, the height of the differential operator. This advantage is due to the fact that the sequence of canonical polynomials enjoys the permanence characteristic.

Acknowledgements

The financial support of the Public Authority for Applied Education and Training (PAAET), Kuwait, during this research is greatly appreciated.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

  Schafke, R. and Schmidt, D. (1980) The Connection Problem for General Linear Ordinary Differential Equations at Two Regular Singular Points with Applications in the Theory of Special Functions. SIAM Journal on Mathematical Analysis, 11, 848-862. https://doi.org/10.1137/0511076  Bender, C.M. and Dunne, G.V. (1996) Quasi-Exactly Solvable Systems and Orthogonal Polynomials. Journal of Mathematical Physics, 37, 6-11. https://doi.org/10.1063/1.531373  Sasaki, R., Yang, W.L. and Zhang, Y.Z. (2009) Bethe Ansatz Solutions to Quasi-Exactly Solvable Difference Equations. Symmetry, Integrability and Geometry: Methods and Applications, 5, 104. https://doi.org/10.3842/SIGMA.2009.104  Gomez-Ullate, D., Kamran, N. and Milson, R. (2009) An Extended Class of Orthogonal Polynomials Defined by a Sturm-Liouville Problem. Journal of Mathematical Analysis and Applications, 359, 352-367. https://doi.org/10.1016/j.jmaa.2009.05.052  Pan, F., Klauder, J.R. and Draayer, J.P. (1999) Quasi-Exactly Solvable Cases of an N-Dimensional Symmetric Decatic Anharmonic Oscillator. Physics Letters A, 262, 131-136. https://doi.org/10.1016/S0375-9601(99)00651-9  Moroz, A. and Miroshnichenko, A.E. (2018) Constraint Polynomial Approach—An Alternative to the Functional Bethe Ansatz Method? arXiv:1807.11871v1.  Xie, Q.-T. (2012) New Quasi-Exactly Solvable Double-Well Potentials. Journal of Physics A: Mathematical and Theoretical, 45, Article ID: 175302. https://doi.org/10.1088/1751-8113/45/17/175302  El-Daou, M.K., AlZanki, T.H. and Al-Mutawa, N.S. (2019) Quasi-Exactly Solvable Differential Models: A Canonical Polynomials Approach. American Journal of Computational Mathematics, 9, 48-60. https://doi.org/10.4236/ajcm.2019.92004  El-Daou, M.K. (2012) Duality between the Generalized Canonical Polynomials and Some Special Determinants. International Journal of Applied Mathematics, 25, 811-824.  Lanczos, C. (1956) Applied Analysis. Prentice-Hall, Englewood Cliffs, NJ.  Ortiz, E.L. (1969) The Tau Method. SIAM Journal on Numerical Analysis, 6, 480-492. https://doi.org/10.1137/0706044  Froes Bunchaft, M.E. (1997) Some Extensions of the Lanczos-Ortiz Theory of Canonical Polynomials in the Tau Method. Mathematics of Computation, 66, 609-621. https://doi.org/10.1090/S0025-5718-97-00816-8  El-Daou, M.K. and Ortiz, E.L. (1994) A Recursive Formulation of Collocation in Terms of Canonical Polynomials. Computing, 52, 177-202. https://doi.org/10.1007/BF02238075  El-Daou, M.K. and Ortiz, E.L. (1998) A Recursive Formulation of Galerkin Method in Terms of the Tau Method: Bounded and Unbounded Domains. Computers & Mathematics with Applications, 35, 83-94. https://doi.org/10.1016/S0898-1221(98)00098-4  El-Daou, M.K. and Al-Hamad, K.M. (2012) Computation of the Canonical Polynomials and Applications to Some Optimal Control Problems. Numerical Algorithms, 61, 545-566. https://doi.org/10.1007/s11075-012-9550-5  Alexandros, G. (2016) Paraskevopoulos. Complete Bases of Canonical Polynomials and the Full Recursive Formulation of the Spectral-Tau Method. https://arxiv.org/abs/1611.09135  Namasivayam, S. and Ortiz, E.L. (1981) Best Approximation and the Numerical Solution of Partial Differential Equations with the Tau Method. Portugaliae mathematica, 40, 97-119.     customer@scirp.org +86 18163351462(WhatsApp) 1655362766  Paper Publishing WeChat 