Initializations

Off[General::"spell1"]

Series Solution about an Ordinary Point

Here we show the details of an example of obtaining the series solution about an ordinary point for a variable coefficient differential equation. This will be done by considering the example differential equation

y''-5x y = 0.

and searching for a series solution about x=0.

Solution

Classify the Point you are expanding about

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) = 0   and    q(x) = -5x

which means that x=0 is an ordinary point, since the functions p(x) and q(x) are both analytic at x=0.

Substitute into the DE and determine the recursion relation

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

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

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

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

The recurence relation is therefore:

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

Solving the Recursion Relation

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

a[0]=(Subscript[a, 0]);
a[1]=(Subscript[a, 1]);
a[2]=0;
a[n_]:=a[n]=((5a[n-3])/(n(n-1)));

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

n=10;
y =(∑_(i=0)^n)a[i](x^(i ))
Collect[y,{(Subscript[a, 0]),(Subscript[a, 1])}]

Pattern in the first solution

Let's look for the pattern in (Subscript[y, 1])=1+((5 (x^3))/6)+((5 (x^6))/36)+((25 (x^9))/2592)+.... If we can find it great--if not, we can still use what we have above as an approximation to the general solution. Let's see if Mathematica can do something useful in finding the pattern (patterns are sometimes easier to find when we do the calculation by hand). We get this solution by setting a[0]=1 and a[1]=0

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

n=20
Table[{i,a[i]},{i,0,n}]//TableForm

From this, we can ascertain the pattern as:

(Subscript[a, 3i])= ((5^i)/((3^i)(Underoverscript[i!∏, (k=1), i])(3k-1))), i=0,1,2,3,4,...

Let's check if we got it right:

n=5;
y1 =1+ Sum[((5^i)/((3^i)i!(∏_(k=1)^i)(3k-1)))(x^(3i)),{i,1,n}]

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

y1 =1+ Sum[((5^i)/((3^i)i!(∏_(k=1)^i)(3k-1)))(x^(3i)),{i,1,∞}]

And, a final check:

FullSimplify[D[y1,{x,2}]-5x y1==0]

Pattern in the second solution

Let's look for the pattern in (Subscript[y, 1])=x+((5 (x^4))/12)+((25 (x^7))/504)+((25 (x^10))/9072)+.... If we can find it great--if not, we can still use what we have above as an approximation to the general solution. Let's see if Mathematica can do something useful in finding the pattern (patterns are sometimes easier to find when we do the calculation by hand). We get this solution by setting a[0]=0 and a[1]=1

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

n=20
Table[{i,a[i]},{i,0,n}]//TableForm

From this, we can ascertain the pattern as:

(Subscript[a, 3i+1])= ((5^i)/((∏_(k=1)^i)(3k+1)(∏_(k=1)^i)(3k))), i=1,2,3,4,...

Let's check if we got it right:

n=5;
y2 =x+ Sum[ ((5^i)/(((∏_(k=1)^i)(3k+1))((∏_(k=1)^i)(3k))))(x^(3i+1)),{i,1,n}]

Note:

((∏_(k=1)^4)(3k+1))((∏_(k=1)^4)(3k))
(∏_(k=1)^4)(3k+1)(∏_(k=1)^4)(3k)

So always use brackets with products!

Let's see if Mathematica knows this function:

y2=x+ Sum[ ((5^i)/(((∏_(k=1)^i)(3k+1))((∏_(k=1)^i)(3k))))(x^(3i+1)),{i,1,∞}]

And, a final check:

FullSimplify[D[y2,{x,2}]-5x y2==0]

5.3: Series Solution about Ordinary Point involving Product of Series

Example

Here we want to solve a differential equation that requires the use of series and products of series. The sample differential equation is

Clear[y,y1,y2,x]
L[y_]:=D[y,{x,2}]+Exp[x]D[y,x]-y==0

L[v[x]]

Notice that we have a product of the solution v(x) and the exponential function. We will need to expand the exponential function is a series so we can collect powers of x.

We can use the result for series multiplication to help us here

(∑_(n=0)^∞)(Subscript[a, n])(((x-(Subscript[x, 0])))^n)×(∑_(n=0)^∞)(Subscript[b, n])(((x-(Subscript[x, 0])))^n)=(∑_(n=0)^∞)((∑_(k=0)^n)(Subscript[a, k])(Subscript[b, n-k]))(((x-a))^n).

Let's assume the solution looks like:

v(x)=(∑_(n=0)^∞)(Subscript[c, n])(x^n).

Then series manipulations leads to:

(d/(d x))v(x)=(∑_(n=0)^∞)n (Subscript[c, n])(x^(n-1))= (∑_(n=0)^∞)(n+1) (Subscript[c, n+1])(x^n).

((d^2)/(d (x^2)))v(x)= (∑_(n=0)^∞)(n+1)(n+2) (Subscript[c, n+2])(x^n).

(e^x)(d/(d x))v(x) =(∑_(n=0)^∞)(1/(n!))(x^n)× (∑_(n=0)^∞)(n+1) (Subscript[c, n+1])(x^n) = (∑_(n=0)^∞)((∑_(k=0)^n)(1/((n-k)!))(k+1) (Subscript[c, k+1]))(x^n).

Putting everything back into the differential equation leads to

(∑_(n=0)^∞)(n+1)(n+2) (Subscript[c, n+2])(x^n)+ (∑_(n=0)^∞)((∑_(k=0)^n)(1/((n-k)!))(k+1) (Subscript[c, k+1]))(x^n)-(∑_(n=0)^∞)(Subscript[c, n])(x^n)=0

(∑_(n=0)^∞)((n+1)(n+2) (Subscript[c, n+2])+ (∑_(k=0)^n)(1/((n-k)!))(k+1) (Subscript[c, k+1])-(Subscript[c, n]))(x^n)=0

(n+1)(n+2) (Subscript[c, n+2])+ (∑_(k=0)^n)(1/((n-k)!))(k+1) (Subscript[c, k+1])-(Subscript[c, n])=0,    n=0,1,2,3,...

The nice thing about doing this analysis is that we have access to the (Subscript[c, n])directly; we can immediately get the solution to any order we desire, and if we are extremely lucky we can recognize a pattern in the solution. Let's use Mathematica to help with evaluating the (Subscript[c, n]):

Clear[c];
num=10;
c[0]:=c0;
c[1]:=c1;
c[n_]:=((- (∑_(k=0)^(n-2))(1/((n-2-k)!))(k+1) c[k+1]+c[n-2])/((n-1)(n)));
(∑_(n=0)^num)c[n](x^n);
Expand[%];
eq1 = Collect[%,{c0,c1}]
y1[x_]=Coefficient[eq1,c0]
y2[x_]=Coefficient[eq1,c1]

It is difficult to recognize a pattern. If we want more terms we can very easily calculate them.

Here we see how the truncated solution approximates the complete solution (for (Subscript[y, 1])(x) in this case):

D[y1[x],{x,2}]+Exp[x]D[y1[x],x]-y1[x]
Series[%,{x,0,10}]

This should be zero, but it is not due to the truncation. Adding more terms to the solution will improve its accuracy in the interval of convergence for the series.

Example With Earlier Truncation

Sometimes the formal method we used above to get the (Subscript[c, n]) is cumbersome; this example shows how you can proceed by truncating earlier. It isn't as mathematically beautiful as the example above, but it will work.

This example shows the manipulations you need to perform to determine the series solution to y'' - cos(x) y = 0. If we want the fundamental set of soutions to have four terms each, we will need at least eight terms (possibly more) in the assumed series solution. This is because we expect the solution to break into two parts, one with arbitrary coefficient (Subscript[a, 0]) and one with arbitraty coefficient (Subscript[a, 1]). If some coeffcients are zero, we may need to go to higher terms.

Since cosine has a convergent Taylor Series for all x, the point x=0 is an ordinary point. Our series solutions will converge for all x.

Clear[a,x,y,i]
n=8;
y = Sum[a[i]x^i,{i,0,n}]

ddy = D[D[y,x],x]

Notice that the second derivative only has terms to order six. Keeping in mind that our goal is to calculate (Subscript[a, n]), n=0,1,2,...,8 we should not be concerned about this.

We will need to expand the cosine to the same order:

cs = Normal[Series[Cos[x],{x,0,n}]]

We now can substitute into the differential equation:

eq1 = ddy - cs y

We now must expand and collect in powers of x:

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

Notice that we get higher order terms than (x^8); these terms are in fact incomplete, since they are missing some contributions since we truncated the series. To see this, think of expanding to order 1, and then to order 2:

n0=1;
y0 = Sum[a[i]x^i,{i,0,n0}];
ddy0 = D[D[y0,x],x];
cs0 = Normal[Series[Cos[x],{x,0,n0}]];
eq0 = ddy0 - cs0 y0

If we set this equal to zero, and then set the coefficients of x equal to zero, we have a[0] = 0 and a[1] = 0! What happened? Let's look at keeping more terms in the expansion:

n0=2;
y0 = Sum[a[i]x^i,{i,0,n0}];
ddy0 = D[D[y0,x],x];
cs0 = Normal[Series[Cos[x],{x,0,n0}]];
eq0 = ddy0 - cs0 y0;
eq0 = Expand[eq0];
eq0 = Collect[eq0,x]

These equations tell us that 2a[2]- a[0] = 0 and a[1] = 0! See how including extra terms has changed the results? We now have all the terms we need to calculate the (x^0) term accurately, but the (x^1) term is still not accurate.

Returning to our original train of thought, we now must set each coefficient of x equal to zero, and use those relations to solve for the a[n], n=0,1,2,...,8; a[0] and a[2] will be arbitrary:

Do[ sol = Solve[Coefficient[eq1,x,i-2]==0,a[i]];
a[i] = a[i] /. First[sol], {i,2,n}]

y

Collect[y,{a[0],a[1]}]

Let's double check our result, by redoing it to higher order, and ensuring the coefficients don't change:

Clear[a,x,y,i]
n=10;
y = Sum[a[i]x^i,{i,0,n}];
ddy = D[D[y,x],x];
cs = Normal[Series[Cos[x],{x,0,n}]];

eq1 = ddy - cs y;
eq1 = Expand[eq1];
eq1 = Collect[eq1,x];
Do[ sol = Solve[Coefficient[eq1,x,i-2]==0,a[i]];
a[i] = a[i] /. First[sol], {i,2,n}]

Collect[y,{a[0],a[1]}]

We see that the terms to order (x^8) did not change, and we can be confident we determined them accurately. Another way to check our accuracy is to substitute our series solution into the differential equation and see what order of x is left:

Expand[eq1]

Notice the highest power of x in the above is 9; therefore we can be confident that our series solution is accurate to order 8.

Example: Better than DSolve?

Notice that this is a very powerful technique, for with a simple modification it can determine a series solution for differential equations which Mathematica cannot solve:

Clear[v,x]
DSolve[v''[x] +Exp[x]v'[x]- Cos[x] v[x]==0,v[x],x]

Clear[a,x,y,i]

n=10;
y = Sum[a[i]x^i,{i,0,n}];
dy = D[y,x];
ddy = D[D[y,x],x];
cs = Normal[Series[Cos[x],{x,0,n}]];
es = Normal[Series[Exp[x],{x,0,n}]];

eq1 = ddy +es dy- cs y;
eq1 = Expand[eq1];
eq1 = Collect[eq1,x];
Do[ sol = Solve[Coefficient[eq1,x,i-2]==0,a[i]];
a[i] = a[i] /. First[sol], {i,2,n}]

Collect[y,{a[0],a[1]}]

Verify to what power of x our solution is valid:

Expand[eq1]

We see that the solution is valid to order 8; notice that in this case we had expanded originally to order 10, but our solution is only valid to order 8.

IVP Example

Here is the general solution from the example above:

ygeneral[x_]=Collect[y,{a[0],a[1]}]
Clear[y,y0]

If we have some initial conditions, we can determine the solution to an initial value problem.

condition1[y_,a_,y0_]:=y==y0/.x->a
condition2[y_,a_,y0_]:=D[y,{x,1}]==y0/.x->a

Here is what the commands do in general:

condition1[y[x],a,y0]
condition2[y[x],a,y0]

Here are some specific initial conditions:

condition1[y[x],0,1]
condition2[y[x],0,1/2]

We want to use these conditions on the general solutions, and solve for the a[0] and a[1]:

eq1=condition1[ygeneral[x],0,1]
eq2=condition2[ygeneral[x],0,1/2]

ygeneral[x_]=Simplify[ygeneral[x]/.a[0]->1/.a[1]->1/2]

Plot[ygeneral[x],{x,-2,2},AxesLabel->{"x","y"},ImageSize->{500,Automatic},PlotStyle->{Thick}]

This solution is only valid in a region about (Subscript[x, 0])=0, since we truncated the series. To estimate the region the solution is valid in, you should do the calculation again with earlier truncation, and see where that solution agrees with the one above.

5.4: Example of Failure of Method When (Subscript[x, 0]) is Not an Ordinary Point

What happens if we are not working around an ordinary point? Will the Taylor series method fail? Let's explore!

Clear[y,a,x]
L[y_]:=x^2*D[y,{x,2}]-2y
L[y[x]]

In this case, x=0 is a regular singular point since p(x)=0, q(x)=-2/(x^2) are not analytic at x=0, but  x p(x)=0, (x^2)q(x)=-2 are analytic at x=0.

Let's substitute a truncated series about (Subscript[x, 0])=0  into the differential equation, and see what we can learn:

n=5
yassumed[x_]:= Sum[a[i]x^i,{i,0,n}]
yassumed[x]

Substitute into the differential equation:

deq1=L[yassumed[x]]

Collect powers of x

deq2 = Expand[deq1]
deq3 = Collect[deq2,x]

Strip off the coefficients, to get equations the a[n] must satisfy

eqs=Table[Coefficient[deq3,x,i]==0,{i,0,n}]

vars=Table[a[i],{i,0,n}]

sol=Solve[eqs,vars]

Notice that all the a[n] must be zero except a[2], which is the only arbitrary constant.

The solution we find is therefore

yassumed[x_]=yassumed[x]/.Flatten[sol]

The method only found one solution. The problem is the other solution looks like 1/x, which is not analytic at (Subscript[x, 0])=0, and our assumed form for the series solution was analytic at (Subscript[x, 0])=0.

Summary:

The method for ordinary points will not work for singular points because the solution is often not analytic about a singular point, and so no Taylor series solution exists about that point.

No Taylor Series Solution exists in general, but we will stll get series solutions!! (Not all series are Taylor series)

5.4: Example of Failure of Method When (Subscript[x, 0]) is Not an Ordinary Point: Hand Solution

Let's revisit the same problem we just finished, but work everything through by hand. We will arrive at the same conclusions.

Clear[y,a,x]
L[y_]:=x^2*D[y,{x,2}]-2y
L[y[x]]==0

Here is our assumed Tayor series solution:

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

Substitute into the differential equation:

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

For this to be true for all values of x, we must have each coefficient of a power of  x equal to zero:

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

Notice that this is NOT what we have seen before with recursion relations! But mathematically, it still makes sense. Let's continue, evaluating at n=0,1,2,3,...:

n=0: -2 (Subscript[a, 0])= 0, so (Subscript[a, 0])= 0.
n=1: -2 (Subscript[a, 1])= 0, so (Subscript[a, 1])= 0.
n=2: 2(Subscript[a, 2])-2 (Subscript[a, 2])= 0, so 0= 0. This means (Subscript[a, 2])is arbitrary.

For higher n>2, we have (cnostant)(Subscript[a, n])=0, so (Subscript[a, n])=0. Therefore, the only arbitrary constant is (Subscript[a, 2]), and the only solution we find is

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

Which is the same as we found above.

Summary:

The method for ordinary points will not work for singular points because the solution is often not analytic about a singular point, and so no Taylor series solution exists about that point.

No Taylor Series Solution exists in general, but we will stll get series solutions!! (Not all series are Taylor series)