unit31.html

unit31.mws

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 31 -- Power Series Solutions

Prof. Douglas B. Meade

Industrial Mathematics Institute

Department of Mathematics

University of South Carolina

Columbia, SC 29208

URL: http://www.math.sc.edu/~meade/

E-mail: meade@math.sc.edu

-------------------------------------------------------------------

>

Outline of Unit 31

31.A Taylor Series Methods for Solving IVPs
31.B Power Series Solutions
31.C Ordinary and Singular Points
31.C-1 Example 1: Polynomial Coefficients
31.C-2 Example 2: Rational Coefficients
31.C-3 Example 3: Legendre's Equation
31.D Method of Frobenius
31.D-1 Example 1: Two Real Roots with a Noninteger Difference
31.E Legendre's Equation
31.F Bessel's Equation

>

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

> with( orthopoly ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

31.A Taylor Series Methods for Solving IVPs

Consider the initial value problem

> ode1 := (x+y(x)+1) * diff( y(x), x ) = 1;

> ic1 := y(0) = 1;

>

If the solution to this IVP has a Taylor series expansion at the initial point, , then approximations to the solution can be obtained from the corresponding Taylor polynomials. For example, the quintic Taylor polynomial is

> N := 5:

> Order := N+1:

> form1 := y(x) = convert( series( y(x), x=0 ), polynom );

>

To completely determine this polynomial, it is necessary to compute the value of , , ..., . The first of these values is simply the initial condition.

> Dy0[0] := ic1;

>

The second can be obtained by substituting the initial condition into the differential equation:

> Dode1[0] := simplify( convert( ode1, D ) );

> Dy0[1] := isolate( subs( x=0, ic1, Dode1[0] ), D(y)(0) );

>

To compute the value of , differentiate the ODE with respect to and substitute and the lower-order derivatives to obtain

> Dode1[1] := simplify( convert( diff( ode1, x ), D ) ):

> Dy0[2] := isolate( subs( x=0, Dy0[i] $ i=0..1, Dode1[1] ), (D@@2)(y)(0) );

>

This process can be continued to determine the value of the remaining derivatives of the solution at the initial point.

> for n from 2 to N-1 do

> Dode1[n] := simplify( convert( diff( ode1, x$n ), D ) );

> Dy0[n+1] := isolate( subs( x=0, Dy0[i] $ i=0..n, Dode1[n] ), (D@@(n+1))(y)(0) );

> end do:

> unassign( 'n' );

> print( Dy0[n] $ n=0..N );

The approximate solution obtained by this method is

> sol1 := eval( form1, [ Dy0[i] $ i=0..N ] );

>

The same solution can be obtained directly from Maple when type=series is included in the argument list for dsolve :

> infolevel[dsolve] := 3:

> sol1a := dsolve( { ode1, ic1 }, y(x), type=series );

> infolevel[dsolve] := 0:

dsolve/series/single:   initial conditions    [[], {1}]

unknown:   non-singular case

dsolve/series/sysol:   Newton iteration...

dsolve/series/sysol:   latest approx is   [1/2*X+1]

dsolve/series/sysol:   latest approx is   [1/2*X+1-3/16*X^2+7/64*X^3]

dsolve/series/sysol:   latest approx is   [1/2*X+1-3/16*X^2+7/64*X^3-79/1024*X^4+1237/20480*X^5]

>

Note the presence of the "order" term in solution. This indicates that the result is an approximation. To plot this solution it is necessary to eliminate the order term; the simplest way to do this is to convert the expression to a polynomial:

> sol1b := convert( sol1a, polynom );

>

The degree of Maple's series solution is determined by the value of the global variable Order . Thus, the ninth-degree Taylor polynomial approximate solution to this IVP can be obtained with

> Order := 10:

> convert( dsolve( { ode1, ic1 }, y(x), type=series ), polynom );

>

Note that this method extends to IVPs for higher-order differential equations provided it is possible to express the highest order derivative as a differentiable function of the lower order derivatives of the unknown function (and the independent variable). For example,

> ode2 := ( x^2 + 1 ) * diff( y(x), x$2 ) + x * diff( y(x), x ) + y(x) = 0;

> Order := 12:

> ic2 := y(0)=1, D(y)(0)=-1;

> sol2 := convert( dsolve( { ode2, ic2 }, y(x), type=series ), polynom );

> Order := 6:

>

The computation of Taylor series expansions is typically limited by the ability to compute and manipulate the derivatives of the ODE. A tool such as Maple is very helpful in this regard, but other methods can be more effective in many situations.

>

31.B Power Series Solutions

Another problem with the Taylor polynomials is that they are approximate solutions. The exact solution can be obtained only when a general pattern for the coefficients in the expansion can be determined. This general pattern to the coefficients can often be found directly by substituting a general power series into the ODE.

Consider the problem of finding the general solution to the first-order (linear) ODE

> ode3 := diff( y(x), x ) - 2 * y(x) = 0;

>

The general form of a power series solution, based at the origin, is

> form3 := y(x) = Sum( a(k)*x^k, k=0..infinity );

>

When this solution is substituted into the ODE and simplified, the result is :

> q1 := eval( ode3, form3 );

> q2 := combine( q1 );

>

While Maple is not quite up to manipulating infinite power series to find the general relation between the coefficients, a little insight and skill can be used to obtain this information. Because this equation is first order and, as seen above, only consecutive pairs of coefficients are related, it suffices to insert the general three-term sum into the ODE and manipulate the expression to obtain

> tform3 := y(x) = Sum( a(n)*x^n, n=k-1..k+1 );

> q1 := eval( ode3, tform3 );

> q2 := simplify( combine( value( q1 ), power ) );

>

The recurrence equation between the coefficients can now be found to be

> recur_eqn3 := map( coeff, q2, x^k );

>

The (general) solution to this recurrence equation can be found with Maple's rsolve command

> recur_sol3 := rsolve( recur_eqn3, {a} );

$recur_sol3 := {a(k) = 2^k*a(0)/GAMMA(k+1)}$

>

The general power series solution to the ODE is

> sol3 := subs( recur_sol3, form3 );

>

To express the function, GAMMA , in terms of factorials ,

> convert( sol3, factorial );

>

Observe that attempts to verify that these functions are solutions to the ODE via odetest requires a little additional work :

> q3 := odetest( sol3, ode3 );

> value( q3 );

>

Because of the simplicity of this example, it is possible to identify the function with this power series:

> value( sol3 );

>

The constant takes the place of the integration constant in the other solution methods. Note that . With this insight, this result is seen to be equivalent to the explicit solution found using dsolve :

> dsolve( ode3, y(x) );

>

There are many examples for which rsolve will not be able to solve the recurrence relation. In that case, it will be necessary to construct a finite set of explicit equations and solve for the first few coefficients :

> recur_sol3t := solve( { recur_eqn3 $ k=0..5 }, { a(k) $ k=1..6 } );

$recur_sol3t := {a(6) = 4/45*a(0), a(2) = 2*a(0), a(...$

> collect( eval( value(subs(infinity=6,form3)), recur_sol3t ), a(0) );

>

31.C Ordinary and Singular Points

Consider the homogeneous second-order linear ODE in normalized form

> ode4 := diff( y(x), x$2 ) + PP(x) * diff( y(x), x ) + Q(x) * y(x) = 0;

>

(Note that the orthopoly package uses the name P for the Legendre polynomials; hence the use of the name PP for the coefficient in the ODE.)

If initial conditions are given at :

> ic4 := y(x0) = y0, D(y)(x0) = y1;

>

then a power series solution should be based at the point :

> form4 := Sum( a(k) * ( x - x0 )^k, k=0..infinity );

>

The point is called an ordinary point of the ODE if the functions and can be expressed as power series based at , i.e. , and are real analytic at . If at least one of and is not real analytic at , then is called a singular point of the ODE. A singular point is called a regular singular point if and are real analytic at , otherwise is an irregular singular point.

>

31.C-1 Example 1: Polynomial Coefficients

When the coefficient functions are polynomials, e.g. ,

> coeff1 := { PP(x) = 1 - x^3, Q(x) = x^2 };

$coeff1 := {PP(x) = 1-x^3, Q(x) = x^2}$

> eval( ode4, coeff1 );

>

all points are ordinary points.

>

31.C-2 Example 2: Rational Coefficients

When the coefficient functions are rational functions, e.g. ,

> coeff2 := { PP(x) = 2/(1 - x)^2, Q(x) = x^(-2)/(x-2) };

$coeff2 := {Q(x) = 1/(x^2*(x-2)), PP(x) = 2*1/((-x+1...$

> eval( ode4, coeff2 );

>

all discontinuities of or of are singular points:

> singpts := `union`(op( map( c -> discont( rhs( c ), x ), coeff2 ) ));

>

To decide if the singular points are regular or irregular, the definition can be applied or Maple's regularsp command (from the DEtools package) can be used:

> regularsp( eval( ode4, coeff2 ), x, y(x) );

>

and so must be an irregular singular point.

>

31.C-3 Example 3: Legendre's Equation

Legendre's equation of order is

> ode5 := ( 1 - x^2 ) * diff( y(x), x$2 ) - 2*x * diff( y(x), x ) + p*(p+1) * y(x) = 0;

>

where is a nonnegative constant. Written in normalized form, the coefficients and are

> q4 := eval( lhs(ode5), [ diff( y(x), x$2 ) = d2y, diff( y(x), x ) = dy, y(x)=y ] );

> coeff3 := { PP(x) = coeff( q4, dy )/coeff( q4, d2y ), Q(x) = coeff( q4, y )/coeff( q4, d2y ) };

$coeff3 := {PP(x) = -2*x/(1-x^2), Q(x) = p*(p+1)/(1-...$

>

The singular points are

> singpts := `union`(op( map( c -> discont( rhs( c ), x ), coeff3 ) ));

>

and both of these points are regular singular points

> regularsp( ode5, x, y(x) );

>

31.D Method of Frobenius

The method for finding a power series solution presented in Unit 31, Section B can be applied at any ordinary point of an ODE. The Method of Frobenius is guaranteed to find at least one nontrivial solution in a neighborhood of a regular singular point. When the Method of Frobenius does not produce a second solution, reduction of order can be used (see Unit 19). For an irregular singular point, the Method of Frobenius can be tried but may not produce a solution in all cases.

The Method of Frobenius produces solutions of the form

> form5 := y(x) = Sum( a(k) * (x-x0)^(k+r), k=0..infinity );

form5 := y(x) = Sum(a(k)*(x-x0)^(k+r),k = 0 .. infi...

>

When is a nonnegative integer, the Frobenius solution is a power series, but does not have to be an integer -- it can be a negative or complex number. When a Frobenius solution is inserted into the ODE, the condition for the vanishing of the lowest order term is a quadratic polynomial in . The general form of the indicial equation is

r^2+(p[0]-1)*r+q[0] = 0

where and are the constant terms in the power series expansions of and q(x) = x^2*Q(x) , respectively.

The full Method of Frobenius is long and exceedingly complicated. In the remainder of this worksheet the basic ideas will be explored in a number of examples.

>

31.D-1 Example 1: Two Real Roots with a Noninteger Difference

Consider the ODE

> ode6 := 2*x*(1+x) * diff( y(x), x$2 ) + (3+x) * diff( y(x), x ) - x * y(x) = 0;

>

The regular singular points for this problem are

> regularsp( ode6, x, y(x) );

>

The indicial equation at the regular singular point is

> ie1 := indicialeq( ode6, x, 0, y(x) );

>

The indicial roots are

> ir1 := solve( ie1, x );

>

Because these roots do not differ by an integer, each indicial root should lead to a nontrivial Frobenius solution

> formF := y(x) = Sum( c(k) * x^(k+r), k=0..infinity ):

> formF1 := eval( formF, [ c=a, r=max(ir1) ] );

> formF2 := eval( formF, [ c=b, r=min(ir1) ] );

>

To obtain the recurrence relations for the coefficients in the two solutions, work with truncated versions of the solution

> tformF1 := eval( y(x) = Sum( a(m) * x^(m+r), m=k-2..k+2 ), r=max(ir1) );

> tformF2 := eval( y(x) = Sum( b(m) * x^(m+r), m=k-2..k+2 ), r=min(ir1) );

> q51 := simplify( combine( value( eval( ode6, tformF1 ) ), power ) );

> q52 := simplify( combine( value( eval( ode6, tformF2 ) ), power ) );

> reqn1 := collect( map( coeff, q51, eval( x^(k+r), r=max(ir1) ) ),

> [ a(k+i) $ i=-3..3 ] );

> reqn2 := collect( map( coeff, q52, eval( x^(k+r), r=min(ir1) ) ),

> [ b(k+i) $ i=-3..3 ] );

>

While these recurrence relations cannot be solved explicitly by Maple

> rsolve( reqn1, {a} );

> rsolve( reqn2, {b} );

>

The first few terms are easily computed

> q61 := isolate( reqn1, a(k+1) );

> q62 := isolate( reqn2, b(k+1) );

> unassign( 'a', 'b' );

> for k from 1 to 8 do

> assign( { q61, q62} );

> print( 'a'(k+1) = a(k+1), 'b'(k+1) = b(k+1) );

> end do;

> unassign( 'k' );

>

and so the (approximate) solutions to the ODE are

> q71 := value( subs( infinity=9, formF1 ) );

> q72 := value( subs( infinity=9, formF2 ) );

>

> sol61 := sort( collect( q71, [ a(0), a(1) ] ) );

> sol62 := sort( collect( q72, [ b(0), b(1) ] ) );

>

To conclude, compare this solution to the series solution obtained by dsolve with type=series :

> Order := 10:

> sol6 := dsolve( ode6, y(x), type=series );

> Order := 6:

>

To see the equivalence of these solutions, observe that the sum of the two linearly independent solutions obtained from the Method of Frobenius agrees with the dsolve solution for a specific choice of , , , and :

> q8 := rhs(sol61) + rhs(sol62) = convert( rhs(sol6), polynom );

> solve( identity(q8,x), {a(0),a(1),b(0),b(1)} );

${a(1) = 0, a(0) = _C2, b(0) = _C1, b(1) = -_C1}$

>

31.E Legendre's Equation

Recall, from Section 31.C-3, that Legendre's equation of order (a nonnegative constant),

> ode5;

(1-x^2)*diff(y(x),`$`(x,2))-2*diff(y(x),x)*x+p*(p+1...

>

has two regular singular points

> regularsp( ode5, x, y(x) );

>

Because , is an ordinary point, there should be no trouble simply finding the general solution in the form of a power series :

> form5 := y(x) = Sum( a(k) * x^k, k=0..infinity );

form5 := y(x) = Sum(a(k)*x^k,k = 0 .. infinity)

>

To obtain the recurrence relation for the coefficients in this solution, work with a truncated version of the solution

> tform5 := y(x) = Sum( a(m) * x^m, m=k-2..k+2 );

tform5 := y(x) = Sum(a(m)*x^m,m = k-2 .. k+2)

> q9 := simplify( combine( value( eval( ode5, tform5 ) ), power ) );

> reqn3 := collect( map( coeff, q9, x^k ), [ a(k+i) $ i=-3..3 ] );

>

While this recurrence relation cannot be solved explicitly by Maple

> rsolve( reqn3, {a} );

$rsolve((p-k^2-k+p^2)*a(k)+(2+3*k+k^2)*a(k+2) = 0,{a...$

>

The first few terms are easily computed

> q10 := factor( isolate( reqn3, a(k+2) ) );

>

> unassign( 'a' );

> for k from 0 to 6 do

> assign( { q10 } );

> print( 'a'(k+2) = a(k+2) );

> end do;

> unassign( 'k' );

>

and so the (approximate) solutions to the ODE are

> q11 := sort( collect( value( subs( infinity=8, form5 ) ), [ a(0), a(1) ] ) );

>

To confirm this computation, compare this solution to the series solution obtained by dsolve with type=series :

> Order := 9:

> sol5 := collect( convert( dsolve( ode5, y(x), type=series ), polynom ), [ y(0), D(y)(0) ] );

> Order := 6:

>

> q12 := rhs(q11) = rhs(sol5):

> solve( identity(q12,x), {a(0),a(1)} );

>

While we have not proven it, the series associated with and the series associated with , i.e. , the two functions

> sol51 := y[1](x) = eval( rhs(q11), [ a(0)=1, a(1)=0 ] );

> sol52 := y[2](x) = eval( rhs(q11), [ a(0)=0, a(1)=1 ] );

>

form a fundamental set of solutions to Legendre's equation of order .

>

Additional inspection of these solutions (or the recurrence relations) when is an integer reveals that either the series of even terms ( ) or the series of odd terms ( ) terminates. That is, one fundamental solution is a polynomial. This polynomial, normalized to have value 1 at , is the Legendre polynomial of order .

> for k from 0 to 5 do

> if type(k,even) then

> q := eval( rhs(sol51), p=k );

> else

> q := eval( rhs(sol52), p=k );

> end if;

> P[k] = q / eval( q, x=1 );

> end do;

>

The Legendre polynomials form an orthogonal family of functions in the sense that

Int(P[m](x)*P[n](x),x = -1 .. 1) = 0 for all nonnegative integers .

Maple provides the Legendre polynomials in the orthopoly package :

> for k from 0 to 5 do

> P( k, x );

> end do;

> unassign( 'k' ):

>

> plot( [ P(k,x) $ k=0..6 ], x=-1..1, legend=[ seq( "P"||k, k=0..6 ) ] );

[Maple Plot]

>

The orthogonality of the (first few) Legendre polynomials is confirmed by the diagonal structure in the matrix of inner products of pairs of Legendre polynomials

> matrix( 9, 9, (m,n) -> int( P(m,x)*P(n,x), x=-1..1 ) );

matrix([[2/3, 0, 0, 0, 0, 0, 0, 0, 0], [0, 2/5, 0, ...

>

The LegendreP and LegendreQ commands provide high level access to the Legendre functions of the first and second kind, respectively. The series command can be used to obtain a series expansion for the appropriate Legendre function.

> series( LegendreP(1/2,x), x=0 );

series(2*GAMMA(3/4)^2/(Pi^(3/2))+1/2*sqrt(Pi)/(GAMM...

>

31.F Bessel's Equation

Bessel's equation of order is

> ode7 := x^2 * diff( y(x), x$2 ) + x * diff( y(x), x ) + ( x^2 - nu^2 ) * y(x) = 0;

>

where is a nonnegative constant. All points are ordinary points except for the regular singular point at the origin.

> regularsp( ode7, x, y(x) );

>

Thus, a series solution at the origin can be found with the Method of Frobenius. The indicial equation at the origin is

> ie2 := indicialeq( ode7, x, 0, y(x) );

>

with indicial roots

> ir2 := solve( ie2, x );

>

Because is a regular singular point for Bessel's equation, solutions are sought via the Method of Frobenius

> formF := y(x) = Sum( c(k) * x^(k+r), k=0..infinity );

> formF3 := eval( formF, [ c=a, r=nu ] );

> formF4 := eval( formF, [ c=b, r=-nu ] );

>

To obtain the recurrence relations for the coefficients in the two solutions, work with truncated versions of the solution

> tformF3 := eval( y(x) = Sum( a(m) * x^(m+r), m=k-2..k+2 ), r=nu );

> tformF4 := eval( y(x) = Sum( b(m) * x^(m+r), m=k-2..k+2 ), r=-nu );

> q131 := simplify( combine( value( eval( ode7, tformF3 ) ), power ) );

> q132 := simplify( combine( value( eval( ode7, tformF4 ) ), power ) );

> reqn3 := collect( map( coeff, q131, eval( x^(k+r), r=nu ) ), [ a(k+i) $ i=-3..3 ] );

> reqn4 := collect( map( coeff, q132, eval( x^(k+r), r=-nu ) ), [ b(k+i) $ i=-3..3 ] );

>

While these recurrence relations cannot be solved explicitly by Maple

> rsolve( reqn3, {a} );

> rsolve( reqn4, {b} );

$rsolve(a(k-2)+(2*k*nu+k^2)*a(k) = 0,{a})$

$rsolve(b(k-2)+(-2*k*nu+k^2)*b(k) = 0,{b})$

>

The first few terms are easily computed

> q141 := factor( isolate( reqn3, a(k) ) );

> q142 := factor( isolate( reqn4, b(k) ) );

>

> unassign( 'a', 'b' );

> for k from 2 to 9 do

> assign( { q141, q142} );

> print( 'a'(k) = a(k), 'b'(k) = b(k) );

> end do;

> unassign( 'k' );

>

and so the (approximate) solutions to the ODE are

> q151 := value( subs( infinity=9, formF3 ) );

> q152 := value( subs( infinity=9, formF4 ) );

>

> sol71 := sort( collect( q151, [ a(0), a(1) ] ) );

> sol72 := sort( collect( q152, [ b(0), b(1) ] ) );

>

To confirm these computations, compare these solution to the series solution obtained by dsolve with type=series :

> Order := 10:

> sol7 := dsolve( ode7, y(x), type=series );

> Order := 6:

>

> q16 := rhs(sol71) + rhs(sol72) = convert( rhs(sol7), polynom ):

> solve( identity(q16,x), {a(0),a(1),b(0),b(1)} );

${b(0) = _C2, b(1) = 0, a(1) = 0, a(0) = _C1}$

>

Note that these series solutions are not well defined when is an integer. When is not an integer, a fundamental set of solutions to Bessel's equation of order nu is { , }

Define the Bessel's function (of the first kind) of order and ( nonnegative) to be

> J[nu](x) = eval( rhs(sol71), [ a(0)=1, a(1)=0 ] );

> J[-nu](x) = eval( rhs(sol72), [ b(0)=1, b(1)=0 ] );

>

While it is tempting to state that is a fundamental set of solutions to Bessel's equation of order , note that and that is not well-defined when is an even positive integer. If is not a positive integer, then is a fundamental set of solutions to Bessel's equation of order . If is a positive integer, a fundamental set of solutions is where , the Bessel function of the second kind of order , is defined for noninteger orders to be a linear combination of the linearly independent functions and :

and, for integer orders ( ) :

Thus, the general solution to Bessel's equation of (nonnegative) order can be written as

Direct access to the Bessel functions of the first and second kind is provided with Maple's BesselJ and BesselY commands.

Two important properties of Bessel functions are that and . This property is essentially a result of the fact that is obtained from the indicial root > 0 and from the indicial root < 0. These properties can be observed in the following plots :