Series Solutions about Ordinary Points

What we want to be able to accomplish when we solve variable coefficient linear differential equations:

Find a closed form for two solutions in the fundamental set of solutions: (Subscript[y, 1]):=(∑_(n=0)^∞) (((Subscript[a, n])(x-(Subscript[x, 0])))^n) and (Subscript[y, 2]):=(∑_(n=0)^∞) (((Subscript[a, n])(x-(Subscript[x, 0])))^n)

Identify the closed form as known functions

In general, neither of these may be possible! Frequently, we can identify a pattern, but not what function it is that the sum represents. Sometimes we can't even figure out the pattern!

Example  y''[x]+x y[x]=0 about (Subscript[x, 0])=0

Recursion Relations & Truncated Solution

Since p(x)=0 and q(x)=x, which are both analytic for all real x, the point (Subscript[x, 0])=0 is an ordinary point.

We expect to find two solutions of the form (∑_(n=0)^∞) (Subscript[a, n])(x^n) (the fundamental set of solutions).

Let's start with that assumption and see what we find out!

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

Substitute into the DE
(∑_(n=0)^∞) n(n-1)(Subscript[a, n])(x^(n-2))+x (∑_(n=0)^∞) (Subscript[a, n])(x^n)=0

Multiply through the x
(∑_(n=0)^∞) n(n-1)(Subscript[a, n])(x^(n-2))+(∑_(n=0)^∞) (Subscript[a, n])(x^(n+1))=0

Relabel to get all terms with an (x^m)
(∑_(m=-2)^∞) (m+2)(m+1)(Subscript[a, m+2])(x^m)+(∑_(m=1)^∞) (Subscript[a, m-1])(x^m)=0

Get all summations to start at m=1
0+0+2(Subscript[a, 2])(x^0)+(∑_(m=1)^∞) (m+2)(m+1)(Subscript[a, m+2])(x^m)+(∑_(m=1)^∞) (Subscript[a, m-1])(x^m)=0
2(Subscript[a, 2])(x^0)+(∑_(m=1)^∞)[ (m+2)(m+1)(Subscript[a, m+2])+(Subscript[a, m-1])](x^m)=0

For this to be true for all values of x
2(Subscript[a, 2])=0
(m+2)(m+1)(Subscript[a, m+2])+(Subscript[a, m-1])=0, m=1,2,3,4,...

The recurrence relations are
(Subscript[a, 0])=unspecified
(Subscript[a, 1])=unspecified
(Subscript[a, 2])=0
(Subscript[a, m+2])=((-(Subscript[a, m-1]))/((m+2)(m+1))), m=1,2,3,4,...

which we can rewrite as
(Subscript[a, 0])=unspecified
(Subscript[a, 1])=unspecified
(Subscript[a, 2])=0
(Subscript[a, n])=((-(Subscript[a, n-3]))/(n(n-1))), n=3,4,...

At this point, we can generate truncated series approximations to the solutions using Mathematica and the recursion relations we have found.

Since it is just as easy to get the fundamental set of solutions as the solution to the IVP, I will focus on the fundamental set of solutions. We can easily use the initial conditions to determine specific values of a0 and a1 if we want to.

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

Here is the solution:

(Subscript[y, gen])[x_,k_]:=(∑_(n=0)^k) a[n](x^n)
(Subscript[y, gen])[x,20]

Notice that we have a series solution times a0 and a series solution times a1; these two series solutions form a fundamental set of solutions.

Collect[(Subscript[y, gen])[x,20],{a0,a1}]
y1[x_]=Coefficient[(Subscript[y, gen])[x,20],a0]
y2[x_]=Coefficient[(Subscript[y, gen])[x,20],a1]

truncatedsolplot=Plot[{y1[x],y2[x]},{x,-5,10},PlotRange->{{-5,10},{-5,5}},PlotStyle->{{Thick,Dashed}},BaseStyle->20]

Find Pattern in Recursion Relations

The first step is to write out the first few (Subscript[a, i]).

It's a good idea to leave the (Subscript[a, 0]) and (Subscript[a, 1]) unspecified at this stage even if you have an IVP since you will then find a fundamental set of solutions.

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

(Subscript[a, 0])=unspecified
(Subscript[a, 1])=unspecified
(Subscript[a, 2])=0

n=3: (Subscript[a, 3])=((-(Subscript[a, 0]))/(3·2))=((-(Subscript[a, 0])·1)/(3!))
n=4: (Subscript[a, 4])=((-(Subscript[a, 1]))/(4·3))=((-(Subscript[a, 1])·2)/(4·3·2))=((-(Subscript[a, 1])·2)/(4!))
n=5: (Subscript[a, 5])=0

n=6: (Subscript[a, 6])=((-(Subscript[a, 3]))/(6·5))=(((Subscript[a, 0])·1)/(6·5·3·2))=(((Subscript[a, 0])·1·4)/(6·5·4·3·2))=(((Subscript[a, 0])·1·4)/(6!))
n=7: (Subscript[a, 7])=((-(Subscript[a, 4]))/(7·6))=((-5)/(7·6·5))((-(Subscript[a, 1]))/(4!))=(((Subscript[a, 1])·2·5)/(7!))
n=8: (Subscript[a, 8])=0

n=9: (Subscript[a, 9])=((-(Subscript[a, 6]))/(9·8))=((-7)/(9·8·7))(((Subscript[a, 0])·1·4)/(6!))=((-(Subscript[a, 0])·1·4·7)/(9!))
n=10: (Subscript[a, 10])=((-(Subscript[a, 7]))/(10·9))=((-8)/(10·9·8))(((Subscript[a, 1])·5)/(7!))=((-(Subscript[a, 1])·2·5·8)/(10!))
n=11: (Subscript[a, 11])=0

These are in groups of three, due to the zeros that are appearing.

We need to define a variable k that we can sum over that will contain the groups of three, rather than using n.

You can make a variety of choices, but the following will get the groups of three part:

y(t)=(∑_(n=0)^∞)(Subscript[a, n])(x^n)
=(∑_(k=0)^∞)((Subscript[a, 3k])(x^(3k))+(Subscript[a, 3k+1])(x^(3k+1))+(Subscript[a, 3k+2])(x^(3k+2)))

k=0
3k=0: (Subscript[a, 3k])=unspecified
3k+1=1: (Subscript[a, 3k+1])=unspecified
3k+2=2: (Subscript[a, 3k+2])=0

k=1
3k=3:  (Subscript[a, 3k])=((-1)/(3!))(Subscript[a, 0])
3k+1=4: (Subscript[a, 3k+1])=((-2)/(4!))(Subscript[a, 1])
3k+2=5: (Subscript[a, 3k+2])=0

k=2
3k=6: (Subscript[a, 3k])=((1·4)/(6!))(Subscript[a, 0])
3k+1=7: (Subscript[a, 3k+1])=((2·5)/(7!))(Subscript[a, 1])
3k+2=8: (Subscript[a, 3k+2])=0

k=3
3k=9: (Subscript[a, 3k])=((-1·4·7)/(9!))(Subscript[a, 0])
3k+1=10: (Subscript[a, 3k+1])=((-2·5·8)/(10!))(Subscript[a, 1])
3k+2=11: (Subscript[a, 3k+2])=0

Notice that this means we can sum over k=0,1,2,3,4,... and we will have (Subscript[a, 3k+2])=0.

We just need to find the patterns for (Subscript[a, 3k]) and (Subscript[a, 3k+1]).

Start by working with the most complicated one first, and using a product for things we can't immediately identify as factorials.

k=3
3k=9: (Subscript[a, 3k])=(((((-1))^k)(∏_(j=0)^(k-1))(3j+1))/((3k)!))(Subscript[a, 0])
3k+1=10: (Subscript[a, 3k+1])=(((((-1))^k)(∏_(j=0)^(k-1))(3j+2))/((3k+1)!))(Subscript[a, 1])
3k+2=11: (Subscript[a, 3k+2])=0

Holy crap! I think we've got the pattern! You should check this works for k=0,1,2,3.

The products can be written in terms of gamma functions, but we don’t need to know that since we will use Mathematica to simplify.

(∏_(j=0)^(k-1))(3j+1)

(∏_(j=0)^(k-1))(3j+2)

Notice that for k=0 the product is a bit weird since it starts at zero and ends at minus one, but Mathematica works it out OK.

(∏_(j=0)^(0-1))(3j+2)

To be safe, we can just pull out the k=0 term in the sum.

If we want to determine the fundamental set of solutions, leave (Subscript[a, 0]) and (Subscript[a, 1]) unspecified and collect terms involving (Subscript[a, 0]) and terms with (Subscript[a, 1]):

(Subscript[a, 3k])=(((((-1))^k)(∏_(j=0)^(k-1))(3j+1))/((3k)!))(Subscript[a, 0]), k>0
(Subscript[a, 3k+1])=(((((-1))^k)(∏_(j=0)^(k-1))(3j+2))/((3k+1)!))(Subscript[a, 1]), k>0
(Subscript[a, 3k+2])=0

The general solution is
y(t)
=(∑_(n=0)^∞)(Subscript[a, n])(x^n)
=(∑_(k=0)^∞)((Subscript[a, 3k])(x^(3k))+(Subscript[a, 3k+1])(x^(3k+1))+(Subscript[a, 3k+2])(x^(3k+2)))
=(Subscript[a, 0])+(Subscript[a, 1])x+(∑_(k=1)^∞)((Subscript[a, 3k])[k]*(x^(3k))+(Subscript[a, 3k+1])[k]*(x^(3k+1)))
=(Subscript[a, 0])+(Subscript[a, 1])x+(∑_(k=1)^∞)((((((-1))^k)(∏_(j=0)^(k-1))(3j+1))/((3k)!))(Subscript[a, 0])*(x^(3k))+(((((-1))^k)(∏_(j=0)^(k-1))(3j+2))/((3k+1)!))(Subscript[a, 1])*(x^(3k+1)))
=(Subscript[a, 0]){1+(∑_(k=1)^∞)((((((-1))^k)(∏_(j=0)^(k-1))(3j+1))/((3k)!))*(x^(3k)))}+(Subscript[a, 1]){x+(∑_(k=1)^∞)((((((-1))^k)(∏_(j=0)^(k-1))(3j+2))/((3k+1)!))*(x^(3k+1)))}
=(Subscript[a, 0])(Subscript[y, 1])+(Subscript[a, 1])(Subscript[y, 2])

A Fundamental Set of Solutions is
(Subscript[y, 1])=1+(∑_(k=1)^∞)((((((-1))^k)(∏_(j=0)^(k-1))(3j+1))/((3k)!))*(x^(3k)))
(Subscript[y, 2])=x+(∑_(k=1)^∞)((((((-1))^k)(∏_(j=0)^(k-1))(3j+2))/((3k+1)!))*(x^(3k+1)))

Determine the Underlying Functions in the Fundamental Set of Solutions

Mathematica can often help determine the underlying function in each case:

y1[x_]=Simplify[1+(∑_(k=1)^∞)((((((-1))^k)(∏_(j=0)^(k-1))(3j+1))/((3k)!))*(x^(3k)))]
y2[x_]=Simplify[x+(∑_(k=1)^∞)((((((-1))^k)(∏_(j=0)^(k-1))(3j+2))/((3k+1)!))*(x^(3k+1)))]

Let’s compare to the truncated solutions we found earlier.

solplot=Plot[{y1[x],y2[x]},{x,-5,10},PlotRange->{{-5,10},{-5,5}},PlotStyle->{Thick},BaseStyle->20];

Show[truncatedsolplot,solplot]

Let's see if Mathematica can solve this IVP

DSolve[{y''[x]+x y[x]==0},y[x],x]

Mathematica’s solution looks different, but the Airy functions it returns are the same are the same as the solutions we found with our series, up to a constant factor (what is important is that they form a fundamental set of solutions).

For a given IVP, both fundamental sets of solutions return the same function, just expressed differently.

DSolve[{y''[x]+x y[x]==0,y[0]==2,y'[0]==-1},y[x],x]

yIVP[x_]=a0*y1[x]+a1*y2[x]/.{a0->2,a1->-1}

Plot[(1/6) (-3 (((-1))^(2/3)) (3^(1/3)) AiryAi[(((-1))^(1/3)) x] Gamma[(1/3)]+(((-1))^(2/3)) (3^(5/6)) AiryBi[(((-1))^(1/3)) x] Gamma[(1/3)]+6 (3^(2/3)) AiryAi[(((-1))^(1/3)) x] Gamma[(2/3)]+6 (3^(1/6)) AiryBi[(((-1))^(1/3)) x] Gamma[(2/3)]),{x,-5,10},PlotRange->{{-5,10},{-5,5}},PlotStyle->{Thick},BaseStyle->20]

Plot[{yIVP[x]},{x,-5,10},PlotRange->{{-5,10},{-5,5}},PlotStyle->{Thick},BaseStyle->20]