Maple_A1.mw

Appendix A1. Power Series and Special Functions

Series Solutions

Maple's dsolve procedure can be used to obtain series solutions to differential equations. Add the optional equation

type=series

The unknown function must also be inserted directly after the IVP.

Example  Obtain a series solution to the following differential equation (Ledder, A4, Exercise 2)

y'' + 2x y' - y = 0

Then obtain the first 6 non-zero terms of the series solution satisfying y(0) = 1, y'(0) = 3. Plot the approximate and the exact solutions. (Think of a mass spring system with a "repelling" spring moving in a fluid that is thickening with time.)

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

DE := (diff(y(x), `$`(x, 2)))+2*x*(diff(y(x), x))-y(x) = 0

> gen_soln := dsolve( DE, y(x), type=series);

gen_soln := y(x) = (series(y(0)+D(y)(0)*x+1/2*y(0)*x^2-1/6*D(y)(0)*x^3-1/8*y(0)*x^4+1/24*D(y)(0)*x^5+O(x^6),x,6))

Lacking initial conditions, Maple finds the series solution near x = 0. This is why the general solution formula is given in terms of the initial values y(0) and D(y)(0). The "order 6" solution is obtained by default. The term

O(x^6)

represents an error term that, as x approaches 0, also goes to 0 at least as fast as x^6 does. To obtain approximate solution values (or curves) convert the output to a polynomial.

> convert(gen_soln,polynom);

y(x) = y(0)+x*D(y)(0)+1/2*y(0)*x^2-1/6*D(y)(0)*x^3-1/8*y(0)*x^4+1/24*D(y)(0)*x^5

If initial values are given, they are substituted into the series. To obtain the first 6 non-zero terms ask for the order 7 solution by inserting the equation

order=7

> soln := dsolve( {DE,y(0)=1,D(y)(0)=3}, y(x), type=series, order=7);

soln := y(x) = (series(1+3*x+1/2*x^2-1/2*x^3-1/8*x^4+1/8*x^5+7/240*x^6+O(x^7),x,7))

> y_approx := convert(rhs(soln),polynom);

y_approx := 1+3*x+1/2*x^2-1/2*x^3-1/8*x^4+1/8*x^5+7/240*x^6

The next entry generates an exact solution formula. It is rather complicated, given in terms of Bessel functions, exponential functions, and the Gamma function.

> exact_soln := dsolve( {DE,y(0)=1,D(y)(0)=3} );

exact_soln := y(x) = 1/4*2^(1/2)*(Pi*2^(1/2)+6*GAMMA(3/4)^2)*exp(-1/2*x^2)*x^(3/2)*(BesselI((-1)/4, 1/2*x^2)+BesselI(3/4, 1/2*x^2))/GAMMA(3/4)-1/2*2^(1/2)*exp(-1/2*x^2)*x^(3/2)*(BesselK(1/4, 1/2*x^2)-...exact_soln := y(x) = 1/4*2^(1/2)*(Pi*2^(1/2)+6*GAMMA(3/4)^2)*exp(-1/2*x^2)*x^(3/2)*(BesselI((-1)/4, 1/2*x^2)+BesselI(3/4, 1/2*x^2))/GAMMA(3/4)-1/2*2^(1/2)*exp(-1/2*x^2)*x^(3/2)*(BesselK(1/4, 1/2*x^2)-...

Compare the exact and approximate solution curves. The approximation is the blue dashed curve (linestyle=3).

> plot( [rhs(exact_soln), y_approx], x=0..3, 0..8, color=[red,blue], linestyle=[1,3]);

[Plot]

Want a recursion relation: Use powsolve

When the differential equation has polynomial coefficients and x = 0 is not a singular point, the recursion relation can be obtained by applying the powsolve procedure in the powseries package.

> with(powseries);

[compose, evalpow, inverse, multconst, multiply, negative, powadd, powcos, powcreate, powdiff, powexp, powint, powlog, powpoly, powsin, powsolve, powsqrt, quotient, reversion, subtract, template, tpsf...[compose, evalpow, inverse, multconst, multiply, negative, powadd, powcos, powcreate, powdiff, powexp, powint, powlog, powpoly, powsin, powsolve, powsqrt, quotient, reversion, subtract, template, tpsf...

> pow_soln := powsolve( {DE, y(0)=1, D(y)(0)=3} );

pow_soln := proc (powparm) local nn, t1; option `Copyright (c) 1990 by the University of Waterloo. All rights reserved.`; table( [( 0 ) = 1, ( 1 ) = 3, ( _k ) = -(-5+2*_k)*a(_k-2)/(_k*(_k-1)) ] ) if t...

The output is a procedure that can be used to generate power series solutions.

Apply the tpsform procedure ("to power series form") to pow_soln as follows to generate the 10th order solution formula

tpsform(pow_soln,x,10)

> tpsform(pow_soln,x,10);

series(1+3*x+1/2*x^2-1/2*x^3-1/8*x^4+1/8*x^5+7/240*x^6-3/112*x^7-11/1920*x^8+13/2688*x^9+O(x^10),x,10)

Apply pow_soln directly to the expression _k as shown below to obtain the recursion relation.

> pow_soln(_k);

-(-5+2*_k)*a(_k-2)/(_k*(_k-1))

Interpret this output as

a[k] = -(-5+2*k)*a[k-2]/(k*(k-1))

For numbers and graphs, you may be better off with dsolve/numeric

The dsolve/numeric procedure and odeplot will also produce a good picture of the exact solution curve.

> numeric_soln := dsolve( {DE, y(0)=1, D(y)(0)=3}, y(x), numeric);

numeric_soln := proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits :...

> with(plots):
display( odeplot( numeric_soln, [x,y(x)], x=-3..3, color=red),

        plot( y_approx, x=-3..3, -8..8, color=blue, linestyle=3) );

[Plot]

Here is now the value of the numeric solution compares to the exact value, say at x = 3.

> Numeric,numeric_soln(3);
Exact,eval(exact_soln,x=3): evalf(%);

Numeric, [x = 3., y(x) = 7.52974309807350828, diff(y(x), x) = 1.29511634649786854]

Exact, y(3) = 7.529743410

Special functions

Maple has, as built in procedures, all of the special functions that are found in most elementary differential equations textbooks. For example, the Bessel functions of the first and second kind of order n are denoted BesselJ(n,x) and BesselY(n,x) respectively. Compare the plots below to the the ones appearing in Chapter 3, Figures 3.7.3 and 3.7.4 in Ledder.

> plot( [BesselJ(0,x),BesselJ(1,x)], x=0..20, color=red, linestyle=[1,3]);

[Plot]

> plot( [BesselY(0,x),BesselY(1,x)], x=0..20, -1..0.6, color=red, linestyle=[1,3]);

[Plot]

To see a list of all the functions that are initially known to Maple (without loading any special packages) go here.

> ?inifcn

If you would like to get more information about the functions that are "built in" to Maple go see the FunctionAdvisor.

> ?FunctionAdvisor