Finding Power Series Solutions Differential Equations

Lesson31.mw

• Outline of Lesson 31
• Initialization
• 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 4: Two Real Roots with a Noninteger Difference
• 31.E Legendre's Equation
• 31.F Bessel's Equation

---

ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Lesson 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

Copyright  2001  by Douglas B. Meade

All rights reserved

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

Outline of Lesson 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 4: Two Real Roots with a Noninteger Difference

31.E Legendre's Equation

31.F Bessel's Equation

Initialization

Warning, the name changecoords has been redefined

31.A Taylor Series Methods for Solving IVPs

Consider the initial value problem

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

If the solution of this IVP has a Taylor series expansion at the initial point, , then the Taylor polynomials generated from this solution are approximations of the solution. For example, the formal quintic Taylor polynomial for the function is

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

To completely determine this polynomial, it is necessary to compute the numbers , , , .

If is to be the solution of the IVP, then the first number to compute is simply the initial condition,

The second number can be obtained by substituting the initial condition into the differential equation, and solving for , which is thereby found to be

>	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 values already computed for 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.  The result is

>	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) );

>	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 .  The result is

>	infolevel[dsolve] := 3:

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

>	infolevel[dsolve] := 0:

`DEtools/convertsys:`, `converted to first-order system Y'(x) = f(x,Y(x)\
)`, `   namely (with Y' represented by YP)`

`DEtools/convertsys:`, `correspondence between Y[i] names and original f\
unctions:`

`dsolve/series/ordinary:`, `converted to first-order system Y'(x) = f(x,\
Y(x))`, `   namely (with Y' represented by YP)`

`dsolve/series/ordinary:`, `correspondence between Y[i] names and origin\
al functions:`

`dsolve/series/ordinary:`, `vector Y of initial conditions at x0 =`, 0, \
[1]
`dsolve/series/ordinary:`, `trying Newton iteration`
`dsolve/series/sysol:`, `latest approx is`, table( [( 1 ) = 1+1/2*x ] )
`dsolve/series/sysol:`, `latest approx is`, table( [( 1 ) = 1+1/2*x-3/16\
*x^2+7/64*x^3 ] )
`dsolve/series/sysol:`, `latest approx is`, table( [( 1 ) = 1+1/2*x-3/16\
*x^2+7/64*x^3-79/1024*x^4+1237/20480*x^5-24817/491520*x^6 ] )

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 the 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

>	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, an 11th-degree polynomial approximation for the solution of the second-order IVP

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

can be obtained with the Maple commands

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

The computation of Taylor series expansions is typically limited by the complexity and tedium of computing and manipulating 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 (represented by the complete infinite power series), 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.  This is particularly true of linear ODEs.

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 );

or even

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 differential equation is first order and, as seen above, only consecutive pairs of coefficients are related, it suffices to insert the general three-term sum

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

into the ODE and manipulate the expression to obtain

>	q1 := eval( ode3, tform3 );

or even

>	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, and it is

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

The general power series solution to the ODE is

>	sol3 := subs( recur_sol3, form3 );

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

>	convert( sol3, factorial );

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

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 .  This solution is

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. For example, the first few equations generated by the recurrence relation in this example are

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

Their solution is found to generate the polynomial

>	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, namely,

>	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 the initial conditions given at are

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

then a power series solution based at the point would be of the form

>	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 (RSP) if and are real analytic at , otherwise is an irregular singular point (ISP).

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 };

so the ODE becomes

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) };

then the ODE becomes

and all discontinuities of or of are singular points.  Maple finds these singular points to be

>	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.  This latter approach results in

>	regularsp( eval( ode4, coeff2 ), 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 ) };

The singular points are

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

and both of these points are regular singular points, as deduced by

>	regularsp( ode5, y(x) );

31.D Method of Frobenius

The method for finding a power series solution presented in Lesson 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 Lesson 27 ). 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 );

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 real number or even a 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

where and are the constant terms in the power series expansions of and , 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 4: 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, 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.  The forms of these two solutions would be

>	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 solutions, namely,

>	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) );

Substitution into the ODE leads to the equations

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

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

from which are extracted the recurrence relations

>	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, as verified by the attempts

the first few terms are easily computed from

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

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

and are found to be

>	assign( { q61, q62} );

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

and so the (approximate) solutions to the ODE are

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

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

or even

>	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 , namely,

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

To see the equivalence of these solutions,  observe that the sum of the two linearly independent solutions obtained from the Method of Frobenius, namely,

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

agrees with the dsolve solution for a specific choice of , , , and , as seen via

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

31.E Legendre's Equation

Recall, from Section 31.C-3 , that

Legendre's equation of order (a nonnegative constant), has two regular singular points in the finite plane, as shown by

>	regularsp( ode5, y(x) );

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

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

constructed at .

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

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

Substitution into the ODE yields

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

so the recurrence relation is

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

Maple can now solve this recurrence relation, finding

>	ak := rsolve( reqn3, {a} );

and so two independent (approximate) solutions to the ODE are

>	sol1 := collect(simplify(convert(add(ak*x^k,k=[0,2,4,6,8]),factorial)),x,factor); sol2 := collect(simplify(convert(add(ak*x^k,k=[1,3,5,7,9]),factorial)),x,factor);

To confirm these computations, obtain by dsolve with type=series , the series solution that satisfies the initial conditions and .  This solution is

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

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 up to the order of the terms computed.

>	simplify(sol1+sol2-rhs(sol5));

While not proven here, the series associated with and the series associated with , i.e. , the two functions

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. Thus, one fundamental solution is a polynomial. This polynomial, normalized to have value 1 at , is the Legendre polynomial of order .

The first few of these polynomials are

>	q := eval( sol1, p=k );

>	q := eval( sol2, p=k );

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

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

for all nonnegative integers .

Maple provides the Legendre polynomials in the orthopoly package, using the name for .  Maple's list of Legendre polynomials begins with

A graph of these first Legendre polynomials is seen in Figure 31.1.

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

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, computed efficiently with

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

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.  For example, the series for begins

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

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. This is verified with

>	regularsp( ode7, 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.  Thus, there should be two formal power series solutions of the form

>	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 the following 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 );

Substitution into Bessel's equation leads to

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

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

from which the recurrence relations

>	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 ] );

or even

>	factor( isolate( reqn3, a(k) ) ); factor( isolate( reqn4, b(k) ) );

are obtained.

Once again, Maple can now solve the recurrence equations, giving the solutions

>	ak := rsolve( reqn3, {a} );

>	bk := rsolve( reqn4, {b} );

The two Frobenius solutions are therefore

>	sol1 := x^nu   *collect(simplify(convert(add(ak*x^k,k=0..9),factorial)),x,factor); sol2 := x^(-nu)*collect(simplify(convert(add(bk*x^k,k=0..9),factorial)),x,factor);

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

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

Noting that only even powers of appear within the parentheses, modify the Frobenius solutions to read

>	sol11 := eval(sol1,a(1)=0); sol22 := eval(sol2,b(1)=0);

then compare the expressions that are multiplied by and .  If the constants are factored out, or equivalently set to 1, the differences in these expressions are found to be

>	simplify(eval(rhs(sol7),{_C2=0,_C1=1})-eval(sol11,a(0)=1)); simplify(eval(rhs(sol7),{_C2=1,_C1=0})-eval(sol22,b(0)=1));

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(sol11, a(0)=1); J[-nu](x) = eval(sol22, b(0)=1);

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 the following 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

These limits essentially result from being obtained from the indicial root > 0 and from the indicial root < 0. These properties can be observed in Figure 31.2

>	plot( [ BesselJ(1,x), BesselY(1,x) ], x=0..10, y=-2..1,

>	legend=[ "J[1]", "Y[1]" ], title="Figure 31.2" );

and Figure 31.3.

>	plot( [ BesselJ(7/3,x), BesselY(7/3,x) ], x=0..10, y=-2..1,

>	legend=[ "J[7/3]", "Y[7/3]" ], title="Figure 31.3" );

As a final comment, the Bessel functions of half-integer orders can be expressed in unexpectedly simple forms. For example it turns out that

>	for k from -5 to 5 by 2 do

>	J[k/2](x) = BesselJ( k/2, x );

[Back to ODE Powertool Table of Contents ]

Finding Power Series Solutions Differential Equations

Source: https://www.maplesoft.com/applications/view.aspx?SID=4719&view=html&L=E

Share :