Chapter 5
Series Solutions of
Ordinary Differential Equations

5.2 The Method of Frobenius

Symbolic Solution

The steps leading to the indicial equation in (5.19) and the recursion relations in (5.21) can be carried out symbolically in Mathematica and expressions for the expansion coefficients a[n] can be obtained, however complex the recursion relations.  This may prove useful if any of the functions P, Q or R are not simple algebraic expressions, resulting in complicated recursion relations that are difficult to solve in closed form.

R[x_,n_]:=Sum[r[k]x^k,{k,0,n}];
P[x_,n_]:=Sum[p[k]x^k,{k,0,n}];
Q[x_,n_]:=Sum[q[k]x^k,{k,0,n}];
y[x_,s_,n_]:=x^s Sum[a[k]x^k,{k,0,n}];

Equation (5.16) is then obtained by using Collect to group the terms according to powers of x.  The resulting expression  (the left-hand side of (5.16)) is called Frobenius:

Frobenius[n_]:=Collect[R[x,n] D[y[x,s,n],{x,2}] + 
     x^(-1)P[x,n] D[y[x,s,n],x] +
     x^(-2)Q[x,n] y[x,s,n],x];

The indicial equation and the recursion relations can now be extracted from Frobenius by collecting the coefficient of the appropriate power of x.  For example, the indicial equation is obtained from the coefficients of x^{s-2}, as in (5.16):

Coefficient[Simplify[Frobenius[0]],x^(s-2)]

Similarly, the recursion relation for a[n] in (5.22) can be determined from the coefficient of x^{s-2+n} in Frobenius by using Solve to obtain an expression for a[n] in terms of a[n-1],...,a[0].  For example, solving for a[1] yields an expression equivalent to that obtained by solving the first of Equations (5.23):

Solve[Coefficient[Simplify[Frobenius[1]],
     x^(s-1)]==0,a[1]];

Simplify[%]

Example 5.2

The steps used earlier in Mathematica for the general case of the method of Frobenius can be applied to solving the differential equation in (5.29).  The expression Frobenius is set up as before by defining the coefficient functions R, P and Q and the solution y:

R[x_]:=1-x^2;
P[x_]:=-2x^2;
Q[x_]:=2x^2;
y[x_,s_,n_]:=(x^s)Sum[a[k]x^k,{k,0,n}];

Frobenius[n_]:=Collect[R[x] D[y[x,s,n],{x,2}] + 
    x^(-1)P[x] D[y[x,s,n],x] + x^(-2)Q[x] y[x,s,n],x];

Frobenius can now be used either to generate the indicial equation and the recursion relations for any given n, as in (5.31) and (5.32).  Coefficient is invoked to obtain individual terms in Frobenius, e.g.

Factor[Coefficient[Expand[Frobenius[0]],x^(s-2)]]

Factor[Coefficient[Expand[Frobenius[1]],x^(s-1)]]

The solution to the differential equation can be obtained by first constructing the recursion relations,

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

The constants a[0] and a[1] have been left as arbitrary.  The series for the solution can now be displayed to any desired order.  For example,

y[x,0,10]

which is the same as (5.37).

5.3 Constructing Complementary Solutions

Example 5.3

The modifications to the Mathematica procedure used in Example 5.2 to account for the s-dependence in the a[n] are quite minimal.  We begin as before by constructing the series representation of the solution directly from the recursion relations, but now we retain the s-dependence in the recursion relations and declare s as an argument of the expansion coefficients:

a[0,s]=a0;
a[n_,s_]:=If[n>0,a[n-1,s]/(n+s)^2]

y[x_,s_,n_]:=x^s Sum[a[k,s]x^k,{k,0,n}]

The solution (5.52) produced by the direct application of the method of Frobenius is then obtained by evaluating y[x,s,n] with s set equal to zero.  Thus, up to terms of 10th order in x, we obtain the solution directly from the recursion relations as in Example 5.2:

y[x,s,10]/.s->0

The denominators corresponding to x^n in this series are readily verified to be equal to (n!)^2. The complementary solution (5.55) is obtained by performing the
construction in (5.44).  Again retaining terms to 10th order in y[x,s,n], we obtain

D[y[x,s,10],s]/.s->0

The verification that (5.55) yields the same expression as that shown above is left as an exercise.

Example 5.4

The procedures used in Mathematica to construct the complementary solution in Example 5.3 can be used here with little modification.  Here, we will use Mathematica to check that the complementary solution is indeed a solution of (5.63).  We first construct (5.77):

S[n_]:=Sum[1/i,{i,2,n}]

y[x_,n_]:=-x^2 Exp[-x] Log[x] + (1 + x - x^2) +
           x^2 Sum[(-1)^k S[k] (x^k/k!),{k,2,n}]

This solution is now substituted into the differential equation and the command 
Series is used to expand each term in the equation to the specified order of the expansion of the equation as a whole.  In the example below, the solution is expanded to 10th order.

Series[x D[y[x,10],{x,2}]+(x-1)D[y[x,10],{x,1}]+
     y[x,10],{x,0,10}]