# 1.27 Lesson 27

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

```