# Lesson 24. Power series, continued
# -------------------------------------------
> restart;
# We saw last time how to determine the Taylor series for
# an implicit function. Actually Maple's "taylor" can
# calculate a series for a "RootOf" expression directly.
> eqn:= x^4 + y^4 = 2*x*y;
4 4
eqn := x + y = 2 x y
> taylor(RootOf(eqn,y), x = 1, 5);
2 3 4 5
1 - (x - 1) - 7 (x - 1) - 49 (x - 1) - 449 (x - 1) + O((x - 1) )
# But there's something strange here. We haven't used the
# condition y(1)=1. There are actually two real y values
# that satisfy the equation when x=1.
> with(plots):\
display({implicitplot(eqn, x=-1.5 .. 1.5, y=-1.5 .. 1.5),\
plot([[1,0],[1,1.5]])});
** Maple V Graphics **
# How did Maple know that it should take one root of the
# equation rather than the other?
# The answer is that it didn't know! It chose one root, which
# happened to be the one we wanted. In another case it might
# not have chosen the "right" one. How can we influence its
# choice?
#
# One way is an additional input to "RootOf" that gives the
# approximate value of the root in question.
> taylor(RootOf(eqn,y,1), x = 1, 5);
2 3 4 5
1 - (x - 1) - 7 (x - 1) - 49 (x - 1) - 449 (x - 1) + O((x - 1) )
> taylor(RootOf(eqn,y,1/2), x = 1, 5);
/ 17 2 12 \ /467 300 2 538\ 2
%1 + |---- %1 + 8/11 %1 + ----| (x - 1) + |--- %1 + --- %1 + ---| (x - 1)
\ 11 11 / \121 121 121/
/33771 21796 2 40262\ 3
+ |----- %1 + ----- %1 + -----| (x - 1)
\ 1331 1331 1331/
/4064781 2210234 2 3410300 \ 4 5
+ |------- + ------- %1 + ------- %1| (x - 1) + O((x - 1) )
\ 14641 14641 14641 /
3 2
%1 := RootOf(_Z + _Z + _Z - 1, 1/2)
> evalf(");
2
.5436890127 + 2.146135923 (x - 1.) + 7.277537948 (x - 1.)
3 4 5
+ 48.88487606 (x - 1.) + 448.8944622 (x - 1.) + O((x - 1.) )
# What is the radius of convergence of our series (let's take the
# first one)? The coefficients seem to grow rapidly.
> ts:= taylor(RootOf(eqn,y,1),x = 1, 20);
2 3 4 5
ts := 1 - (x - 1) - 7 (x - 1) - 49 (x - 1) - 449 (x - 1) - 4627 (x - 1)
6 7 8
- 51199 (x - 1) - 594085 (x - 1) - 14264587/2 (x - 1)
9 10 11
- 175702481/2 (x - 1) - 2207953761/2 (x - 1) - 28195399867/2 (x - 1)
12 13
- 182414675533 (x - 1) - 2386447732823 (x - 1)
14 15
- 31513572925963 (x - 1) - 419492527683129 (x - 1)
16 17
- 22492227943132537/4 (x - 1) - 303339206031843395/4 (x - 1)
18 19
- 4113049885482030151/4 (x - 1) - 56037827987699400917/4 (x - 1)
20
+ O((x - 1) )
> convert(ts, polynom);
2 3 4 5 6
2 - x - 7 (x - 1) - 49 (x - 1) - 449 (x - 1) - 4627 (x - 1) - 51199 (x - 1)
7 8 9
- 594085 (x - 1) - 14264587/2 (x - 1) - 175702481/2 (x - 1)
10 11 12
- 2207953761/2 (x - 1) - 28195399867/2 (x - 1) - 182414675533 (x - 1)
13 14
- 2386447732823 (x - 1) - 31513572925963 (x - 1)
15 16
- 419492527683129 (x - 1) - 22492227943132537/4 (x - 1)
17 18
- 303339206031843395/4 (x - 1) - 4113049885482030151/4 (x - 1)
19
- 56037827987699400917/4 (x - 1)
> subs(x=t+1,");
2 3 4 5 6 7 8
1 - t - 7 t - 49 t - 449 t - 4627 t - 51199 t - 594085 t - 14264587/2 t
9 10 11 12
- 175702481/2 t - 2207953761/2 t - 28195399867/2 t - 182414675533 t
13 14 15
- 2386447732823 t - 31513572925963 t - 419492527683129 t
16 17
- 22492227943132537/4 t - 303339206031843395/4 t
18 19
- 4113049885482030151/4 t - 56037827987699400917/4 t
> L:= [coeffs(")];
L := [1, -1, -7, -49, -449, -4627, -51199, -594085, -14264587/2, -175702481/2,
-2207953761/2, -28195399867/2, -182414675533, -2386447732823,
-31513572925963, -419492527683129, -22492227943132537/4,
-303339206031843395/4, -4113049885482030151/4, -56037827987699400917/4]
> ratios:= [seq(L[jj]/L[jj+1], jj=1..nops(L)-1)];
49 449 4627 51199 1188170 14264587
ratios := [-1, 1/7, 1/7, ---, ----, -----, ------, --------, ---------,
449 4627 51199 594085 14264587 175702481
175702481 2207953761 28195399867 182414675533 2386447732823
----------, -----------, ------------, -------------, --------------,
2207953761 28195399867 364829351066 2386447732823 31513572925963
31513572925963 1677970110732516 22492227943132537
---------------, -----------------, ------------------,
419492527683129 22492227943132537 303339206031843395
303339206031843395 4113049885482030151
-------------------, --------------------]
4113049885482030151 56037827987699400917
> evalf(");
[-1., .1428571429, .1428571429, .1091314031, .09703911822, .09037285884,
.08618127036, .08329508594, .08118603061, .07957706547, .07830900684,
.07728380347, .07643774176, .07572761548, .07512308527, .07460221882,
.07414876645, .07375043203, .07339773923]
# The limit of these ratios, if it exists, would be the radius
# of convergence. It does appear plausible that the limit
# exists and is about .07.
> plot([seq([nn, ratios[nn]], nn=4..19)]); \
** Maple V Graphics **
# Note what happens on the implicit plot for x about 1.07.
# Let's find this point more exactly.
> subs(x=x(y), eqn);
4 4
x(y) + y = 2 x(y) y
> diff(", y);
3 / d \ 3 / d \
4 x(y) |---- x(y)| + 4 y = 2 |---- x(y)| y + 2 x(y)
\ dy / \ dy /
> solve(", diff(x(y),y));
3
4 y - 2 x(y)
- -------------
3
4 x(y) - 2 y
> solve({ eqn, 4*y^3-2*x = 0}, {x,y});
8 8 3
{y = 0, x = 0}, {y = RootOf(16 _Z - 3), x = 2 RootOf(16 _Z - 3) }
> xmax:= evalf(2* (3/16)^(3/8));
xmax := 1.067592398
# Taylor series can also be used for many other purposes.
# For example, integration. What is the arc length of the
# ellipse x^2 + y^2/b^2 = 1?
> ye:= b*sqrt(1-x^2);
2 1/2
ye := b (1 - x )
> L:= 4*Int(sqrt(diff(ye,x)^2 + 1), x=0..1);
1
/ / 2 2 \1/2
| | b x |
L := 4 | |------ + 1| dx
| | 2 |
/ \1 - x /
0
# Maple can't express the integral as a closed form. But it
# can write the integrand as a Taylor series, and integrate
# each term of the series. Various series are possible.
# But it's best to first express the integral as a proper one.
> student[changevar](x=sin(t),L,t);
1/2 Pi
/ / 2 2 \1/2
| | b sin(t) |
4 | |- ------------- + 1| cos(t) dt
| | 2 |
/ \ - 1 + sin(t) /
0
> simplify(");
1/2 Pi
/ / 2 2 2 2\1/2
| | - b + b cos(t) - cos(t) |
4 | |- ---------------------------| cos(t) dt
| | 2 |
/ \ cos(t) /
0
> simplify(",symbolic);
1/2 Pi
/
| 2 2 2 2 1/2
4 I | (- b + b cos(t) - cos(t) ) dt
|
/
0
> L:= 4*Int(sqrt(b^2-b^2*cos(t)^2+cos(t)^2), t=0 .. Pi/2);
1/2 Pi
/
| 2 2 2 2 1/2
L := 4 | (b - b cos(t) + cos(t) ) dt
|
/
0
> value(");
1/2 Pi
/
| 2 2 2 2 1/2
4 | (b - b cos(t) + cos(t) ) dt
|
/
0
> taylor(L,b=1);
2 3 17 4
2 Pi + Pi (b - 1) + 1/8 Pi (b - 1) - 1/16 Pi (b - 1) + --- Pi (b - 1)
512
19 5 6
- ---- Pi (b - 1) + O((b - 1) )
1024
# Note that b = 1 makes the ellipse into a circle.
