# Lesson 27. Asymptotic Series
# -----------------------------
> restart;
# Today we look at series that approximate a quantity as a variable (say n) goes to +infinity.
# In some cases the "asympt" command does this for us.
> q:= 1/(n^2+1) - n/(n^3+3);
1 n
q := ------ - ------
2 3
n + 1 n + 3
> asympt(q,n);
1 3 1
- ---- + ---- + O(----)
4 5 6
n n n
# What is "asympt" doing? Change variables to t=1/n which goes to 0 as n -> infinity.
> subs(n=1/t, q);
1 1
-------- - ------------
1 / 1 \
---- + 1 t |---- + 3|
2 | 3 |
t \ t /
> qt:= normal(");
4
t (3 t - 1)
qt := -------------------
2 3
(1 + t ) (1 + 3 t )
# This has a Taylor series in powers of t.
> taylor(",t);
4 5 6
- t + 3 t + O(t )
# Now change variables back to n.
> subs(t=1/n,");
1 3 1
- ---- + ---- + O(----)
4 5 6
n n n
# 1/q also has an asymptotic expression.
> asympt(1/q, n);
4 3 2
- n - 3 n + O(n )
> 1/qt;
2 3
(1 + t ) (1 + 3 t )
-------------------
4
t (3 t - 1)
> taylor(",t);
Error, does not have a taylor expansion, try series()
# It doesn't have a Taylor series because it has a singularity at t=0 (look at the t^4 in the
# denominator). But if we multiplied the expression by t^4 we would get a series.
> taylor(t^4/qt, t);
2 3 4 5 6
- 1 - 3 t - 10 t - 33 t - 99 t - 300 t + O(t )
> convert(",polynom)/t^4;
2 3 4 5
- 1 - 3 t - 10 t - 33 t - 99 t - 300 t
------------------------------------------
4
t
> expand(subs(t=1/n,"));
4 3 2 300
- n - 3 n - 10 n - 33 n - 99 - ---
n
# But not all asymptotic expressions arise this way.
> int(exp(-t)/t, t=x .. infinity);
Ei(1, x)
> asympt(", x, 10);
1 2 6 24 120 720 5040 40320 1
1/x - ---- + ---- - ---- + ---- - --- + --- - ---- + ----- + O(---)
2 3 4 5 6 7 8 9 10
x x x x x x x x x
-------------------------------------------------------------------
exp(x)
# Let's see if we can reproduce this series.
> J:= exp(x)*Int(exp(-t)/t, t=x .. infinity);
infinity
/
| exp(- t)
J := exp(x) | -------- dt
| t
/
x
> with(student):
> changevar(t=x+u,J,u);
infinity
/
| exp(- x - u)
exp(x) | ------------ du
| x + u
/
0
> J:= expand(");
infinity
/
| 1
J := | -------------- du
| exp(u) (x + u)
/
0
> series(1/(x+u),u);
1 1 2 1 3 1 4 1 5 6
1/x - ---- u + ---- u - ---- u + ---- u - ---- u + O(u )
2 3 4 5 6
x x x x x
# This series is valid when |u| < x. It isn't always so in
# this integral. But let's not worry about rigour.
> subs(1/(x+u)=convert(",polynom),J);
2 3 4 5
u u u u u
infinity 1/x - ---- + ---- - ---- + ---- - ----
/ 2 3 4 5 6
| x x x x x
| -------------------------------------- du
| exp(u)
/
0
> value(");
5 4 3 2
x - x + 24 x - 120 + 2 x - 6 x
----------------------------------
6
x
> expand(");
1 24 120 2 6
1/x - ---- + ---- - --- + ---- - ----
2 5 6 3 4
x x x x x
# Let's take the sum all the way to infinity.
> 1/(x+u)=Sum((-u)^k/x^(k+1), k=0..infinity);
infinity
----- k
1 \ (- u)
----- = ) --------
x + u / (k + 1)
----- x
k = 0
# A useful integral:
> assume(m>0): int(u^m*exp(-u), u=0..infinity);
GAMMA(m~) m~
# This is equal to m! because GAMMA(m) = (m-1)!.
> J = Sum((-1)^k * k!/x^(k+1), k=0..infinity);
infinity infinity
/ ----- k
| 1 \ (-1) k!
| -------------- du = ) --------
| exp(u) (x + u) / (k + 1)
/ ----- x
0 k = 0
# Do you believe this series?
# For what x does it converge?
# Nevertheless, it is useful!
> J:= intparts(J,(x+u)^(-1));
infinity
/
| 1
J := 1/x - | --------------- du
| 2
/ (x + u) exp(u)
0
> J:= intparts(J,(x+u)^(-2));
infinity
/
1 | 2
J := 1/x - ---- + | --------------- du
2 | 3
x / (x + u) exp(u)
0
> for jj from 3 to 8 do J:= intparts(J,(x+u)^(-jj)) od;
infinity
/
1 2 | 6
J := 1/x - ---- + ---- - | --------------- du
2 3 | 4
x x / (x + u) exp(u)
0
infinity
/
1 2 6 | 24
J := 1/x - ---- + ---- - ---- + | --------------- du
2 3 4 | 5
x x x / (x + u) exp(u)
0
infinity
/
1 2 6 24 | 120
J := 1/x - ---- + ---- - ---- + ---- - | --------------- du
2 3 4 5 | 6
x x x x / (x + u) exp(u)
0
infinity
/
1 2 6 24 120 | 720
J := 1/x - ---- + ---- - ---- + ---- - --- + | --------------- du
2 3 4 5 6 | 7
x x x x x / (x + u) exp(u)
0
infinity
/
1 2 6 24 120 720 | 5040
J := 1/x - ---- + ---- - ---- + ---- - --- + --- - | --------------- du
2 3 4 5 6 7 | 8
x x x x x x / (x + u) exp(u)
0
1 2 6 24 120 720 5040
J := 1/x - ---- + ---- - ---- + ---- - --- + --- - ----
2 3 4 5 6 7 8
x x x x x x x
infinity
/
| 40320
+ | --------------- du
| 9
/ (x + u) exp(u)
0
> for nn from 1 to 5 do S[nn]:= sum(op(k,J), k=1..nn) od;
S[1] := 1/x
1
S[2] := 1/x - ----
2
x
1 2
S[3] := 1/x - ---- + ----
2 3
x x
1 2 6
S[4] := 1/x - ---- + ---- - ----
2 3 4
x x x
1 2 6 24
S[5] := 1/x - ---- + ---- - ---- + ----
2 3 4 5
x x x x
> for nn from 6 to 30 do S[nn]:= sum((-1)^k*k!/x^(k+1), k=0..nn) od:
> evalf(subs(x=1,
> [exp(x)*Ei(1,x), seq(S[nn],nn=1..8)]));
[.5963473622, 1., 0, 2., -4., 20., 620., -4420., 35900.]
> evalf(subs(x=5,
> [exp(x)*Ei(1,x), seq(S[nn],nn=1..20)]));
[.1704221762, .2000000000, .1600000000, .1760000000, .1664000000, .1740800000,
.1756160000, .1627136000, .1833574400, .1461985280, .2205163520,
.05701713920, .4494152499, -.5708198380, 2.285838408, -6.284136330,
21.13978283, -72.10154232, 263.5672282, -1011.974100, 4090.191212]
> evalf(subs(x=10,
> [exp(x)*Ei(1,x), seq(S[nn],nn=1..30)]));
[.09156333393, .1000000000, .09000000000, .09200000000, .09140000000,
.09164000000, .09159200000, .09154160000, .09158192000, .09154563200,
.09158192000, .09154200320, .09158990336, .09152763315, .09161481144,
.09148404401, .09169327191, .09133758448, .09197782185, .09076137084,
.09319427285, .08808517863, .09932518591, .07347316917, .1355180093,
-.01959409109, .3836973700, -.7051895750, 2.343693871, -6.498068123,
20.02721786]
# Let f(x) be a smooth function (as many derivatives as we want)
# and F(x) an "antidifference", i.e.
> eq1:= f(x) = F(x+1) - F(x):
> map(t -> sum(subs(x=j, t), j=1..5), eq1) ;
f(1) + f(2) + f(3) + f(4) + f(5) = - F(1) + F(6)
> F(n+1) = F(1) + Sum(f(j), j=1 .. n);
/ n \
|----- |
| \ |
F(n + 1) = F(1) + | ) f(j)|
| / |
|----- |
\j = 1 /
> F(x+1) = exp(D)(F)(x);
F(x + 1) = exp(D)(F)(x)
> f(x) = exp(D)(F)(x) - F(x);
f(x) = exp(D)(F)(x) - F(x)
> taylor(exp(t)-1,t);
2 3 4 5 6
t + 1/2 t + 1/6 t + 1/24 t + 1/120 t + O(t )
> convert(",polynom);
2 3 4 5
t + 1/2 t + 1/6 t + 1/24 t + 1/120 t
> f(x) = D(F)(x) + 1/2*(D@@2)(F)(x) + 1/6*(D@@3)(F)(x) + 1/24*(D@@4)(F)(x) + 1/120*(D@@5)(F)(x) + err;
(2) (3) (4)
f(x) = D(F)(x) + 1/2 D (F)(x) + 1/6 D (F)(x) + 1/24 D (F)(x)
(5)
+ 1/120 D (F)(x) + err
> Q:= m -> sum((D@@(k+m))(F)(n)/k!, k = 1 .. 7-m);
7 - m
----- (k + m)
\ D (F)(n)
Q := m -> ) --------------
/ k!
-----
k = 1
> f(n) = Q(0);
(2) (3) (4)
f(n) = D(F)(n) + 1/2 D (F)(n) + 1/6 D (F)(n) + 1/24 D (F)(n)
(5) (6) (7)
+ 1/120 D (F)(n) + 1/720 D (F)(n) + 1/5040 D (F)(n)
> D(f)(n) = Q(1);
(2) (3) (4) (5)
D(f)(n) = D (F)(n) + 1/2 D (F)(n) + 1/6 D (F)(n) + 1/24 D (F)(n)
(6) (7)
+ 1/120 D (F)(n) + 1/720 D (F)(n)
> int(f(x),x=0..n) = Q(-1)+K;
n
/
| (2) (3)
| f(x) dx = F(n) + 1/2 D(F)(n) + 1/6 D (F)(n) + 1/24 D (F)(n)
|
/
0
(4) (5) (6)
+ 1/120 D (F)(n) + 1/720 D (F)(n) + 1/5040 D (F)(n)
(7)
+ 1/40320 D (F)(n) + K
> taylor(u/(exp(u)-1), u,9);
2 4 6 8
1 - 1/2 u + 1/12 u - 1/720 u + 1/30240 u + O(u )
> int(f(x),x=0..n) - 1/2*f(n) + 1/\
12*D(f)(n) - 1/720*(D@@3)(f)(n) + 1/30240*(D@@5)(f)(n)=\
Q(-1)+K - 1/2*Q(0) + 1/12*Q(1) - 1/720 * Q(3) + 1/30240*Q(5);
n
/
| (3) (5)
| f(x) dx - 1/2 f(n) + 1/12 D(f)(n) - 1/720 D (f)(n) + 1/30240 D (f)(n) =
|
/
0
K + F(n)
K + F(n)
