Series Solution near a Regular Singular Point: Roots Differing by an Integer

Here we outline two ways of obtaining the second solution from the first solution for a differential equation that has a regular singlar point and roots of the indicial equation that differ by an integer. This will be done by considering the example differential equation

x y'' + 3 y' - y = 0, x>0.

Getting the first solution

First, we set the problem up and obtain the series solution for the first case. In any method we need to determine one solution first! We have

p(x) = (3/x)   and    q(x) = ((-1)/x)

which means that x=0 is a singular point. The functions x p(x) = 3 and (x^2) q(x) = - x are both analytic at x=0, so x=0 is a regular singular point.

Assume a solution looks like y= (∑_(n=0)^∞)(Subscript[a, n])(x^(n+r)), and substitute into the DE:

x y'' + 3 y' - y = 0 
x (∑_(n=0)^∞)(n+r)(n+r-1)(Subscript[a, n])(x^(n+r-2))+ 3  (∑_(n=0)^∞)(n+r)(Subscript[a, n])(x^(n+r-1)) - (∑_(n=0)^∞)(Subscript[a, n])(x^(n+r)) = 0 
(∑_(n=0)^∞)(n+r)(n+r-1)(Subscript[a, n])(x^(n+r-1))+(∑_(n=0)^∞)(3n+3r)(Subscript[a, n])(x^(n+r-1)) -(∑_(n=0)^∞)(Subscript[a, n])(x^(n+r)) = 0 
((∑_(n=0)^∞)(n+r)(n+r-1)(Subscript[a, n])(x^(n-1))+(∑_(n=0)^∞)(3n+3r)(Subscript[a, n])(x^(n-1)) -(∑_(n=0)^∞)(Subscript[a, n])(x^n) )(x^r)= 0 
(∑_(n=0)^∞)(n+r)(n+r-1)(Subscript[a, n])(x^(n-1))+(∑_(n=0)^∞)(3n+3r)(Subscript[a, n])(x^(n-1)) -(∑_(n=0)^∞)(Subscript[a, n])(x^n) = 0 , since x>0
(∑_(n=0)^∞)(n+r)(n+r-1)(Subscript[a, n])(x^(n-1))+(∑_(n=0)^∞)(3n+3r)(Subscript[a, n])(x^(n-1)) - (∑_(n=1)^∞)(Subscript[a, n-1])(x^(n-1)) = 0 
(r(r-1)(Subscript[a, 0])+3r (Subscript[a, 0]))(x^(-1))+(∑_(n=1)^∞)(n+r)(n+r-1)(Subscript[a, n])(x^(n-1))+(∑_(n=1)^∞)(3n+3r)(Subscript[a, n])(x^(n-1)) - (∑_(n=1)^∞)(Subscript[a, n-1])(x^(n-1)) = 0 
(r(r-1)(Subscript[a, 0])+3r (Subscript[a, 0]))(x^(-1))+(∑_(n=1)^∞)((n+r)(n+r-1)(Subscript[a, n])+(3n+3r)(Subscript[a, n])-(Subscript[a, n-1]))(x^(n-1))= 0

Now, this has top be true for all coefficients of powers of  x, so we require

r(r-1)(Subscript[a, 0])+3r (Subscript[a, 0])=0
(n+r)(n+r-1)(Subscript[a, n])+(3n+3r)(Subscript[a, n])-(Subscript[a, n-1])= 0, n=1,2,3,4,...

So we can choose the indicial equation to be (and also require (Subscript[a, 0])≠0):

r(r-1)+3r =0 (indicial equation)
r=0,r=-2

These two roots differ by an integer.

Largest Root of Indicial Equation

We pick (Subscript[r, 1])=0 since it is the largest root. The recurence relation is therefore:

n(n-1)(Subscript[a, n])+3n (Subscript[a, n])-(Subscript[a, n-1])= 0, n=1,2,3,4,... 
(Subscript[a, n])= ((Subscript[a, n-1])/(n(n+2))), n=1,2,3,4,...

Here is the recurrence relation coded into Mathematica (remember, (Subscript[a, 0]) is arbitrary since the recursion relation doesn't determine it):

a[0]=1;
a[n_]:=a[n]=a[n-1]/(n(n+2));

Now, we can get the first few terms of the solution:

n=5;
y1 = Sum[a[i]x^i,{i,0,n}]

Let's look for the pattern. If we can find it great--if not, we can still use what we have above as an approximation to the solution. Let's see if Mathematica can do something useful in finding the pattern (until now we have looked for patterns solely by hand):

Clear[a]
a[0]=1;
a[n_]:=a[n]=a[n-1]*HoldForm[1/(n(n+2))];

n=5
Table[{n,a[n]},{n,0,4}]//TableForm

From this, we can ascertain the pattern as:

(Subscript[a, i])= ((2!)/(i!(i+2)!)), i=1,2,3,4,...

Let's check if we got it right:

n=5;
y1 = Sum[((2!)/(i!(i+2)!))x^i,{i,0,n}]

We absolutely rock! Let's see if Mathematica knows this function:

y1 = Sum[((2!)/(i!(i+2)!))x^i,{i,0,Infinity}]

And, a final check:

FullSimplify[x*D[y1,{x,2}]+3D[y1,x]-y1==0]

The Other Root of Indicial Equation

What if we tried to use the other indicial root? What would happen?

For (Subscript[r, 2])=-2 the recurence relation is therefore:

(n-2)(n-3)(Subscript[a, n])+(3n-6)(Subscript[a, n])-(Subscript[a, n-1])= 0, n=1,2,3,4,... 
(Subscript[a, n])= ((Subscript[a, n-1])/((n-2)(n-3)+(3n-6))), n=1,2,3,4,... 
= ((Subscript[a, n-1])/((n-2)(n-3)+3(n-2)))
= ((Subscript[a, n-1])/(n(n-2))), n=1,2,3,4,...

Here is the recurrence relation coded into Mathematica :

Clear[a]
a[0]=1;
a[n_]:=a[n]=a[n-1]/(n(n-2));

Now, we can get the first few terms of the solution:

n=5;
y1 = Sum[a[i]x^i,{i,0,n}]

Division by zero means we can't find a second solution!

Series substitution to find the second solution

We can find a second solution by direct series substitution, just as we did above. I'll show each step in the process here, since the Log function complicates things slightly. First, we need to define the second solution in terms of the first (for now we will not use the first solution, we will insert it later; this second solution comes from pg 291 Eq 24 in Boyce and DiPrima, and I am just assuming a=1 to simplify life!):

Clear[y1,dy1,ddy1,dy2,y2,dy2,ddy2,a,b,eq]
n=7;
y2 = y[x] Log[x] +(x^(-2)) Sum[b[i](x^i),{i,0,n}]

Calculate derivatives:

dy2 = D[y2,x]

ddy2 = D[y2,{x,2}]

Substitute into the DE:

eq = Collect[Simplify[x ddy2 +3 dy2-  y2],Log[x]]

The Log[x] term drops out since it multiplies the L[y1]=0. So we are left with:

eq = -((b[0])/(x^2))-((b[1])/(x^2))-((b[1])/x)-b[2]+3 b[3]-x b[3]+8 x b[4]-(x^2) b[4]+15 (x^2) b[5]-(x^3) b[5]+24 (x^3) b[6]-(x^4) b[6]+35 (x^4) b[7]-(x^5) b[7]+((2 y[x])/x)+2 y'[x]

Now we must differentiate the original solution y=(Subscript[y, 1]), and insert it in the differential equation:

y[x] = 1+(x/3)+((x^2)/24)+((x^3)/360)+((x^4)/8640)+((x^5)/302400)
y'[x] = D[y[x],x]

eq

Now we must collect powers of x:

eq = Expand[eq];
eq = Collect[eq,x]

Now we must solve for all the b[n], and we will set the arbitrary constant b[2] to be 29/36 for comparison with the second method below:

b[1] = 2;
b[0] = -2;
b[2] = 29/36;
Solve[ {- b[2]+3 b[3]+ 4/3==0,(1/4)-b[3]+8 b[4]==0,(1/45)-b[4]+15 b[5]==0,(1/864)-b[5]+24 b[6]==0,(1/25200)-b[6]+35 b[7]==0},{b[3],b[4],b[5],b[6],b[7]}]

y2 = y2 /.{b[3]->-(19/108),b[4]->-(23/432),b[5]->-(163/32400),b[6]->-(401/1555200),b[7]->-(3239/381024000)}

This is the second solution to the differential equation.

seriessubstitution = y2;

Reduction of Order to find the second solution

This method is simpler to implement than the last. It involved using reduction of order, so the most difficult aspect of the calculation is the squaring and inversion of a series, which is easy to do with a computer algebra system.

Clear[y1,y2,dy2,ddy2,a,b,eq]

y2 = y1[x] Integrate[ Exp[ -Integrate[p[x],x]]/y1[x]^2,x]

y1[x] = 1+(x/3)+((x^2)/24)+((x^3)/360)+((x^4)/8640)+((x^5)/302400)

y1sq[x] = Expand[y1[x]^2]

iny1sq[x] = Normal[Series[1/y1sq[x],{x,0,5}]]

y2 = y1[x] Integrate[Exp[-Integrate[3/x,x]] iny1sq[x],x]

Are these two solutions the same?

Simplify[Series[seriessubstitution,{x,0,3}]]

Simplify[Series[4 y2,{x,0,3}]]

I know which method I prefer, which one do you prefer?