4.1.5 Example 5 case one

Solve Bessel ode (from Kovacic original paper)

\begin{align} y^{\prime \prime }-\frac {4\left ( m^{2}-x^{2}\right ) -1}{4x^{2}}y & =0\tag {1}\\ y^{\prime \prime }+ay^{\prime }+by & =0\nonumber \end{align}

Hence

\begin{align*} a & =0\\ b & =\frac {4\left ( m^{2}-x^{2}\right ) -1}{4x^{2}}\end{align*}

It is first transformed to the following ode by eliminating the first derivative

\begin{equation} z^{\prime \prime }=rz \tag {2}\end{equation}

Using what is known as the Liouville transformation given by

\begin{equation} y=ze^{\frac {-1}{2}\int adx} \tag {3}\end{equation}

Where it can be found that \(r\) in (2) is given by

\begin{align} r & =\frac {1}{4}a^{2}+\frac {1}{2}a^{\prime }-b\nonumber \\ & =\frac {4\left ( m^{2}-x^{2}\right ) -1}{4x^{2}} \tag {4}\end{align}

Hence the DE we will solve using Kovacic algorithm is Eq (2) which is

\begin{equation} z^{\prime \prime }=\frac {4\left ( m^{2}-x^{2}\right ) -1}{4x^{2}}z \tag {5}\end{equation}

Step 0 We need to find which case it is. \(r=\frac {s}{t}\) where

\begin{align*} s & =4\left ( m^{2}-x^{2}\right ) -1\\ t & =4x^{2}\end{align*}

The free square factorization of \(t\) is \(t=\left [ 1,x\right ] \). Hence

\begin{equation} m=2 \tag {6}\end{equation}

Since \(m\) is number of elements in the free square factorization. in this case we set

\begin{align*} t_{1} & =1\\ t_{2} & =x \end{align*}

Now

\begin{align*} O\left ( \infty \right ) & =\deg \left ( t\right ) -\deg \left ( s\right ) \\ & =2-2\\ & =0 \end{align*}

There is one pole at \(x=0\) of order 2. Looking at the cases table giving up, reproduced here

case	allowed pole order for \(r=\frac {s}{t}\)	allowed \(O\left ( \infty \right ) \) order	\(L\)
1	\(\left \{ 0,1,2,4,6,8,\cdots \right \} \)	\(\left \{ \cdots ,-8,-6,-4,-2,0,2,3,4,5,6,7,\cdots \right \} \)	\(\left [ 1\right ] \)
2	\(\left \{ 2,3,5,7,9,\cdots \right \} \)	no condition	\(\left [ 2\right ] \)
3	\(\left \{ 1,2\right \} \)	\(\left \{ 2,3,4,5,6,7,\cdots \right \} \)	\(\left [ 4,6,12\right ] \)

Shows that only case 1,2 are possible. \(L=\left [ 1,2\right ] \).

Step 1

This step has 4 parts (a,b,c,d).

part (a) Here the fixed parts \(e_{fixed},\theta _{fixed}\) are calculated using

\begin{align*} e_{fixed} & =\frac {1}{4}\left ( \min \left ( O\left ( \infty \right ) ,2\right ) -\deg \left ( t\right ) -3\deg \left ( t_{1}\right ) \right ) \\ \theta _{fixed} & =\frac {1}{4}\left ( \frac {t^{\prime }}{t}+3\frac {t_{1}^{\prime }}{t_{1}}\right ) \end{align*}

Using \(O\left ( \infty \right ) =0,t=4x^{2},t_{1}=1\) the above gives

\begin{align*} e_{fixed} & =\frac {1}{4}\left ( \min \left ( 0,2\right ) -2-3\left ( 0\right ) \right ) \\ & =\frac {1}{4}\left ( 0-2\right ) \\ & =-\frac {1}{2}\\ \theta _{fixed} & =\frac {1}{4}\left ( \frac {\frac {d}{dx}\left ( 4x^{2}\right ) }{4x^{2}}+3\left ( 0\right ) \right ) \\ & =\frac {1}{2x}\end{align*}

part (b)

Here the values \(e_{i},\theta _{i}\) are found for \(i=1\cdots k_{2}\) where \(k_{2}\) is the number of roots of \(t_{2}\). In other words, the number of poles of \(r\) that are of order \(2\).

\[ r=\frac {4\left ( m^{2}-x^{2}\right ) -1}{4x^{2}}\]

These will be the zeros of \(t_{2}\) in the above square free factorization of \(t\). From above we found that

\[ t_{2}=x \]

Label these zeros of \(t_{2}\) as \(c_{1},c_{2},\cdots ,c_{k_{2}}\). The zeros of \(t_{2}\) are \(\left \{ 0\right \} \). Therefore \(k_{2}=1\). Hence

\[ M=1 \]

Now we iterate over each zero \(c_{i}\) times finding \(e_{i}\) and \(\theta _{i}\) from each. These are found to be (following formula in paper) to be

\begin{align*} b_{1} & =m^{2}-\frac {1}{4}\\ e_{1} & =\sqrt {1+4b}=\sqrt {1+4\left ( m^{2}-\frac {1}{4}\right ) }=2m\qquad m>0 \end{align*}

Where \(b_{1}\) is the coefficient of \(\frac {1}{\left ( x-c_{1}\right ) ^{2}}\) in the partial fractions decomposition of \(r\) which is \(r=-1+\left ( m^{2}-\frac {1}{4}\right ) \frac {1}{x^{2}}\).  And

\[ \theta _{1}=\frac {e_{1}}{x-c_{1}}=\frac {2m}{x}\]

Part (c)

This part applied only to case 1. It is used to generate \(e_{i},\theta _{i}\) for poles of \(r\) order \(4,6,8,\cdots ,M\) if any exist. Since non exist here. This is skipped. Hence \(M\) stays \(1\).

Part(d)

Now we need to find \(e_{0},\theta _{0}\). If \(O\left ( \infty \right ) >2\) then \(e_{0}=1,\theta _{0}=0\). But if \(O\left ( \infty \right ) =2\) then \(\theta _{0}=0\) and \(e_{0}=\sqrt {1+4b}\) where \(b\) is the coefficient of \(\frac {1}{x^{2}}\) in the Laurent series expansion of \(r\) at \(\infty \). Since \(O\left ( \infty \right ) =0\) here then none of these cases applies. For case 1 \(\left ( n=1\right ) \) following the method in the paper we find

\begin{align*} e_{0} & =0\\ \theta _{0} & =2i \end{align*}

Hence now we have

\begin{align*} e & =\left \{ 0,2m\right \} \\ \theta & =\left \{ 2i,\frac {2m}{x}\right \} \end{align*}

The above are arranged such that \(e_{0}\) is the first entry. Same for \(\theta \). This to keep the same notation as in the paper. The above complete step 1, which is to generates the candidate \(e^{\prime }s\) and \(\theta ^{\prime }s\). In step 2, these are used to generate trials \(d\) and \(\theta \) and find from them \(P\left ( x\right ) \) polynomial if possible.

Step 2

In this step, we now have all the \(e_{i},\theta _{i}\) values found above in addition to \(e_{fix},\theta _{fix}\).

Starting with \(n=1\). And since we have \(M=1\) then there are \(\left ( n+1\right ) ^{M+1}=2^{2}=4\) sets \(s\) to try. The first set \(s\) is

\[ s=\left \{ \frac {-n}{2},\frac {-n}{2}\right \} =\left \{ \frac {-1}{2},\frac {-1}{2}\right \} \]

Now we generate trial \(d\) using

\begin{equation} d=\left ( n\right ) \left ( e_{fix}\right ) +s_{0}e_{0}-\sum _{i=1}^{M}s_{i}e_{i}\nonumber \end{equation}

Since \(M=1\) then the above becomes

\begin{equation} d=\left ( n\right ) \left ( e_{fix}\right ) +s_{0}e_{0}-s_{1}e_{1} \tag {7}\end{equation}

If \(d\geq 0\) then we go and find a trial \(\Theta \). We need to have both \(d,\Theta \) to go to the next step.  \(\Theta \) is found using

\begin{equation} \Theta =\left ( n\right ) \left ( \theta _{fix}\right ) +\sum _{i=0}^{M}s_{i}\theta _{i} \tag {8}\end{equation}

Hence the first trial \(d\) is (using Eq (7)) and recalling that \(e_{fix}=-\frac {1}{2},\theta _{fixed}=\frac {1}{2x}\) gives

\begin{align*} d & =\left ( 1\right ) \left ( -\frac {1}{2}\right ) +\left ( \frac {-1}{2}\right ) \left ( 0\right ) -\left ( \frac {-1}{2}\right ) \left ( 2m\right ) \\ & =m-\frac {1}{2}\end{align*}

We now have to assume something about \(m\) to be to continue otherwise we will not be able to decide if \(d\) is integer and \(d\geq 0\). If we assume that \(m\) is half of all positive odd integers (\(\frac {1}{2},\frac {3}{2},\frac {5}{2},\cdots \)) then \(d\geq 0\). We can also assume that \(m\) is half of all negative odd integers and the other \(s\) set will match. So under the assumption that \(m\) is half of all positive odd integers the above \(d\) can be used for the next step.  To continue, we assume \(m\) takes some specific value to simplify the steps. Let \(m=\frac {3}{2}\) from now on. Hence \(d=1\). Therefore \(e,\theta \) become

\begin{align*} e & =\left \{ 0,3\right \} \\ \theta & =\left \{ 2i,\frac {3}{x}\right \} \end{align*}

Using Eq (8) gives (using \(M=1\))

\begin{align*} \Theta & =\left ( n\right ) \left ( \theta _{fix}\right ) +s_{0}\theta _{0}+s_{1}\theta _{1}\\ & =\left ( 1\right ) \left ( \frac {1}{2x}\right ) +\left ( \frac {-1}{2}\right ) \left ( 2i\right ) +\left ( \frac {-1}{2}\right ) \left ( \frac {3}{x}\right ) \\ & =-\frac {1}{x}\left ( ix+1\right ) \end{align*}

Now that we have good trial \(d\) and \(\Theta \), then step 3 is called to generate \(\omega \) if possible.

Step 3

The input to this step is the integer \(d=1\) and \(\Theta =-\frac {1}{x}\left ( ix+1\right ) \) found from step 2 and also \(r=\frac {4\left ( n^{2}-x^{2}\right ) -1}{4x^{2}}\) which comes from \(z^{\prime \prime }=rz\). This step is broken into these parts. First we find the \(p_{-1}\left ( x\right ) \) polynomial. If we are to solve for its coefficients, then next we build the minimal polynomial from the \(p_{i}\left ( x\right ) \) polynomials constructed during finding \(p_{-1}\left ( x\right ) \). The minimal polynomial \(p_{\min }\left ( x\right ) \) will be a function of \(\omega \). Next we solve for \(\omega \) from \(p_{\min }\left ( x\right ) =0\). If this is successful, then we have found \(\omega \) and the first solution to the ode \(y^{\prime \prime }=ry\) is \(e^{\int \omega dx}\,\). Below shows how this is done.

We start by forming a polynomial

\begin{align*} p\left ( x\right ) & =x^{d}+a_{d-1}x^{d-1}+\cdots +a_{0}\\ & =x+a_{0}\end{align*}

The goal is to solve for the \(a_{0}\) coefficient. Now depending on case  number \(n\), we do the following. Since we are in case \(n=1\) then

\begin{align*} p_{1} & =-p\\ & =-x-a_{0}\\ p_{0} & =-p_{1}^{\prime }-\Theta p_{1}\\ & =-\left ( -x-a_{0}\right ) ^{\prime }-\left ( -\frac {1}{x}\left ( ix+1\right ) \right ) \left ( -x-a_{0}\right ) \\ & =1-\left ( -\frac {1}{x}\left ( ix+1\right ) \right ) \left ( -x-a_{0}\right ) \\ & =-\frac {1}{x}\left ( ix^{2}+ia_{0}x+a_{0}\right ) \\ p_{-1} & =-p_{0}^{\prime }-\Theta p_{0}-\left ( 1\right ) \left ( 1\right ) rp_{1}\\ & =-\left ( -\frac {1}{x}\left ( ix^{2}+ia_{0}x+a_{0}\right ) \right ) ^{\prime }-\left ( -\frac {1}{x}\left ( ix+1\right ) \right ) \left ( -\frac {1}{x}\left ( ix^{2}+ia_{0}x+a_{0}\right ) \right ) -\frac {4\left ( n^{2}-x^{2}\right ) -1}{4x^{2}}\left ( -x-a_{0}\right ) \\ & =-\frac {2}{x}\left ( ia_{0}-1\right ) \end{align*}

Now we try to solve for \(a_{i}\) using \(p_{-1}\left ( x\right ) =0\). This gives \(a_{0}=-i\).  Hence this implies

\[ p\left ( x\right ) =x-i \]

Since this is case \(n=1\) then

\begin{align*} \omega & =\frac {p^{\prime }}{p}+\Theta \\ & =\frac {\left ( x-i\right ) ^{\prime }}{x-i}-\frac {1}{x}\left ( ix+1\right ) \\ & =\frac {-i\left ( ix-x^{2}+1\right ) }{\left ( -x+i\right ) x}\\ & =-\frac {\left ( ix^{3}+1\right ) }{x^{3}+x}\end{align*}

Before using this, we will verify it is correct. For case 1 the above should satisfy

\[ \omega ^{\prime }+\omega ^{2}=r \]

Let us see if this is the case or not.

\begin{align*} \frac {d}{dx}\left ( \frac {-i\left ( ix-x^{2}+1\right ) }{\left ( -x+i\right ) x}\right ) +\left ( \frac {-i\left ( ix-x^{2}+1\right ) }{\left ( -x+i\right ) x}\right ) ^{2} & =\frac {4\left ( \left ( \frac {3}{2}\right ) ^{2}-x^{2}\right ) -1}{4x^{2}}\\ \frac {\left ( -2ix^{3}+3x^{2}+1\right ) }{x^{2}\left ( x^{2}+1\right ) ^{2}}+\frac {\left ( ix^{2}+x-i\right ) ^{2}}{x^{2}\left ( x-i\right ) ^{2}} & =-\frac {1}{x^{2}}\left ( x^{2}-2\right ) \\ -\frac {1}{x^{2}}\left ( x^{2}-2\right ) & =-\frac {1}{x^{2}}\left ( x^{2}-2\right ) \end{align*}

Verified. Hence \(\omega =\frac {-i\left ( ix-x^{2}+1\right ) }{\left ( -x+i\right ) x}=-\frac {\left ( ix^{3}+1\right ) }{x^{3}+x}\) will give the solution to the ode \(y^{\prime \prime }-\frac {4\left ( m^{2}-x^{2}\right ) -1}{4x^{2}}y=0\) when \(m=\frac {3}{2}\). Since solution \(\omega \) is found and verified, then first solution to the ode is

\begin{align*} z & =e^{\int \omega dx}\\ & =e^{\int -\frac {\left ( ix^{3}+1\right ) }{x^{3}+x}dx}\\ & =\frac {1}{x}e^{-ix}\left ( x-i\right ) \end{align*}

Hence first solution to given ODE is

\begin{align*} y & =ze^{\frac {-1}{2}\int adx}\\ & =\frac {1}{x}e^{-ix}\left ( x-i\right ) e^{-\frac {1}{2}\int 0dx}\\ & =\frac {1}{x}e^{-ix}\left ( x-i\right ) \end{align*}

One difficulty in implementation of Kovacic algorithm using an ode with a parameter \(m\) like in this Bessel ode example, is that it makes it hard to decide if \(d\geq 0\) or not. So in practice, it is better to use this algorithm for specific values of any parameters that can be involved.