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