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