# 1.23 Lesson 23

```# Lesson 23.  Power Series
# ------------------------
> restart:
# A power series about the point x=c is a series of the form
> Sum(a[n]*(x-c)^n, n=0..infinity);

infinity
-----
\                  n
)     a[n] (x - c)
/
-----
n = 0
# A power series has a radius of convergence R (possibly 0 or
# infinity), such that the series converges absolutely when
# |x-c|<R and diverges when |x-c|>R.
# The Taylor series of f(x) about the point x=c is the power
# series
> Sum((D@@n)(f)(c)/n!*(x-c)^n, n=0..infinity);

infinity
-----    (n)              n
\      D   (f)(c) (x - c)
)     -------------------
/               n!
-----
n = 0
# If f is analytic at x=c, this series converges to f(x) for x
# in some interval around x=c.
# Maple has the "taylor" command to find a given number of terms
# of the Taylor series of an expression.
> S:=taylor(exp(x), x=2, 5);

2                     3
S := exp(2) + exp(2) (x - 2) + 1/2 exp(2) (x - 2)  + 1/6 exp(2) (x - 2)

4            5
+ 1/24 exp(2) (x - 2)  + O((x - 2) )
# The "5" tells Maple to what order to compute the series, i.e.
# up to but not including the term in (x-2)^5.
# The result of "taylor" looks like an ordinary polynomial with
# O((x-2)^5), but it's really something a bit different: a
# special "series" data structure.
> whattype(S);

series
# For just about any manipulations that you want to do with a
# series, you'll need the Taylor polynomial rather than this
# "series" structure.  You can't just subtract the O((x-2)^5),
# you must use "convert":
> Sp:= convert(S, polynom);

2                     3
Sp := exp(2) + exp(2) (x - 2) + 1/2 exp(2) (x - 2)  + 1/6 exp(2) (x - 2)

4
+ 1/24 exp(2) (x - 2)
# Let's look at the Taylor series for exp about x=0
# and its convergence to exp.  It's convenient to define a
# function that will calculate the n'th degree Taylor polynomial
# at a given point.
> P:= (n,t) -> subs(x=t, convert(taylor(exp(x), x=0, n+1), polynom));

P := (n,t) -> subs(x = t, convert(taylor(exp(x), x = 0, n + 1), polynom))
> P(4,t);

2        3         4
1 + t + 1/2 t  + 1/6 t  + 1/24 t
> plot({ exp(x), seq(P(nn, x), nn=1 .. 10) },
>       x=-6 .. 2, y=-2 .. 8, colour = black);

** Maple V Graphics **

# For x > 0, it's more convenient to look at the remainders.
> plot({ seq(exp(x)-P(nn, x), nn=1 .. 12) },
>     x = -6 .. 6, y = -3 .. 3, colour=black);

** Maple V Graphics **

# An animation is another possibility.  There is an
# "animate" command, but here it's simplest to use
# "display" from the "plots" package, with the option
# "insequence=true".
> with(plots,display):
> display([ seq(
>   plot(exp(x)-P(nn,x), x=-7 .. 7, y=-3 .. 3),
>       nn=1 .. 16) ], insequence=true);
# It's almost (but not quite) true that the curves for x > 0
# march off to the right at a constant rate of 1/E per step.
> plot({ seq((exp(x+nn/E) - P(nn, x+nn/E)), nn = 1 .. 20) },\
x=-1 .. 1.5, y = -1 .. 2, color=black);

** Maple V Graphics **

# The y value here turns out to be approximately proportional to
# 1/sqrt(n).  Here is an animation where we multiply by
# that factor sqrt(n).
> display([ seq(\
plot((exp(x+nn/E) - P(nn, x+nn/E))*sqrt(nn), x=-1 .. 1.5, y = -1 .. 2), nn=1..25) ], insequence=true);
# This can be analyzed as follows.  We are looking at
> R:= (n,x) -> Sum((x+n/E)^k/k!, k=n+1..infinity);

infinity
-----            k
\      (x + n/E)
R := (n,x) ->     )     ----------
/          k!
-----
k = n + 1
# At least for x > -n/E, the terms are all positive, so R(n,x)
# is greater than the first term, which is
> FirstTerm:= (x+n/E)^(n+1)/(n+1)!;

(n + 1)
(x + n/E)
FirstTerm := ----------------
(n + 1)!
# On the other hand, for x < n/E, each successive term is
# multiplied by (x+n/E)/k < 2/E, so
> 'R(n,x)' < (x+n/E)^(n+1)/(n+1)!*  sum((2/E)^k, k=0 .. infinity);

(n + 1)
(x + n/E)        E
R(n, x) < ------------------
(n + 1)! (- 2 + E)
# Thus, up to a factor that's between 1 and E/(E-2), R(n,x)
# is the first term.  As n -> infinity, Stirling's formula
# says n! is asymptotically
> asympt(n!, n,2);

1/2   1/2
2    Pi            1/2   1/2      1/2          3/2
---------- + 1/12 2    Pi    (1/n)    + O((1/n)   )
1/2
(1/n)
---------------------------------------------------
n
exp(n) (1/n)
> asympt(FirstTerm,n,3);

1/2                     1/2
exp(x E - 1) 2    exp(-2) exp(2) (1/n)
1/2 ----------------------------------------- + O(1/n)
1/2
Pi
# So let's look at R(n,x) divided by this.  The result should
# be approximately constant (or at least between 1 and E/(E-2)).
> display([ seq(\
plot((exp(x+nn/E) - P(nn, x+nn/E))*sqrt(2*nn*Pi)*exp(1-x*E), x=-1 .. 1.5, y = -1 .. 2), nn=1 .. 25) ], insequence=true);
#
# We'd get a very different picture for, e.g., arctan(x).  The
# Taylor series in this case:
> P:= subs(exp=arctan, eval(P));

P := (n,t) -> subs(x = t, convert(taylor(arctan(x), x = 0, n + 1), polynom))
> P(10,t);

3        5        7        9
t - 1/3 t  + 1/5 t  - 1/7 t  + 1/9 t
> plot({ seq(arctan(x)-P(2*nn-1, x), nn=1 .. 10) }, x=-2..2, \
y=-1..1, colour=black);

** Maple V Graphics **

# Here the Taylor series only converges on the interval
# [-1 .. 1], and outside that interval P(n,x) is useless as an
# approximation to arctan(x).
#
# Various operations can be done to obtain new series from old
# series:
# the basic operations of arithmetic, as well as substitution,
# differentiation and integration.
#
# Example: Starting with "well-known" series, obtain the degree
# 10 Taylor polynomial for ln(1+x^2) sin(cos(x))) in powers of x.
#
#
> taylor(1/(1+t), t=0);

2    3    4    5      6
1 - t + t  - t  + t  - t  + O(t )
> int(", t);

2        3        4        5        6      7
t - 1/2 t  + 1/3 t  - 1/4 t  + 1/5 t  - 1/6 t  + O(t )
# This is the Taylor series for ln(1+t) in powers of t.
> s1:= subs(t=x^2, convert(", polynom));

2        4        6        8        10        12
s1 := x  - 1/2 x  + 1/3 x  - 1/4 x  + 1/5 x   - 1/6 x
> s2:= taylor(cos(x), x=0, 11);

2         4          6            8              10      11
s2 := 1 - 1/2 x  + 1/24 x  - 1/720 x  + 1/40320 x  - 1/3628800 x   + O(x  )
# Now I want sin(cos(x)), which is approximately sin(s2).
# But it would be wrong to use the power series of sin x in
# powers of x: we need it in powers of x-1.
> s3:= taylor(sin(t), t=1, 11);

2                     3
s3 := sin(1) + cos(1) (t - 1) - 1/2 sin(1) (t - 1)  - 1/6 cos(1) (t - 1)

4                       5                       6
+ 1/24 sin(1) (t - 1)  + 1/120 cos(1) (t - 1)  - 1/720 sin(1) (t - 1)

7                         8
- 1/5040 cos(1) (t - 1)  + 1/40320 sin(1) (t - 1)

9                           10            11
+ 1/362880 cos(1) (t - 1)  - 1/3628800 sin(1) (t - 1)   + O((t - 1)  )
> subs(t=convert(s2, polynom), convert(s3,polynom));

2                3                 4
sin(1) + cos(1) %1 - 1/2 sin(1) %1  - 1/6 cos(1) %1  + 1/24 sin(1) %1

5                  6                   7
+ 1/120 cos(1) %1  - 1/720 sin(1) %1  - 1/5040 cos(1) %1

8                     9                      10
+ 1/40320 sin(1) %1  + 1/362880 cos(1) %1  - 1/3628800 sin(1) %1

2         4          6            8              10
%1 :=     - 1/2 x  + 1/24 x  - 1/720 x  + 1/40320 x  - 1/3628800 x
> taylor(", x=0, 11);

2                               4
sin(1) - 1/2 cos(1) x  + (1/24 cos(1) - 1/8 sin(1)) x

6   /   209                       \  8
+ (7/360 cos(1) + 1/48 sin(1)) x  + |- ----- cos(1) + 1/960 sin(1)| x
\  40320                      /

/  1259             193        \  10      12
+ |------- cos(1) - ------ sin(1)| x   + O(x  )
\3628800          241920       /
> taylor(s1*convert(", polynom), x=0, 11);

2                                4                                6
sin(1) x  + (- 1/2 cos(1) - 1/2 sin(1)) x  + (7/24 cos(1) + 5/24 sin(1)) x

/  121                    \  8   /143           4999       \  10      12
+ |- --- cos(1) - 1/6 sin(1)| x  + |--- sin(1) + ----- cos(1)| x   + O(x  )
\  720                    /      \960          40320       /
# Find the first six terms of the Taylor series for f(x) in
# powers of x-1, if y=f(x) satisfies
# the equation x^4 + y^4 = 2 x y with y(1) = 1.
#
# We expect the answer to be 1 plus a sum of coefficients times
# powers of x-1.
> ys:= 1 + sum(a[n] * (x-1)^n, n=1..5);

ys :=

2               3               4               5
1 + a[1] (x - 1) + a[2] (x - 1)  + a[3] (x - 1)  + a[4] (x - 1)  + a[5] (x - 1)
# Plug that in to the equation.
> x^4 + ys^4 = 2 * x * ys;

4
x  + (

2               3               4               5
1 + a[1] (x - 1) + a[2] (x - 1)  + a[3] (x - 1)  + a[4] (x - 1)  + a[5] (x - 1)

)^4 = 2 x (

2               3               4               5
1 + a[1] (x - 1) + a[2] (x - 1)  + a[3] (x - 1)  + a[4] (x - 1)  + a[5] (x - 1)

)
# Take the difference between the two sides, and take the Taylor
# series (discarding terms in (x-1)^n for n > 5).
>  g:= taylor(op(1,") - op(2,"), x=1, 6);

2                               2
g := (2 a[1] + 2) (x - 1) + (6 a[1]  + 2 a[2] + 6 - 2 a[1]) (x - 1)

2                                                        3
+ (4 + 4 (a[1]  + 2 a[2]) a[1] + 4 a[2] a[1] + 2 a[3] - 2 a[2]) (x - 1)  +

2
(1 + 4 a[3] a[1] + 2 a[2]  + 2 a[4] + 4 (2 a[2] a[1] + 2 a[3]) a[1]

2          2                  4
+ (a[1]  + 2 a[2])  - 2 a[3]) (x - 1)  + (2 a[5] + 4 a[1] a[4]

2
+ 4 a[2] a[3] + 4 a[1] (2 a[3] a[1] + a[2]  + 2 a[4])

2                                                   5
+ 2 (a[1]  + 2 a[2]) (2 a[2] a[1] + 2 a[3]) - 2 a[4]) (x - 1)

6
+ O((x - 1) )
# The "coeff" command can be used to extract any coefficient.
>  coeff(g, x-1, 2);

2
6 a[1]  + 2 a[2] + 6 - 2 a[1]
# The coefficient of (x-1)^n for each n from 1 to 5 should give
# us an equation.
>  eqs:= { seq( coeff(g, x-1, nn) = 0, nn=1..5) };

2
eqs := {6 a[1]  + 2 a[2] + 6 - 2 a[1] = 0,

2
1 + 4 a[3] a[1] + 2 a[2]  + 2 a[4] + 4 (2 a[2] a[1] + 2 a[3]) a[1]

2          2
+ (a[1]  + 2 a[2])  - 2 a[3] = 0,

2
2 a[5] + 4 a[1] a[4] + 4 a[2] a[3] + 4 a[1] (2 a[3] a[1] + a[2]  + 2 a[4])

2
+ 2 (a[1]  + 2 a[2]) (2 a[2] a[1] + 2 a[3]) - 2 a[4] = 0,

2
4 + 4 (a[1]  + 2 a[2]) a[1] + 4 a[2] a[1] + 2 a[3] - 2 a[2] = 0,

2 a[1] + 2 = 0}
# Now solve this set of equations for the variables a[1] to a[5].
>  solve(eqs, {a[1], a[2], a[3], a[4], a[5]});

{a[1] = -1, a[2] = -7, a[3] = -49, a[5] = -4627, a[4] = -449}
# So here is our solution:
> subs(", ys);

2             3              4               5
2 - x - 7 (x - 1)  - 49 (x - 1)  - 449 (x - 1)  - 4627 (x - 1)
>

```