The Power Method with Summa.m

Examples of manipulation of infinite series in ODEs

The main reason that prompted me to write the simple, almost naive, package Summa.m was to allow my lazy self to check, without significant effort, intermediate steps in the application of the power method to ordinary differential equations (ODE). Reindexing, collecting and simplification of sums are quite straigthforward but nonetheless tedious when performed by hand. Incidentally, such manipulations extend to other fields, different from ODE solving and a few examples are given in the notebook "Application Examples". Despite the fact that the package be rather crude (it allows only for simple manipulations, without exerting any sort of control on user input - so beware of what you ask it to do!) there are several possible applications. Starred sections are still partially unfinished.

Initialization

First of all we load the package

<<Summa.m

Here's a list of the functions and tags implemented in the package:

?Summa`*

What's new in version 2: setting the environment

After loading the package, built-in evauation of the function Sum is temporarily disabled. All symbolic and infinite sums (that is, sums with one or both non numeric endpoints) will be immediatly turned into Summa objects, allowing them to escape Mathematica's evaluator. They will nonetheless continue to be formatted as ordinary Sum(s)

expr = Sum[k, {k, 1, n}]TableForm[{{InputForm[expr], OutputForm[expr]}}, 
      TableHeadings -> {None, {"InputForm", "OutputForm"}}]

This feature allows us to input sum(mma)s with the traditional notation used by Sum(s). They will be automatically converted into Summas and won't be evaluated.

(∑_(k = 0)^n) k

InputForm[%]

When we need to enable ordinary Sum evaluation we can do it locally (with respect to a cartain expression)

EnableSumEvaluation[expr]

or globally, restoring the evaluation of Sum for all the input to come, by evalating the following command

EnableSumEvaluation[]

We can see that every new Sum is treated as expected in a standard Mathematica fashion. There's a catch, though: previously inputted Sum[s] were turned into Summa[s] as soon as they were entered. So they need to be converted into Sum[s] by a direct replacement, if we want to see any effect.

Sum[k, {k, 1, n}]

expr /. Summa -> Sum

The following disables Sum evaluation.

DisableSumEvaluation[]

A note on SumCoefficient

Version 3 understands this syntax

SumCoefficient[expr, x^k]

While on version 4 it is possible to use this form:

SumCoefficient[%, x, k]

So, if you are getting error messages while using getCoefficient, try to switch syntax: you are using the wrong one for your version. This will be fixed in future released of Summa. I apologize for the inconvenient.

Power method: ordinary points

Introduction to the power method

The power method is based on the following theorem of existence:

Theorem (existence and convergence of power solutions): Any ordinary differential equation in the form
	(y^((n)))(x)+s(x) (y^((n - 1)))(x)+p(x) y'(x)+q(x) y(x)=f(x)
whose variable coefficients s(x), ..., p(x), q(x) and f(x) allow for a power series representation in a neighbor of the point (Subscript[x, 0]), has a power series solution y(x)
	\(TraditionalForm
      \`y(x) = ∑\+\(k = 0\)\%∞\(\( a\_k\)(x - x\_0)\)\^k\)
whose convergence radius is as big as the smallest convergence radius of the series for the coefficients s(x),..., q(x) and f(x) seen as functions of a complex variable. An ordinary differential equation of the form α(x)(y^((n)))(x)+β(x) (y^((n - 1)))(x)+...+γ(x) y(x)=f(x) can be translated into the previous form by dividing both sides by α(x), provided that the coefficient of the highest derivative is not zero in (Subscript[x, 0]).

Finding a power series solution to a differential equation means finding the general form for the coefficient (Subscript[a, k]). This can be done by plugging the power series for y[x] and its derivatives into the equation (along with the power series for every variable coefficient, carrying out whatever calculation is needed to reduce the equation to a power series equated to zero, and by applying the principle of identity of polynomials to find the conditions to be imposed over the sum's coefficients (Subscript[a, 0]),  (Subscript[a, 1]),  (Subscript[a, 2]),   ...  , (Subscript[a, k]),  ... . The first n relations (with n being the order of the differential equation) will need the specification of the initial conditions. Form the definition of power series it is easily shown that (Subscript[a, 0]) = y((Subscript[x, 0])), (Subscript[a, 1]) = y' ((Subscript[x, 0])), (Subscript[a, 2]) = y'' ((Subscript[x, 0]))/2!, ..., (Subscript[a, j]) = (y^((j)))/j!

1. Solving y' == y

Let's see how the method applies to the trivial ordinary differential equation

eq = y'[x] - y[x] == 0;

This particular ODE's coefficients are simple constants (+1 and -1) so there is nothing to worry about the applicability of the power method. We choose to find a power series in a neighbor of x=0, thus limiting ourselves to a McLaurin series.We plug the power series and its first derivative into the original equation to get

eq /. {y[x] -> (∑_(k = 0)^∞) (Subscript[a, k]) (x^k), 
        y'[x] -> 
          SumD[(∑_(k = 0)^∞) (Subscript[a, k]) (x^k), x]}

Since we want to combine both sums into a single sum, we first need to change the index of summation with SumExponentAlign so that the terms will show the same power of x,and then we extract part of one of the two sums to make every index start form the same point. k=0, in this case. SumSplit is what we need for this task.

SumExponentAlign[%, x, k -> k]

SumSplit[%, 0]

As an aside, note that nothing has changed, apart from the starting index of the first sum. The reason is that, being the derivative of a power series, its first term - the derivative of the constant term of the original series - is zero. Now that both sums have the same endpoints and involve the same power of x we can proceed in merging them into a single structure

SumSimplify[%]

We're almost done. The sum is identically zero iff the coefficients of the power series are zero for every k. We have a relation to be imposed over the coefficents that can be solved for the coefficient with the highest index. On version 3 we use this:

SumCoefficient[%, x^k] == 0

On version 4 this form is better:

SumCoefficient[%, x, k] == 0

Solve[%, (Subscript[a, k + 1])]

Now look at what we've got here: the next coefficient of the series expansion of the solution is given by the last computed coefficient divided by its index. We have to specify one coefficient to get the others, and this is just a way to specifiy which particular solution to pick among the infinitely many solutions of a first order differential equation.

We can express our solution in term of the first coefficient and, since y(0) = (Subscript[a, 0]) this is equivalent to specifying the initial condition (Subscript[y, 0]). Let's see what happens to the first few coefficients and if we can find a closed form expression for the general coefficient (Subscript[a, n])

\(TraditionalForm\`\(	
    a\_0 = \(y\_0
	a\_1 = 
        \(\(a\_0\ \((1\ 1\/1)\) = y\_0\)
	a\_2 = 
          \(\(a\_1\/2 = \(a\_0\/\(2\ 1\) = y\_0\/\(2!\)\)\)
	a\_3 = 
            \(\(a\_2\/3 = \(a\_0\/\(3\ 2\) = y\_0\/\(3!\)\)\)
	…

                	a\_n = \(a\_\(n - 1\)\/n = y\_0\/\(n!\)\)\)\)\)\)\)\)

Our power series solution is thus

(∑_(n = 0)^∞) ((Subscript[y, 0])/(n!)) (x^n) // BringOut

a sum converging for every x and a Sum that Mathematica knows very well. In fact locally enabling the built-in Sum evaluation will turn it into an exponential function, which is the known solution of the differential equation.

EnableSumEvaluation[%]

We have thus shown that the power series method has allowed us to compute the solution to the simple ODE y'==y that is validated by our previous knowledge.

2. Solving y'' == y

Let's use the power method to solve a second order linear ODE. Again we attack a simple problem whose solution is well known.

eq = y''[x] + y[x] == 0;

This ODE's coefficients are simple constants (+1 and -1) and we can apply the power method without problems. We choose to find a power series in a neighbor of x=0. Note that we can explicit the first two terms of the second derivative (they are zero) to get a neater expression

{ys = (∑_(k = 0)^∞) (Subscript[a, k]) (x^k), 
      y2s = SumSplit[
          SumD[(∑_(k = 0)^∞) (Subscript[a, k]) (x^k), {x, 2}], 2]}

If we plug the power series for our tentative solution y(x) and its second derivative  into the equation, we get

eq /. {y[x] -> ys, y''[x] -> y2s}

Aligning the exponents of the powers of x has the effect of making the indices of both sums equal. A simplification will give us a single sum

SumExponentAlign[%, x, k -> k] // SumSimplify

For this equality to be true, it is necessary (and sufficient) that the coefficient of (x^k)be zero for every (non negative) k. Equating this coefficient to zero and solving for (Subscript[a, k + 2])yields the recursion relation we're after

SumCoefficient[%, x^k] == 0

On version 4 we'd use SumCoefficient[%,x,k]==0

Solve[%, (Subscript[a, k + 2])] /. k -> k - 2 // Simplify

This recursion relation will generate every coefficient of the power series that satisfies the differential equation, once the first two coefficients are known. It is easy to see that their values are related to the values of the solutions in x=0, namely y(0) = (Subscript[a, 0])and y' (0) = (Subscript[a, 1]).

Let's see if we can find a closed form expression for the general coefficient (Subscript[a, n]). We can identify two separate sequences whose numerical coefficients are related to the factorial of the index

({{(Subscript[a, 0]) = (Subscript[y, 0]),  , (Subscript[a, 1]) = (Subscript[yp, 0])}, {(Subscript[a, 2]) = 
                -((Subscript[a, 0])/2) = 
                  -((Subscript[y, 0])/2) = -((Subscript[y, 0])/(2!)),  , (Subscript[a, 3]) = 
                -((Subscript[a, 1])/6) = 
                  -((Subscript[yp, 0])/6) = -((Subscript[yp, 0])/(3!))}, {(Subscript[a, 4]) = -((Subscript[a, 2])/12) = ((Subscript[y, 0])/24) = ((Subscript[y, 0])/(4!)),  , (Subscript[a, 5]) = 
                -((Subscript[a, 3])/20) = ((Subscript[yp, 0])/120) = ((Subscript[yp, 0])/(5!))}, {(Subscript[a, 6]) = 
                -((Subscript[a, 4])/30) = 
                  -((Subscript[y, 0])/720) = -((Subscript[y, 0])/(6!)),  , (Subscript[a, 7]) = 
                -((Subscript[a, 5])/42) = 
                  -((Subscript[yp, 0])/5040) = -((Subscript[yp, 0])/(7!))}, {…,  , …}})

The pattern is clearly identified and a little background in calculus allows us to spot the coefficients of the McLaurin expansion of the elementary trigonometric functions Cos and Sin.

({{(Subscript[a, 2 n]) = (((((-1))^n) (Subscript[y, 0]))/((2 n)!)),  , (Subscript[a, 2 n + 1]) = 
                (((((-1))^n) (Subscript[yp, 0]))/((2 n + 1)!))}})

We can then write the power series solution in the following way:

\(TraditionalForm\`\(	
    y(x) = 
      y\_0\ \(∑\+\(n = 0\)\%∞\(\((\(-1\))\)\^n\ 
                x\^\(2\ n\)\)\/\(\((2\ n)\)!\)\) + 
        yp\_0\ \(∑\+\(n = 0\)\%∞\(\((\(-1\))\)\^n\ 
                x\^\(2\ n + 1\)\)\/\(\((2\ n + 1)\)!
                \)\)
\( = 
            \(y\_0\) \(cos(x)\) + \(yp\_0\) \(sin(x)\)\)\)\)

Taken separately, the two sums are particular solutions of the differential equation; the first one is the solution with initial conditions y(0) = 1and y' (0) = 0, while the second one is a solution for y(0) = 0 and y' (0) = 1. Since the first power series contains only even powers of x and the latter only odd powers, they are also  linearly independent. If we consider (Subscript[y, 0])and (Subscript[yp, 0])as arbitrary constants, then we have found the general solution to our ODE.

Mathematica can verify for us that we have found the correct solution. We need to enable built-in Sum evaluation before calling DSolve, though. By wrapping EnableSumEvaluation around DSolve we can show that the particular solution with initial onditions y(0) = (Subscript[y, 0])and y' (0) = (Subscript[yp, 0])is

DSolve[{
y''[x] + y[x] == 0, 
          
y[0] == yo, 
          y'[0] == ypo
}, y[x], x] // 
      EnableSumEvaluation

and has the same form as the well known general solution found for a constant coefficient ODE, and by NDSolve

DSolve[y''[x] + y[x] == 0, y[x], x] // EnableSumEvaluation

3. Solving (4-x^2)y'[x]+y[x]==0

Most functions in the Summa package can operate directly over equations in the form lhs==rhs, as we shall see n the next example. We can obviosly work by addressing only the left hand side, as in this example. Here is the lhs of the ODE to solve:

eq=(4-x^2)y'[x]+y[x];

We can find a power series solution as long as we stay clear form the points x=-2 and x=+2. The tentative power series solution and its derivatives, expanded in a neighbor of x=0, are

{ys = Sum[(Subscript[a, n]) x^n, {n, 0, Infinity}], 
        
y1s = SumD[ys, x], 
        
y2s = SumD[y1s, x]} // ColumnForm

By plugging them into the equation we get

eqs=eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}

Let's expand the products and bring the factors inside the sums

eqs=BringIn[SumExpand[eqs]]

It would be better to have all the power in the form (x^n): To achieve this goal we can use the function SumExponentAlign (previously named AlignExponent) in the following way (the third argument is in the form (name of the exponent)->(form into which has to be cast) )

eqs=SumExponentAlign[eqs,x,n->n]

Then we try to get sums that have the same starting value, without tampering with the indexes any more: we use SumSplit (BreakSum) to get the exceeding terms out of the sums. This command will make all the sums starting from 1

SumSplit[eqs,1]

Let's simplify this expression, compacting all the sum into a single one. SumSimplify (the late SimplifySum) does this for us:

res=SumSimplify[%]

Now all we have to do is to apply the principle of identity of polynomials. This can be done as follows: 
We first suppress the Sum[_ , it] structure with SumToTerm (the procedure formerly known as GetElement)

eqRicors=SumToTerm[res]

And then we can equate the coefficients of (x^n) and (x^0) to zero, solving for the unknown coefficients (Subscript[a, n]). It may be necessary to translate the index variable and simplify the result.

{cn} = 
      Solve[SumCoefficient[eqRicors, x^n] == 0, (Subscript[a, n + 1])] /. 
          n -> n - 1 // Simplify

{c1} = Solve[SumCoefficient[eqRicors, x, 0] == 0, (Subscript[a, 1])]

Alternatively we could have used the built-in function Coefficient.

{c1} = 
      Solve[(Coefficient[res, x, 0] /. _Summa -> 0) == 0, (Subscript[a, 1])]

Now we just have to construct the series, given the recursion relation.

\(TraditionalForm\`\(	
    a\_0 = \(\(y(0)\)
	a\_1 = 
        \(\(-\(a\_0\/4\)\)
	a\_n = 
          \(\((n - 2)\)\ a\_\(n - 2\) - a\_\(n - 1\)\)\/\(4\ n\)\)\)\)\)

But that's another story...(see section "Building power series from recursion relations")

4. Solving (4+x^2)y'[x]+y[x]==0

The same thing we did with the right hand side of the equation can be carried on with the whole equation: being essentialy rule-based, Summa manipulations procedures get mapped inside the Equal structure. Let's start with the equation

eq=(4+x^2)y'[x]+y[x]==0

As usual we build the power series expansions around our ordinary point (x=0). We choose to explicitate the first zero coefficients of the derivatives.

\(\(\ Clear[ys, y1s, y2s, a]
{
          ys = ∑\+\(n = 0\)\%∞\( a\_n\) x\^n, 
          
y1s = SumSplit[SumD[ys, x], 1], 
          
y2s = SumSplit[SumD[y1s, x], 2]} // ColumnForm
      \)\)

Substitution into the original equation yields

eqs=eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}

As before we expand and bring the factors in

eqs=BringIn[SumExpand[eqs]]

then we align exponents and starting point by reindexing and element extraction

eqs=SumExponentAlign[eqs,x,n->n]
SumSplit[eqs,1]

Simplification and solving are what's left to do. Here we use another method to extract the coefficient of (x^0)

res=SumSimplify[%]

Solve[SumCoefficient[res, x^n] == 0, (Subscript[a, n + 1])] /. n -> n - 1 // 
      Simplify

Solve[res /. _Summa -> 0, (Subscript[a, 1])]

This is the recursion relation that allow ys to compute every coefficient in the power series of the solution around x=0..

\(TraditionalForm\`\(	
    a\_0 = \(\(y(0)\)
	a\_1 = 
        \(\(-\(a\_0\/4\)\)
	a\_n = 
          \(\((n - 2)\)\ a\_\(n - 2\) + a\_\(n - 1\)\)\/\(4\ n\)\)\)\)\)

5. Solving (x+1)y''[x]+y'[x]+x y[x]==0

eq=(x+1)y''[x]+y'[x]+x y[x]==0

We are looking for a solution in the form of power series. Fuchs' theorem assures us that such a solution exist as long as we stay clear from the point x=-1, that makes the coefficient of the highest derivative zero. The function and its derivatives are

{ys = Sum[(Subscript[a, n]) x^n, {n, 0, Infinity}], y1s = DSum[ys, x], 
        y2s = DSum[y1s, x]} // ColumnForm

We choose to ignore the fact that the derivatives can be further simplified by acounting for their first null terms. Substituting these espressions in the original equation we get:

eqs=eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}

Let's expand the fist term, bringing all the factors inside the sums

eqs=SumExpand[BringIn[eqs]]

We're trying to ge all the sums with terms having the same form for the power of x, namely (Subscript[c, n]) (x^n): to achieve this goal we call SumExponentAlign, in this way (the third argument is in the form (name of the index)->(form into which it has to be cast) )

eqs=SumExponentAlign[eqs,x,n->n]

Then we try to get sums that have the same starting value, without tampering with the indexes any more: we use SumSplit to get the exceeding terms out of the sums. This command will make all the sums starting from 1

SumSplit[eqs,1]

Let's simplify this expression, compacting all the sum into a single one

res=SumSimplify[%]

Now all we have to do is to apply the principle of identity of polynomials. This can be done as follows: 
This suppresses the Sum[ , it] structure giving us an equation than can be used to find the relationship between the coefficients of the power series of the solution

eqRicors=GetElement[res]

And then we can equate the coefficients of (x^0) and (x^n) to zero, solving for the unknown coefficients (Subscript[a, n]). A shift in the index and a further simplification will yield a neat expression for the general coefficient with n>2.

{c1} = Solve[Coefficient[eqRicors[[1]], x, 0] == 0, (Subscript[a, 2])]

{cn} = 
      Solve[Coefficient[eqRicors[[1]], x, n] == 0, (Subscript[a, n + 2])] /. 
          n -> n - 2 // Simplify

We have no formula for the first two coefficients a0 and a1, and that is all right since for a second order ODE the general solution happens to depend upon two parameters. In this case it easy to see, by inspecting the power series for y[x] and y'[x], that (Subscript[a, 0]) = y(0) and (Subscript[a, 1]) = y' (0). Every coefficient of the power series can then by computed according to the following relationships

(Subscript[a, 0]) = y(0)
	(Subscript[a, 1])=y'(0)
	(Subscript[a, 2]) = -((Subscript[a, 1])/2)
	(Subscript[a, n]) = 
            -(((Subscript[a, n - 1]) (((n - 1))^2) + 
                    (Subscript[a, n - 3]))/((n - 1) n))

Building power series from recursion relations

We have applied the Summa procedures to reach a symbolic recursive form for the coefficients of the power series solution of our ODEs. We might also want to make this solution a little bit less ethereal and getting down to numbers. While there are more efficient techiniques to do this right away by starting  form the equation, the following shows how to avoid recursion errors in handling the coefficient expression.
Let's focus on the last equation's solution:

(x + 1) y''(x)+y'(x)+x y(x)=0

The recursion relation for the coefficients of the power series for y(x) is

(Subscript[a, 2]) = -((Subscript[a, 1])/2)
	\(TraditionalForm
      \`a\_n = \(-
          \(\(a\_\(n - 1\)\ \((n - 1)\)\^2 + a\_\(n - 3\)\)\/\(\((n - 1)\)\ n
                \)\)\)\)

We have to specify the value of the first two coefficients. We can do that by imposing initial conditions. Since

y(0) = (Subscript[a, 0])
	y'(0)=(Subscript[a, 1])

Imposing the I.C.y[0=]=1 and y'[0]=0 will result in choosing the first two coefficients as (Subscript[a, 0])= 1 and (Subscript[a, 1])= 0 We have enough information to compute the whole sequence of coefficients

al[0] = 1; \(
\)al[1] = 0; \(
\)al[2] = -al[1]/2; \(
\)al[n_  /; n > 2] := 
      al[n] = 
        -((al[n - 1] (((n - 1))^2) + al[n - 3])/((n - 1) n
                ))

Please note that the general recursion relation stores the computed values so that it won't be necessary to repeat their calculation for every new coefficient. Here are the first ten coefficients

Table[al[k], {k, 0, 9}]

The truncated series is easily computed as a scalar product with the vector of integer powers of x.

f[x_] = % . Table[x^k, {k, 0, 9}]

And we can see that it does indeed satisfy the given initial conditions (passing by the point (0,1) with zero slope). Obviously, this being a series approximation of the actual exact solution, it will be a good substitute for it only in a reasonably small neighbor of x=0. Also remember that the power series solution we have found is not good in a neighbor of the point x=-1
Let's compare it to the exact solution:

g[x_] = 
      y[x] /. DSolve[{
                (x + 1) y''[x] + y'[x] + x y[x] == 0, 
                
y[0] == 1, y'[0] == 0}, y[x], x][
            [1]] // EnableSumEvaluation

Block[{$DisplayFunction = Identity}, 
      
exactSol = 
        Plot[g[x], {x, -1, 1}, Frame -> True]; 
      seriesSol = 
        Plot[f[x], {x, -1, 1}, PlotStyle -> Dashing[{ .01,  .02}]]; 
      
]; 
Show[exactSol, seriesSol];

Plotting the error with full range, will show that x=-1 is indeed problematic:

Plot[g[x] - f[x], {x, -1, 1}, Frame -> True, 
      PlotRange -> All];

Building a truncated power series right from the equation (*)

If we are seeking for a finite truncaton of the infinite series expansion of the solution, we could apply a different technique, far more straightforward and efficient. What follows has been adapted from posts in the MathGroup.

The long way (*)

Consider the differential equation

(t  +  1)*y''[t]  +  y'[ t]  +  t y[t] == 
      0

This is the differential operator associated to that equation

T[y_]  :=  
        (t  +  1)*D[y,  t,  t]  +  D[y,  t]  +  t*y;

We look for a solution whose taylor expansion is

\(TraditionalForm\`\(	y(t) = ∑\+\(k = 0\)\%n a\_k\ t\^k\)\)

If we content ourselves with a maximum of 9 terms, we have

n  =  9;
y = Sum[a[k] t^k, {k, 0, 9}]

If we plug this expression in the differential operator we get

terms  =  Simplify[T[y]]

By equating this to zero and by aspplying the principle of identity of polinomyals we can get a system of equations giving the value of the coefficients

Tycoeffs  =  
        CoefficientList[terms,  t][[Range[n  -  1]]]

(eqns  =  Map[(#  ==  0) &,  Tycoeffs])//ColumnForm

Now solve for the coefficients in the solution satisfying y[0]=1, y'[0]=0. The initial conditions give us two more equations in the coefficients of the series expansion

{y == 1, D[y, t] == 0} /. t -> 0

We append these two further equations to those we already had, to have a system of 10 equations in 10 unknowns, whose solution is

(ycoeffs1 = First[Solve[Join[%, eqns], Table[a[i], {i, 0, n}]]]) // 
      ColumnForm

We can substitute this value in the expression of y to get the series expansion of our solution up to the ninth power of t:

y1[t_] = y /. ycoeffs1

Now we solve the differential equation with the initial conditions y[0]=0, y'[0]=1.This imposes the following condition on the coefficients of the series expansion

{y == 0, D[y, t] == 1} /. t -> 0

(ycoeffs2 = First[Solve[Join[%, eqns], Table[a[i], {i, 0, n}]]]) // 
      ColumnForm

y2[t_]=y /.  ycoeffs2

seriesSols = 
      Plot[{y1[t], y2[t]}, {t, -1, 1.5}, Frame -> True, 
        PlotStyle -> {{Dashing[{ .02,  .01}], Hue[ .66]}}];

The solution given by Mathematica is

\(Clear[y]; 

    y1e[t_] = 
      y[t] /. First[
          DSolve[\ {
              
\((t\  + \ 1)\)*\(\(y'\)'\)[t]\  + \ 
                  \(y'\)[\ t]\  + \ t\ y[t] == 0, 
              
y[0] == 0, 
              \(y'\)[0] == 1
}, y[t], t]]\)

\(Clear[y]; 

    y2e[t_] = 
      y[t] /. First[
          DSolve[\ {
              
\((t\  + \ 1)\)*\(\(y'\)'\)[t]\  + \ 
                  \(y'\)[\ t]\  + \ t\ y[t] == 0, 
              
y[0] == 1, 
              \(y'\)[0] == 0
}, y[t], t]]\)

exactSols = 
      Plot[{y1e[t], y2e[t]}, {t, -1, 1.5}, Frame -> True];

These compares the exact and the approximate solutions. Please note that for t<1 one of the exact solutions is not even real.

Show[seriesSols, exactSols]

The short way (*)

By using Series. To be added in a next release of this notebook.

Frobenius method: regular singular points

0. The Indicial Equation

Here's a generical second order equation with variable coefficients.

eq=y''[x]+p[x]y'[x]+q[x]y[x]==0

Let's expand the coefficients p[x] and q[x] into power series (we use x0=0 to make things simpler)

{qs=Sum[cq[n] x^(n-2),{n,0,Infinity}],
 ps=Sum[cp[n] x^(n-1),{n,0,Infinity}]}//ColumnForm

We do the same for the function y[x] considering a modified series to cope with the eventual singularity

ys = Sum[(Subscript[a, n]) x^(r + n), {n, 0, Infinity}]

Now we can substitute this tentative solution (along with its derivatives) into the original equation, obtaining

eqs=eq/.{
	p[x]->ps,q[x]->qs,
	y[x]->ys,y'[x]->SumD[ys,x],y''[x]->SumD[ys,{x,2}]
	}

This is a rather complicated expression, but let's see what the first term looks like. We split the sums at position 1 in order to make the term for n=0 explicit. By expanding the result and further simplifying the result we get the following:

eqs2=SumSimplify[SumSplit[eqs,1]//SumExpand//BringIn]

A necessary condition for the above equality to be true is that the coefficient of (x^(r - 2)) be zero. This results in the celebrated indicial equation, that poses a condition on the value of r to be used in the potential solution \(TraditionalForm
      \`y[x] = x\^r\ Sum[\(a\_k\) x\^k, {k, 1, ∞}]\)

indicialEq=Coefficient[#,x^(r-2)]& /@ eqs2

The indicial equation solutions are functions of the first coefficients of the power series expansions of p[x] and q[x]

(sol=Solve[indicialEq,r])//ColumnForm

The form of the solution  depends on whether r1 and r2 are distinct, differing by an integer or equal.

1. Distinct solutions 			(r1-r2 not integer)

This differential equation, in a neighborhood of x=0

eq=x y''[x]+(1+x)y'[x]-1/(16x)y[x]==0;

has an indicial equation with solutioni r=1/4 ed r=-1/4. Since they do not differ by an integer, Frobenius method states that  we can find two linear independent solutions in the form

\(TraditionalForm\`\(	
    \(y\_1\)(x) = 
      \(\@x\%4\ \(∑\+\(n = 0\)\%∞ a\_n\ x\^n\)
	
          \(\(y\_2\)(x)\) = 
        \(1\/\@x\%4\) \(∑\+\(n = 0\)\%∞ b\_n\ x\^n\)\)\)\)

Let's find the first one, for r=1/4

\({ys = x^\((1/4)\) Sum[\(a\_n\) x^n, {n, 0, Infinity}], 	
\ 
        y1s = SumD[ys, x], 
\ y2s = SumD[y1s, x]} // ColumnForm\)

The equation has already polynomial coefficients so we don't need to expand them in power series. We plug the power series for our tentative solution and its derivatives into the original equation to get

eqs=eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}

eqs=SumExpand[BringIn[eqs]]

We align the exponents to n-3/4

eqs2=SumExponentAlign[eqs,x,n->n-3/4]

We make all the sum start from 1

eqs2=SumSplit[eqs2,1]

And we squeeze all the sums intoa single one

res=SumSimplify[eqs2]

Then we solve for coefficents

sol = Solve[SumToTerm[res], (Subscript[a, n])][[1]] // Simplify

The second linear independent solution is found in the same way

\({ys = x^\((\(-1\)/4)\) Sum[\(b\_n\) x^n, {n, 0, Infinity}], 	
\ 
        y1s = SumD[ys, x], 
\ y2s = SumD[y1s, x]} // ColumnForm\)

We expand and bring inside the coefficients...

eqs=SumExpand[
	BringIn[
		eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}]
]

Add some manipulation...

eqs2=SumExponentAlign[eqs,x,n->n-5/4]

res=SumSplit[eqs2,1]//SumSimplify

and solve...

sol = Solve[SumToTerm[res], (Subscript[b, n])][[1]]

2. Integer separated solutions 	(r1-r2 integer)

In a neighborhood of 0, the equation

eq=y''[x]+(1-2/(3x))y'[x]+(4/(9x^2)+x)y[x]==0;

has indicial equation with solutions 1/3 and 4/3, hence differing by an integer. The Frobenius method states that we must look for solutions in the form

\(TraditionalForm\`\(	
    \(y\_1\)(x) = 
      \(x\^\(4/3\)\ \(∑\+\(n = 0\)\%∞ a\_n\ x\^n\)
	
          \(\(y\_2\)(x)\) = 
        \(x\^\(1/3\)\) \(∑\+\(n = 0\)\%∞ b\_n\ x\^n\) + 
          c\ \(\(y\_1\)(x)\)\ \(log(x)\)\)\)\)

The first one, for r=4/3

\(\(\ {ys = x^\((4/3)\) Sum[\(a\_n\) x^n, {n, 0, Infinity}], 	
\ 
        y1s = SumD[ys, x], 
\ y2s = SumD[y1s, x]} // ColumnForm\)\)

Yields

eqs=SumExpandAll[
	BringIn[		
		eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}]
	]

(Note that we used ExpandAllSum : with ExpandSum only we would have needed a further expansion to get the same result). Now we can align all the exponents to n+7/3

eqs2=SumExponentAlign[eqs,x,n->n+7/3]

Breaking in 0 brings us to

res=SumSplit[eqs2,0]//SumSimplify

As usual we solve for coefficients. Note that we use GetCoefficient to extract the coefficient of x^(n+7/3) from the sum. And we consider res[[1]] to get the rhs of res.

{eq1=Coefficient[res[[1]],x^(1/3)]==0,
eq2=Coefficient[res[[1]],x^(4/3)]==0,
eqn=SumCoefficient[res[[1]],x^(n+7/3)]==0}//ColumnForm

sol = 
      Solve[{eq1, eq2, eqn}, {(Subscript[a, n + 3]), (Subscript[a, 2]), (Subscript[a, 1])}][[1]]

The second linear indepent solution may contain a Log term

\(TraditionalForm\`\(	
    \(y\_2\)(x) = 
      \(\@x\%3\) \(∑\+\(n = 0\)\%∞ b\_n\ x\^n\) + 
        k\ \(\(y\_1\)(x)\)\ log\ \((x)\)\)\)

We can use  a  brute force approach to insert this solution and its derivatives into the equation. Care must be used to avoid cumbersome expressions

\(ys = 
      k\ Log[x] x^\((4/3)\) Sum[\(a\_n\) x^n, {n, 0, Infinity}] + 
\ \ \ \ 
        x^\((1/3)\) Sum[\(b\_n\) x^n, {n, 0, Infinity}]\)

y1s=SumD[ys,x]//BringOut

y2s=SumD[y1s,x]//BringOut

Here's the equation. Not a very beautiful sight...

eqs=eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}//BringOut

The k Log[x] is quite annoying. Let's see if there is a way to get rid of it

Collect[eqs[[1]], Log[x]]

by inspection we can see that one of the terms is in the form k Log[x]*(eq /. y[x]->y1). Since y1 is a solution of the equation, this is equivalent to k Log[x]*(0)=0. Hence our equation is:

eqs=eqs/.Log[x] _->0

We now proceed as in the usual case

eqs2=SumExpand[eqs]//BringIn

eqs2=SumExponentAlign[eqs2,x,n->n+4/3]

res=SumSplit[eqs2,0]//SumSimplify

Then we solve for the coefficients

{eq1=Coefficient[res[[1]],x^(-2/3)]==0,
eq2=Coefficient[res[[1]],x^(1/3)]==0,
eqn=SumCoefficient[res[[1]],x^(n+4/3)]==0}//ColumnForm

sol = 
      Solve[{eq1, eq2, eqn}, {(Subscript[b, 0]), (Subscript[b, 2]), (Subscript[b, n + 3])}][[1]]

sol/.k->0//Simplify

sol/.k->1//FullSimplify

3. Identical solutions 		(r1-r2=0) (*)

Finally, for

eq=x y''[x]+(2-x)y'[x]+1/(4x)y[x];

we have r=-1/2 and r=0. The solutions are

\(TraditionalForm\`\(	
    \(y\_1\)(x) = 
      \(\(1\/\@x\) 
          \(∑\+\(n = 0\)\%∞\(\( a(\ n)\)\ x\^n\)
	
              \(\(y\_2\)(x)\)\) = 
        ∑\+\(n = 0\)\%∞ b\ \((n)\)\ x\^n + 
          \(\(y\_1\)(x)\)\ \(log(x)\)\)\)\)

The first one is easy

{ys=(x^(-1/2)Sum[a[n]x^n,{n,0,Infinity}]),	
 y1s=SumD[ys,x],
 y2s=SumD[y1s,x]}//ColumnForm

we expand with respect to p and q and bring everything inside a single sum

eqs=SumExpand[
	BringIn[
		eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}]
]

Further manipulation brings to

eqs2=SumExponentAlign[eqs,x,n->n-3/2]

eqs2=SumSplit[eqs2,1]

res=SumSimplify[eqs2]

and solving

sol=Solve[GetCoefficient[res,x^(n-3/2)]==0,a[n]][[1]]

The second independent solution is less easily found, but can be obtained in the same way as in case 2. As is customary, being this a somewhat tedious task,  let's say that it is left as an exercise.

Hypergeometric Equation (*)

This is Gauss' hypergeometric equation. Needless to say, its solutions include the celebrated hypergeometric functions.

gaussEq =  
      x (1 - x) y''[x] + (c - (a + b + 1) x) y'[x] - 
          a b y[x] == 0

In order to solve it with the Frobenius method, we need to find what values of r the tentative power series solution

ys = x^r Sum[(Subscript[a, n]) x^n, {n, 0, ∞}];

can satisfy the equation. We can find it out by substituting ys and its derivative into the equations and by imposing that the coefficient with the lowest power of x be zero.

eqs = 
      gaussEq /. {y[x] -> ys, y'[x] -> SumD[ys, x], 
              y''[x] -> SumD[ys, {x, 2}]} // SumExpandAll // 
        SumFullSimplify

SumSplit[eqs, 1] /. Sum[_, _] -> 0

Collect[# , x]& /@%

The smallest power is (x^(r - 1)). By equating its coefficient to zero we get the indicial equation, whose solutions are r = 0 and r = 1 - c.  We can now try to find a power series solution in either one of the two forms

\(TraditionalForm\`\(	
    ∑\+\(n = 0\)\%∞ a\_n\ x\^n
	x\^\(1 - c\)\ 
        \(∑\+\(n = 0\)\%∞ a\_n\ x\^n\)\)\)

Euler-Cauchy Equation

A little 'amarcord'. Here we solve the Cauchy-Euler equation  by enabling Sum evaluation and using old naming conventions for the procedures defined in Summa.m. We will then work with Summa objects directl as it was done with version 1.x of the Summa.m package.

EnableSumEvaluation[]

Clear[y,x,a];
cauchyEq=x^2 y''[x]-2x y'[x]-10y[x]==0

We can get the indicial equation directly by manipulating the equation to be solved. We look for a solution in the form of a power series about the point (Subscript[x, 0]) = 0. We do not know whether such a solution exists or not, since the theorem of existence is not applicable.

\(\(\ {ys = x^r\ Summa[\(a\_n\) x^n, {n, 0, Infinity}], 
\ 
        y1s = DSum[ys, x], 
\ y2s = DSum[y1s, x]} // ColumnForm\)\)

with r that needs to be determined. Substitution into the equation yields

eqs=cauchyEq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}

Now we begin the usual manipulation of sums in order to reach a form with a single sum, plus some extra terms. We start by clumping all the Sums together

eqs=SimplifySum[eqs]

the first element is in the form c (x^r), as we can see by expliciting every term with index n<1.

eqs2=BreakSum[eqs,1]

hence, for the identity principle, we need r to satisfy the following indicial equation

indicialEq=Coefficient[#,x^r]& /@ eqs2

whose solutions are easily found:

Solve[indicialEq,r]

Solving for r = -2

eqs=ExpandSum[BringIn[
	cauchyEq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}/.r->-2]
]

SimplifySum[eqs]

One obvious solution is y(x)=0, corresponding to a series with the general term (Subscript[a, n]) = 0. A non trivial solution can be found when n=0 and n=7 regardless of the values of (Subscript[a, n]). We can then find a particular solution by specifying two non zero values for the coefficients (Subscript[a, 0]) and (Subscript[a, 7])

y1[x_] = (Subscript[a, 0]) x^(0 - 2) + (Subscript[a, 7]) x^(7 - 2)

Let's find a second solution by

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

Hence our general solution has no logarithm and is

Collect[c1 y1[x] + c2 y2[x], x]

Let's simplify things a bit, by renaming the constants

yg[x_]=C[1] x^5+C[2] x^(-2)

We can verify that this is indeed the general solution given by Mathematica. Note that, being the built-in Sum evaluation process already restored, there is no need to wrap EnableSumEvaluation around DSolve.

DSolve[x^2 y''[x] - 2  x y'[x] - 10  y[x] == 0, 
      y[x], x]

Other examples: Special Functions

Introduction to Sturm Liouville problems and special functions (*)

Many functions of mathematical and physical interest are solution of a particular class of differential equations subjected to specific conditions. The following second order differential equation is the starting point for many such systems. In the following α(x), β(x) and γ(x) are real functions of the real variable x.

α(x) y''(x)+β(x) y'(x)+ γ(x) y(x)=λ y(x)

This equations is better seen as an eigenvalue problem. We introduce a differential operator 𝕃 that, when applied to y(x) produce the left hand side of the above equation

\(TraditionalForm\`\(	
    \([\(α(x)\)\ d\^2\/dx\^2 + 
            \(β(x)\)\ d\/dx + \ 
            γ(x)]\) \(y(x)\) = λ\ \(y(x)\)\)\)

To reach a more compact representation

\(TraditionalForm\`\(	𝕃\ y\  = λ\ y\)\)

The differential operator 𝕃 is an hermitian operator if

\(TraditionalForm\`\(	
    ⟨𝕃\ ϕ, 
          ψ⟩\_w = 
      ⟨ϕ, 
          \(𝕃\&_\) ψ⟩\_w\)\)

for every pair of functions ϕ and ψ. It can be shown (see for example Byron and Fuller) that this is the case in every Hilbert space with an inner product whose weigth function w(x) satisfies

(d/d)x[
              w(x) α(x)] = w(x) α(x)
\(	
            Lim\+\(x -> a\) \([
                  \(\( w(x)\)\ \(α(x)\)\ 
                    \((\(ϕ(x)\)\&_\ ψ' \((x)\) - 
                        \(ϕ' \((x)\)\)\&_\ \(ψ(x)\))\)\) ]\) = 
              Lim\+\(x -> b\) \([
                  \(\( w(x)\)\ \(α(x)\)\ 
                    \((\(ϕ(x)\)\&_\ ψ' \((x)\) - 
                        \(ϕ' \((x)\)\)\&_\ \(ψ(x)\))\)\) ]\)\)

for every pair of functions ϕ and ψ in the Hilbert space. (Here the overbar stands for complex conjugation, while a and b are the endpoints of the interval over which the functions are defined). The eigenvalue problem admits then real eigenvalues and the solutions form an orthonormal basis spanning the Hilbert space charaterized by the weigth function w(x) in its inner product  \(TraditionalForm
      \`⟨·\(, 
            ·\)⟩\_w\).

To make explicit the role of the weigth function w(x) in the determination of the solution of the eigenvalue problem, we can recast the problem in the following form:

(d/d)x[
              w(x) α(x) 
                d((y(x))/d)x] + 
            (γ(x) - λ) w(x) y(x) = 0

This is a slightly more compact form, often found in textbooks:

\(TraditionalForm\`\(	
    \(d\/dx\)[
          \(p(x)\) d\(y(x)\)\/dx] + 
        \(q(x)\)\ \(y(x)\) = λ\ \(w(x)\)\ \(y(x)\)\)\)

The above equation is named after Sturm and Liouville and, along with the specified boundary conditions (which allow the selection of the orthogonal solutions only), constitute a Sturm-Lioville problem. A great deal of special functions used in mathematical physics can be derived by this kind of problem with specific instances of the coefficients α(x), β(x) and γ(x) and of the weigth function w(x).

If we limit ourselves to the search for polynomial solutions to the Sturm-Liouville problem, the following exhausts the alternatives

({{Interval, Polynomial,  , w(x),  , Differential Equation}, { ,  ,  ,  ,  ,  }, {[0, 1], Jacobi,  , (((x^p)(1 - x))^q),  , x(1 - x) y'' - ((p + q + 2) x + p + 1) y' + 
                  n(n + p + q + 1) y = 0}, {[-1, 1], Jacobi,  , (((1 + x))^p) (((1 - x))^q),  , (1 - (x^2)) y'' - 
                  ((p + q + 2) x + q - p) y' + 
                  n(n + p + q + 1) y = 0}, { ,  , Gegenbauer, (((1 - (x^2)))^p),  , (1 - (x^2)) y'' - 2 (p + 1) x y' + 
                  n(n + 2  p + 1) y = 0}, { ,  , Chebychev, (((1 - (x^2)))^(-1/2)),  , (1 - (x^2)) y'' - x y' + (n^2) y = 0}, { ,  , Legendre, 1,  , (1 - (x^2)) y'' - 2  x y' + n(n + 1) y = 0}, {[0, ∞ ), Laguerre,  , (x^p) (e^(-x)),  , x y'' - (x - (p + 1)) y' + n y = 0}, {(-∞, ∞), Hermite,  , (e^(-(x^2))),  , y'' - 2  x y' + 2  n y = 0}})

There are different kinds of special functions that are not polynomials. The ordinary trigonometric functions arising from the differential equation of the harmonic motion are just the simplest example

Hermite Equation

hermiteEq = y''[x] - 2  x y'[x] + 2  n y[x] == 0;

Here we set up a list of replacement rules that exploit the fact that the derivative of a power series starting from zero loses its constant terms

\(ys = Sum[\(a\_k\) x\^k, {k, 0, ∞}]; 

    rules = {y[x] -> ys, 
        
\(y'\)[x] -> SumSplit[SumD[ys, x], 1], 
        
\(\(y'\)'\)[x] -> 
          SumSplit[SumD[ys, {x, 2}], 2]}\)

Plugging this into the orignal equation yields

hermiteEq /. rules // BringIn

As usual we handle the expression with the aim to reduce all sums to a single sum equated to zero. Shifting indices so that every sum will feature the same power of x is the first step.

eqs = SumExponentAlign[%, x, k -> k] // SumSimplifyEach

There are two way to make every  sum start from the same index. Either we explicitate the terms in excess (that happens to be in the second and third sums)...

SumSplit[%, 1] // SumSimplify

(but this approach has a shortcoming: it gives an additional equation for the coefficient of x^0 - namely (Subscript[a, 2]) = -n (Subscript[a, 0])- that is not strictly necessary to our purpose)
...or we consider the fact that the first sum has a term proprtional to k that goes to zero when k=0. It is then immaterial whether its summation index start from zero or from one.

(∑_(k = 1)^∞)k (x^k) (Subscript[a, k])=0+(∑_(k = 1)^∞) k (x^k) (Subscript[a, k])=(∑_(k = 0)^∞) k (x^k) (Subscript[a, k])

We simply have to apply a transformation to the first sum to change the value of the starting value without touching anythig else. We can use MapAt in conjuction with SumPosition to identify the position of the first sum inside the equation.

MapAt[
        (# /. Sum[el_, {k, 1, end_}] :> Sum[el, {k, 0, end}])
            &, eqs, SumPosition[eqs, 1]] // SumSimplify

A power series is identically zero for every x iff its coefficients are all  zero. We can obtain the equation stating the condition on the coefficients by simply suppressing the Sum structure. Solving for the highest coefficient will give us a recursion relation

eqr = SumToTerm[%]\(
\)Solve[Coefficient[eqr[[1]], x, k] == 0, (Subscript[a, k + 2])][[1]]

We have all we need to compute a series solution of Hermite's equation.

(Subscript[a, 0]) = y(0)
	(Subscript[a, 1])=y'(0)
	(Subscript[a, k + 2]) -> 
            ((2 (k - n) )/((k + 1) (k + 2))) 
              (Subscript[a, k])

The presence of the (k-n) factor has the effect of turning the n-th coefficient to zero, no matter the initial conditions chosen. Since the generic coefficient (Subscript[a, k + 2])depends upon the k-th computed coefficient this will turn to zero every next coefficient. If n is odd, this means that odds coefficients beyond (Subscript[a, n])will be zero. If n is even, then all even coefficients beyond (Subscript[a, n])wll be zero. At least one solution of the Hermite equation, for integer n, is a polynomial of degree n, then. We can make all other coefficients zero by setting the initial conditions (the values of (Subscript[a, 0])and (Subscript[a, 1])).

DSolve[{y''[x] - 2  x y'[x] + 2 n y[x] == 0}, y[x], 
        x] // EnableSumEvaluation

Hermite polynomial for even n

\(TraditionalForm\`\(	
    a\_0 = \(1
	a\_1 = 
        0
	a\_\(k + 2\) -> 
          \(\(2\ \((k - n)\)\ \)\/\(\((k + 1)\)\ \((k + 2)\)\)\) a\_k\)\)\)

\(ae[0] = 1; 
ae[1] = 0; 

    ae[k_\  /; k ≥ 2] := 
      \(ae[k] = 
        \(\(2\ \((k - 2 - n)\)\ \)\/\(k \((k - 1)\)\)\) ae[k - 2]\)\)

Here are the first ten coefficients

evenCoeff = Table[ae[k], {k, 0, 10}]

(evenCoeff /. n -> #)&/@{0, 2, 4, 6, 8, 10}

Our series is

(hermiteEven = % . Table[x^k, {k, 0, 10}]) // ColumnForm

We can see that the polynomials obey the initial conditions y(0)=1, y'(0)=0, i.e. they pass with zero slope by the point (0,1)

Plot[Evaluate[hermiteEven], {x, -1, 1}, Frame -> True];

These polynomials differs from those given by Mathematica only by a (signed) normalization factor.

Table[Simplify[hermiteEven[[n + 1]]/HermiteH[2  n, x]], {n, 0, 5}]

Hermite polynomial for odd n

\(TraditionalForm\`\(	
    a\_0 = \(0
	a\_1 = 
        1
	a\_\(k + 2\) -> 
          \(\(2\ \((k - n)\)\ \)\/\(\((k + 1)\)\ \((k + 2)\)\)\) a\_k\)\)\)

\(ao[0] = 0; 
ao[1] = 1; 

    ao[k_\  /; k ≥ 2] := 
      \(ao[k] = 
        \(\(2\ \((k - 2 - n)\)\ \)\/\(k \((k - 1)\)\)\) ao[k - 2]\)\)

Here are the first ten coefficients

oddCoeff = Table[ao[k], {k, 0, 9}]

(oddCoeff /. n -> #)&/@{1, 3, 5, 7, 9}

Our terinating series are

(hermiteOdd = % . Table[x^k, {k, 0, 9}]) // ColumnForm

Se the effect of the different I.C. on these solutions: they all pass by the origin with slope 1.

Plot[Evaluate[hermiteOdd], {x, -1, 1}, Frame -> True];

These polynomials differ from those given by Mathematica only by a normalization factor. The values are the same seen in the case for even n.

Table[
      Simplify[hermiteOdd[[n + 1]]/HermiteH[2  n + 1, x]], {n, 0, 4}]

Legendre Equation (*)

legendreEq = 
      (1 - (x^2)) y''[x] - 2  x y'[x] + n (n + 1) y[x] 
        == 0;

The coefficient of the highest derivative is (1 - (x^2))and is not zero in x=0. We can then look for a power series solution with center in x=0.

eqs = 
      legendreEq /. {
          y[x] -> (∑_(k = 0)^∞) (Subscript[a, k]) (x^k), 
          
y'[x] -> 
            SumSplit[SumD[(∑_(k = 0)^∞) (Subscript[a, k]) (x^k), x], 
              1], 
y''[x] -> 
            SumSplit[
              SumD[(∑_(k = 0)^∞) (Subscript[a, k]) (x^k), {x, 2}], 2]}

Expanding every sum will show sums with different endpoints and different values for the exponent of the general term.

SumExpandAll[eqs] // BringIn

Shifting the indices to make all the power of x alike, will yield the following equation, that can be further simplified by splitting the sums at the highest startpoint, i.e. k=2.

SumExponentAlign[%, x, k -> k]\(
\)eq = SumSplit[%, 2] // SumSimplify

We now have an equation that can be easily solved: applying the principle of identity of polynomials will result in three equations

eqr = SumToTerm[eq]; \(
\)({eq1, eq2, eq3} = {SumCoefficient[eqr[[1]], x, 0] == 0, 
            
SumCoefficient[eqr[[1]], x, 1] == 0, 
            
SumCoefficient[eqr[[1]], x, k] == 0})
         // ColumnForm

Solving w.r.t. the highest coefficient in each equation yields the recursion relation needed to compute every coefficient in the power series solution. We factorize the results to achieve a neater form

(Flatten[
            
{Solve[eq1, (Subscript[a, 2])], Solve[eq2, (Subscript[a, 3])], 
              Solve[eq3, (Subscript[a, k + 2])]}] /. 
          Rule[a_, b_] :> Rule[a, Factor[b]]) // ColumnForm

Even coefficients happen to depend only on the initial value of y, while odd coefficients depends on the initial value of the first derivative y'. We can identify two sequences of coefficients

({{(Subscript[a, 0]) = (Subscript[y, 0]),  , (Subscript[a, 1]) = (Subscript[yp, 0])}, {(Subscript[a, 2]) = -(( n (n + 1))/(2!)) (Subscript[y, 0]),  , (Subscript[a, 3]) = 
                -(((n - 1) (n + 2))/(3!))  (Subscript[yp, 0])}, {(Subscript[a, 4]) = 
                -(((n - 2) (Subscript[a, 2]) (n + 3))/12) = 
                  (((n - 2) n (n + 1) (n + 3))/(4!))  
                    (Subscript[y, 0]),  , (Subscript[a, 5]) = 
                (( (n - 3) (Subscript[a, 3]) (n + 4))/20) = 
                  (( (n - 3) (n - 1) (n + 2) 
                          (n + 4))/(5!)) (Subscript[yp, 0])}, {(Subscript[a, 6]) = 
                -(((n - 4) (Subscript[a, 4]) (n + 5))/30) = 
                  -(((n - 4) (n - 2) n (n + 1) 
                            (n + 3) (n + 5))/(6!))  
                    (Subscript[y, 0]),  , (Subscript[a, 7]) = 
                -(( (n - 5) (Subscript[a, 5]) (n + 6))/42) = 
                  (((n - 5) (n - 3) (n - 1) (n + 2) 
                          (n + 4) (n + 6) )/(7!)) (Subscript[yp, 0])}, {…,  , …}})

An interesting thing to note is that, since each sequence depends upon the previous coefficient, and since for certain integer values of n it is possible to set to zero one of the coefficients, then our series will terminate in a polynomial. In particular, for even n we will have a finite sequence of even coefficients (and an infinite series with odd coefficients), while for odd n we will find a polynomial in odd powers of x (and an infinite series in even powers of x). If we write the generic even and odd coefficients as

\(TraditionalForm
      \`a\_\(2  m\) = 
        \((\(-1\))\)\^m\ 
          \(\(\(\((n - 2  m + 2)\)\ \((n - 2  m + 4)\) ... \)\ n\ 
                  \((n + 1)\) ... \)\ \((n + 2  m - 3)\) 
              \((n + 2  m - 1)\)\)\/\(\((2  m)\)!\)\ \ y\_0\)
	\(TraditionalForm
      \`a\_\(2  m + 1\) = 
        \((\(-1\))\)\^m\ 
          \(\(\(\((n - 2  m + 1)\)\ \((n - 2  m + 3)\) ... \)\ \((n - 1)\)\ 
                  \((n + 2)\) ... \)\ \((n + 2  m - 2)\) 
              \((n + 2  m)\)\)\/\(\((2  m + 1)\)!\)\ \ yp\_0\)

Then we can express our solution in the following way

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

For even n we obtain an even polynomial multiplied by y(0).For odd n we get an odd polynomial multiplied by y'(0)The polynomials are named Legendre polynomial (LegendreP).The other linear independent solution is known as a Legendre function of the second kind and it is customary to denote it with the letter Q (LegendreQ).

DSolve[{
          (1 - (x^2)) y''[x] - 2  x y'[x] + 
              n (n + 1) y[x] == 0}, y[x], x] // EnableSumEvaluation

Moreover, if we specify the initial conditions so that the a0 is zero for odd n and a1 is zero for even n, the particular solution obtained is the legendre polynomial of order n.

Legendre Polynomials for even n

LegendreP[0, x]LegendreP[2, x]

Series[LegendreQ[2, x], {x, 0, 12}]

Legendre Polynomials for odd n

LegendreP[1, x]LegendreP[3, x]

Series[LegendreQ[1, x], {x, 0, 12}]

Airy Equation (*)

airyEq =  y''[x] - x y[x] == 0;

{ys = (∑_(k = 0)^∞) (Subscript[a, k]) (x^k),  
        y2s = SumSplit[SumD[ys, {x, 2}], 2]} // ColumnForm

eqs = 
      airyEq /. {y[x] -> ys, y''[x] -> y2s} // BringIn

SumExponentAlign[eqs, x, k -> k] // Simplify

It's about time to gather everything under the same sum sign. We just have to split the second sum at k=1 to get

SumSplit[%, 1] // SumSimplify

We can now obtain the recursion relation. To speed things up a bit we mix the k-th coefficient with the coefficient of the zero-th power of x

eqr = SumToTerm[%]Solve[Coefficient[eqr[[1]], x, 0] == 0, (Subscript[a, 2])][[1]]Solve[Coefficient[eqr[[1]], x, k] == 0, (Subscript[a, k + 2])][[1]]

This recursion relation enable us to compute every coefficient of the power series of the solution of Airy's equation in a neighbor of x=0.

(Subscript[a, 0]) = y(0)
	(Subscript[a, 1]) = y' (0)
	(Subscript[a, 2]) = 0
	(Subscript[a, k + 2])=(1/((k + 1) (k + 2)))(Subscript[a, k - 1])

We can now build a table of the first few coefficients that will depend upon the coefficients (Subscript[a, 0])and (Subscript[a, 1])and ultimately upon the initial conditions y(0) and y'(0). Since the generic coefficient (Subscript[a, k])(for k>2) is a function of the the (k-3)th coefficient only,

% /. k -> k - 2 // Simplify

we can immediately infer that - being (Subscript[a, 2]) = 0- all the coefficients with index (2+3 n) are zero. The remaining coefficients will ultimately depen upon (Subscript[a, 1])(i.e. y'(0)) and (Subscript[a, 0])(i.e. y(0)). We are going to find two independent series solutions that will be given a name for their importance in certain physical applications.

({{(Subscript[a, 0]) = (Subscript[y, 0]),  , (Subscript[a, 1])=(Subscript[y, 0]^′),  , (Subscript[a, 2]) = 0}, {(Subscript[a, 3]) = ((Subscript[a, 0])/(3 (3 - 1))) = ((Subscript[y, 0])/(3 2)),  , (Subscript[a, 4])=((Subscript[a, 1])/(4 (4 - 1)))=((Subscript[y, 0]^′)/(4 3)),  , (Subscript[a, 5]) = ((Subscript[a, 2])/(5 (5 - 1))) = 0}, {(Subscript[a, 6]) = ((Subscript[a, 3])/(6 (6 - 1))) = ((Subscript[y, 0])/(6 5 3 2)),  , (Subscript[a, 7])=((Subscript[a, 4])/(7 (7 - 1)))=((Subscript[y, 0]^′)/(7 6 4 3)),  , (Subscript[a, 8]) = ((Subscript[a, 5])/(8 (8 - 1))) = 0}, {…,  , …,  , …}, {(Subscript[a, 3 n]) = 
                ((Subscript[y, 0])/(((3 n) (3 n - 1)) … 
                    (6 5) (3 2))),  , (Subscript[a, 3 n + 1])=((Subscript[y, 0]^′)/(((3 n + 1) ( 3 n)) … 
                    (7 6) (4 3))),  , (Subscript[a, 3 n + 2]) = 0}})

We can write our solution in the following way

y(x)=(Subscript[a, 0])+(Subscript[a, 1]) x+(∑_(n = 1)^∞)(( (Subscript[a, 0]) 
                    (x^(3  n)))/(((3 n) (3 n - 1)) 
                    … (6 5) (3 2)))+y'(0)(∑_(n = 1)^∞)(( (Subscript[a, 1]) 
                      (x^(3  n + 1)))/(((3 n + 1) (3 n)) 
                      … (7 6) (4 3)))
=y(0) 
                (1 + (∑_(n = 1)^∞) ((x
                          ^(3  n))/(((3 n) (3 n - 1)) 
                          … (6 5) (3 2))))+y'(0)(x + (∑_(n = 1)^∞) ((x
                          ^(3  n + 1))/(((3 n + 1) (3 n)) 
                          … (7 6) (4 3))))

DSolve[{
y''[x] - x y[x] == 0, 
          
y[0] == yo, 
          
y'[0] == ypo}, y[x], x] // 
      EnableSumEvaluation

Bessel Equation for noninteger μ

This is an instance of Bessel's differential equation

eq=x^2 y''[x]+x y'[x]+(x^2-μ^2)y[x]==0;

the indicial equation has solutions +μ and - μ. If the difference 2μ is not an integer (that is when μ is not an integer or an half-integer) we should look for functions in the form

\(TraditionalForm\`\(	
    \(y\_1\)(x) = 
      \(x\^μ\ \(∑\+\(n = 0\)\%∞ a\_n\ x\^n\)
	
          \(\(y\_2\)(x)\) = 
        x\^\(-μ\)\ \(∑\+\(n = 0\)\%∞ b\_n\ x\^n\)\)\)\)

The first function and its derivatives are

\(\(\ {ys = x^μ\ Sum[a\_n\ x^n, {n, 0, Infinity}], 	
	
        y1s = SumD[ys, x], y2s = SumD[y1s, x]} // ColumnForm\)\)

Let's plug them into the equation, expand all the sums and bring all factors inside the sums.

eqs=SumExpand[	
	BringIn[
		eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}]
	]

Next step: we align the exponents to the highest power

eqs2=SumExponentAlign[eqs,x,n->n+μ+2]

Now we break each sum in order to have them all starting from n=0.

eqs2=SumSplit[eqs2,0]//SumSimplify

And finally we apply the principle of identity of polynomials, solving for the coefficents. We obtain the recursion relation for the coefficient (Subscript[a, n])of the power expansion.

Solve[Coefficient[eqs2[[1]],  x^(μ + 1)]  ==  0,  
        (Subscript[a, 1])][[1]]Solve[SumCoefficient[eqs2,  x^(n + μ + 2)]  ==  0,  
        (Subscript[a, n + 2])][[1]]

If μ is not integer or half integer, the second linear indipendent solution is obtained in the same manner as before

\({ys = x^\((\(-μ\))\) Sum[\(b\_n\) x^n, {n, 0, Infinity}], 	
	
      y1s = SumD[ys, x], y2s = SumD[y1s, x]}\)

substituting in the eq

eqs=SumExpand[
	BringIn[
		eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}]
	]

Alignment to the highest power yields

eqs2=SumExponentAlign[eqs,x,n->n-μ+2]

Have them all to start from the highest starting value:

eqs2=SumSplit[eqs2,0]//SumSimplify

Solving for the coefficients yields

Solve[Coefficient[eqs2[[1]], x^(1 - μ)] == 0, (Subscript[b, 1])][
      [1]]Solve[SumCoefficient[eqs2, x^(n - μ + 2)] == 0, (Subscript[b, n + 2])]
        [[1]]

Odd coefficients are all set to zero

Bessel Equation for a particular case of integer μ (μ=2)

Let's see what happens in the case μ is integer (and hence the difference between the solutions of the indicial equation is an integer)

eq=x^2 y''[x]+x y'[x]+(x^2-2^2)y[x]==0

The first solution is in the form

\({ys = x^\((2)\) Sum[\(a\_n\) x^n, {n, 0, Infinity}], 	
	
        y1s = SumD[ys, x], y2s = SumD[ys, {x, 2}]} // ColumnForm\)

this is the equation, when we replace the above

eqs=SumExpand[	
	BringIn[
		eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}]
	]

These are the manipulation to get to the result.

eqs2=SumExponentAlign[eqs,x,n->n+4]

eqs2=SumSplit[eqs2,0]//SumSimplify

Solve[Coefficient[eqs2[[1]],  x^3]  ==  0,  (Subscript[a, 1])][[1]]
      \(
\)  sol  =  
      Solve[SumCoefficient[eqs2,  x^(n  +  4)]  ==  0,  
          (Subscript[a, n + 2])][[1]]

The second independent solution can be found by introducing a logarithmic term. (Again I want to stress the fact that this notebook is just an exemplification of the package Summa.m, and there is no pretense of mathematical rigor nor - Math forgive me- correctness. The package is just an help to avoid cumbersome calculation by hand. But it may lead to wrong results if the user does not apply the proper mathematical concepts).

\(Clear[a, y1]; 
    
{ys = x^\((\(-2\))\) Sum[\(b\_n\) x^n, {n, 0, Infinity}] + 
            k\ Log[x]\ y1[x], 
\ y1s = DSum[ys, x], 
\ 
        y2s = DSum[ys, {x, 2}]} // ColumnForm\)

eqs=BringIn[
	SumExpand[
		eq/.{y[x]->ys,y'[x]->y1s,y''[x]->y2s}]
	]

Let's try to get rid of the logarithm. We work on the left hand side of the equation in order to be able to apply Mathematica's expression manipulation functions. Collect will pave the way to a nice simplification of the above equation

Collect[%[[1]],{k Log[x],x,x^2},BringOut]

The term in parethesis that multiplies  k log(x) is just the left hand side of Bessel's equation with y(x)=y1(x). But y1(x) is a solution for that equation, and that means that the term in parenthesis is zero.We can simply neglect the logarithmic term, and we gladly isntruct Mathematica to do so:

eqs= %==0 /.Log[x] _->0

We need to substitute to y1'(x) the corresponding expression in term of power series

eqs = 
      eqs /. y1'[x] -> 
          SumD[x^(2) Sum[(Subscript[a, n]) x^n, {n, 0, Infinity}], x]

eqs2=SumExponentAlign[eqs,x,n->n+1]//SumSimplify

eqs2=SumSplit[eqs2,1]//SumSimplify

We now extract an associated equation that expresses the conditions on the coefficient of the power series representing our solution. We then use SumCoefficient to extract the coefficient of the various powers of x and then solve with respect to the

eqs3 = SumToTerm[eqs2]

Solve[SumCoefficient[eqs3[[1]], x, -1] == 0, (Subscript[b, 1])][[1]]
      Solve[SumCoefficient[eqs3[[1]], x, 0] == 0, (Subscript[b, 2])][[1]]Solve[SumCoefficient[eqs3[[1]], x, 1] == 0, (Subscript[b, 3])][[1]]Solve[SumCoefficient[eqs3[[1]], x, n + 1] == 0, (Subscript[b, n + 3])][
      [1]]

% /. n -> n - 3  /. k -> 1 // Simplify