Series Solutions of Differential Equations

Steve Dunbar
Math 221H, Fall 1993

Covers Nagle & Saff, Sections 8.2,8.3

...but the fiull sum of me
Is sum of nothing;

Shakespeare, The Merchant of Venice, Act 3, Scene 2 (an out of context quote from Portia)

A typical series solution

We need to set up a series to workl with.  This is done with the Sum command, which is similar to the sigma notation from mathematics

y = Sum[c[i] x^i, {i,0,12}] + O[x]^13

OutputForm[y]

The term O[x^13] in both the command and the output is a "big-oh term" of order 13.  This is like the dots in mathematical notation for a series, as we usually write it on the board.  The "big-oh" notation is somewhat more precise, as it has additional information than just the dots.  This terms "holds" and "hides" all of the terms of the series which are of power or order 13 and higher.

The big-oh notation also occurs in other parts of mathematics.   The big-oh notation has a mathematical menaing as indicating a certain type of limit occurs.  Mathematically, it means that all terms remaining converge to zero like a power x^13 or higher.  Although the notation and use are common in higher mathematics we will not use it further here, except within Mathematica.

The use of the big-oh notation in Mathematica has a dual purpose beyond the mathematical notation.  You know that a series is different from a polynomial.   For instance, there are infinitely many terms, even if some of them are "hidden" in the dots notation.  These extra terms are significant when doing a differentiaion, an integration, combining two series, and so forth.  We have seen example sof these operations on series in the notes.

The "big-oh" notation is a way of signalling to Mathematica that the object here is a series.  Naturally, these series contain more terms of higher-order which must be considered in differentiation, integration, multiplication, further expansion and so forth.  

Technically, a series is a special type of  data structure in Mathematica.  This type is differnet from the expressions we usually use, and so must  specifially be indicated to Mathematica.

Now using this series, let's work the following typical problem 
(from Zill, Section 6.2, Number 22).

	(x^2 - 1) y'' + x y'  - y = 0
	
Insert the series into the differential equation.  Label the result of the substituion as de1.

de1 = (x^2 -1) D[y, {x,2}] + x D[y,x] - y == 0

OutputForm[de1]

Noticew that all of the terms have been gathered autiomatically.

Notice also that the order, indicated by the "big-oh" notation has changed.  The order of the substitution into the left side of the equation is now 11.  Effectively, the differentiation twice in the first term, lowered the order of the series by 2.  All other series expressions inserted into the differential equation  had their terms  of order 11 or higher absorbed into the order 11 term.

Now we set all coefficients equal to zero.  Mathematica has a  general purpose commend which will do this automatically, all at once!  See Mathematica, Section 3.6.6 for an explanation of Logical Expand in the contexxt of Series.  ( The command hasother uses as well.)

coeffeqns = LogicalExpand[ de1]

Notice all of the double ampersands &&.  This says the conditions are 
-c[0] - 2c[2] == 0 
AND
-6 c[3] == 0
AND 
so on.  It's just like writing all of the conditions on the board in a row.

Now we want to solve for all the coefficients in terms of the first 2 coefficents c[0] and c[1].  We do this with a variant of the Solve command.  We order Mathematica to solve for coefficeints c[2], c[3], c[4], ...c[12].  This will treat c[0] and c[1] as constants and solve for everything in terms of them.  Rather than exhaustively write out all of c[2], c[3], c[4], c[12] in the Solve command, we use the Table command to automatically generate the table of coefficients we want to solve for.  Neat, isn't it ?!

coeffs = Solve[ coeffeqns, Table[ c[i], {i,2,12}]]

Substitute all of these coefficients back into the series.

y = y /. coeffs

OutputForm[y]

Now we want to pick off the two linearly independent solutions by factoring out the corresponding coefficients.  There are many ways to do this, I choose the following way:

y1 = Coefficient[ y,c[1]]

Here's the second solution,   at least the first 7 terms of the series,, out to order 12.  It looks like a series, but I don't recognize it as being familiar..

y2 = Coefficient[ y, c[0]]

A Further Investigation of the Solution

At this stage,  the problem is done for our purposes in 221.  However, when I first saw this solution, the pattern of the coefficients was so nice, I felt that I ought to recognize the second solution.  Fortunately, I have another way to get the second solution, using the first solution by using the reduction of order technique.

Rather than using the METHOD of reduction of order, as I would if  I were working on the blackboard, I will use the canned formula for the reduction of order.  So I use the formula for the for the second solution from the first found in the green box on page 149 of Nagle & Saff.

y1 = x

p = x/(x^2 -1)

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

That's the result of the order of reduction formula for the second solution, starting from the first fairly simple solution.  Let's test it  to see if it really works in the differential equation.

(x^2 -1) D[y2, {x,2}] + x D[y2,x] - y2

Simplify[%]

Yes it works!!

Now the question remains:  how does this simple solution in closed form relate to the solution in series solution obtained in the previous section?  To answer this, we take the Taylor series expansion ofht eclosed form solution and compare to the series solution.

Use the following command  for automatic Taylor series expansions.  See Mathematica, Section 3.6.

Series[y2, {x,0,12}]

OutputForm[%]

Something's wrong!!  I'm getting lots of imaginary constants I where I didn't expect them!!

If I go back, and examine the solution, I begin to see why.  The solution y2 = Sqrt[-1 + x^2]
is NOT defined for x = 0. ( or more precisely, it is imaginary at x = 0).  With a little thought, I see the solution defined in the neighborhood of x = 0 should be y2 = Sqrt[1 - x^2 ].  Try it in the equation to make sure it works.

y2 = Sqrt[1 - x^2]

(x^2 -1) D[y2, {x,2}] + x D[y2,x] - y2

Simplify[%]

So it really works as a solution!  But does the Taylor series really give the same as we derived for the series solution?

Series[y2, {x,0,12}]

OutputForm[%]

Yes!  It compares the same.
What does this mean?  It probably means that in creating the problems, Zill, worked backward from a known solution, and created a non-constant coefficient differential equation.  Then the solution in terms of series is routine.  

Notice that with Mathematica, we have been able to reverse engineer the solution to discover how the problem was probably created.!

Exercises

Exercise 1

(1) (a) Find approximate series for the two linearly independent solutions to  the equation

	x''(t)  + ( 1- t/100) x(t) = 0
	
by using the power-series method.

(b) Find the solution to the IVP with x(0) =1 and x'(0) = 0.
Don't make this too hard, you should almost be able to do this in your head!

Exercise 2

(2) (a) Find approximate series for the two linearly independent solutions to  the equation

	x''(t)  + ( 1- t/100 + t^2/(2*100^2) x(t) = 0
	
by using the power-series method.

(b) Find the solution to the IVP with x(0) =1 and x'(0) = 0.
Don't make this too hard, you should almost be able to do this in your head!

Particular Case of Legendre's Equation

Legendre's equation is the differential equation
	(1 -x^2) y'' - 2 x y' + n(n+1) y = 0
where n is an integer.  I will examine the case when n = 3.

Legendre's equation comes up in many physical situations involving spherical symmetry.

y = Sum[ c[i] x^i, { i, 0,12}] + O[x]^13

Insert into the differential equation

de = (1 -x^2) D[y, {x,2}] - 2 x D[y,x] + 12 y == 0

Get the coefficient equations

coeffeqns = LogicalExpand[ de]

solvedcoeffs = Solve[ coeffeqns, Table[ c[i], {i,2,12}]]

Put these back into the series expansion.

y = y /. solvedcoeffs

Extract out the two linearly independent solutions

Coefficient[ y, c[0]]

Coefficient[ y, c[1]]

This solution is actiually a polynomial, in this case of degree 3.  By solving Legendre's equation equation for any integer n, we obtain the full set of corresponding Legendre polynomials of all degrees.  The Legendre polynomials are anoother important set of mathematical functions with amazing and wonderful properties.  In fact, they are so useful, they are built into Mathematica automatically.

LegendreP[3,x]

This is different from the polynomial obtained earlier, but only by a constant factor.  The constant factor is used to set the normailization for the Legendre polynomials.