# 1.17 Lesson 17

```# Lesson 17.  Integration in Elementary Functions
# -----------------------------------------------
#
# The main idea behind the integration of rational functions, which we saw
# last time, is the following:
#
# 1) Predict the general form that the antiderivative of the given function f
# will take.
# 2) See what equations between the coefficients are required to make F' = f.
# 3) Solve these equations.
#
# In the case of p/q where p and q are polynomials in x, degree(p) <
# degree(q), and q has roots r[i] with multiplicities n[i], the general form
# was the sum of
# a[0,i] ln(x - r[i]) and a[j,i]/(x-r[i])^j for 1 <= j <= n[i]-1.  The
# equation F' = f is the partial fraction decomposition of f.
#
# The same general idea occurs repeatedly in mathematics.  In differential
# equations it is the method of Undetermined Coefficients.  Staying with
# integration, we could use it to find the antiderivative of a function such
# as the following:
> restart: \
f:= x^3 * exp(x) * cos(2*x);

3
f := x  exp(x) cos(2 x)
# Suppose we know (or can guess) that the antiderivative is a sum of terms
# x^i exp(x) cos(2 x) and x^i exp(x) sin(2 x), where i goes from 0 to 3.
> F:= sum((a[i] * cos(2*x) + b[i] * sin(2*x))* x^i * exp(x), i = 0 .. 3);

F := (a[0] cos(2 x) + b[0] sin(2 x)) exp(x)

+ (a[1] cos(2 x) + b[1] sin(2 x)) x exp(x)

2
+ (a[2] cos(2 x) + b[2] sin(2 x)) x  exp(x)

3
+ (a[3] cos(2 x) + b[3] sin(2 x)) x  exp(x)
> d:= normal(diff(F,x) - f);

d := - 2 exp(x) a[0] sin(2 x) + 2 exp(x) b[0] cos(2 x) + exp(x) a[0] cos(2 x)

+ exp(x) b[0] sin(2 x) - 2 x exp(x) a[1] sin(2 x)

+ 2 x exp(x) b[1] cos(2 x) + exp(x) a[1] cos(2 x) + exp(x) b[1] sin(2 x)

+ x exp(x) a[1] cos(2 x) + x exp(x) b[1] sin(2 x)

2                           2
- 2 x  exp(x) a[2] sin(2 x) + 2 x  exp(x) b[2] cos(2 x)

+ 2 x exp(x) a[2] cos(2 x) + 2 x exp(x) b[2] sin(2 x)

2                         2
+ x  exp(x) a[2] cos(2 x) + x  exp(x) b[2] sin(2 x)

3                           3
- 2 x  exp(x) a[3] sin(2 x) + 2 x  exp(x) b[3] cos(2 x)

2                           2
+ 3 x  exp(x) a[3] cos(2 x) + 3 x  exp(x) b[3] sin(2 x)

3                         3                         3
+ x  exp(x) a[3] cos(2 x) + x  exp(x) b[3] sin(2 x) - x  exp(x) cos(2 x)
> coeff(d,x,3);
Error, unable to compute coeff

# That was because d was not a polynomial.  First let's divide out the exp(x).
> d:= normal(d/exp(x));

d := - 2 a[0] sin(2 x) + 2 b[0] cos(2 x) + a[0] cos(2 x) + b[0] sin(2 x)

- 2 x a[1] sin(2 x) + 2 x b[1] cos(2 x) + a[1] cos(2 x) + b[1] sin(2 x)

2
+ x a[1] cos(2 x) + x b[1] sin(2 x) - 2 x  a[2] sin(2 x)

2
+ 2 x  b[2] cos(2 x) + 2 x a[2] cos(2 x) + 2 x b[2] sin(2 x)

2                  2                    3
+ x  a[2] cos(2 x) + x  b[2] sin(2 x) - 2 x  a[3] sin(2 x)

3                    2                    2
+ 2 x  b[3] cos(2 x) + 3 x  a[3] cos(2 x) + 3 x  b[3] sin(2 x)

3                  3                  3
+ x  a[3] cos(2 x) + x  b[3] sin(2 x) - x  cos(2 x)
# Now separate the terms in cos(2 x) and sin(2 x).
> costerms:= subs(cos(2*x)=1, sin(2*x)=0, d);

2
costerms := 2 b[0] + a[0] + 2 x b[1] + a[1] + x a[1] + 2 x  b[2] + 2 x a[2]

2           3           2         3         3
+ x  a[2] + 2 x  b[3] + 3 x  a[3] + x  a[3] - x
> sinterms:= subs(cos(2*x)=0, sin(2*x)=1, d);

2
sinterms := - 2 a[0] + b[0] - 2 x a[1] + b[1] + x b[1] - 2 x  a[2] + 2 x b[2]

2           3           2         3
+ x  b[2] - 2 x  a[3] + 3 x  b[3] + x  b[3]
> eqns:= { seq(coeff(costerms, x, ii), ii = 0 .. 3),\
seq(coeff(sinterms, x, ii), ii = 0 .. 3) };

eqns := {2 b[3] + a[3] - 1, - 2 a[1] + b[1] + 2 b[2], - 2 a[3] + b[3],

2 b[0] + a[0] + a[1], 2 b[1] + a[1] + 2 a[2], 2 b[2] + a[2] + 3 a[3],

- 2 a[0] + b[0] + b[1], - 2 a[2] + b[2] + 3 b[3]}
> solve(eqns, {seq(a[ii], ii = 0 .. 3), seq(b[ii], ii= 0 .. 3)});

144          42            12            66            12
{b[0] = ---, a[0] = ---, b[1] = - ---, a[1] = - ---, b[2] = - ----,
625         625           125           125            25

a[2] = 9/25, a[3] = 1/5, b[3] = 2/5}
> subs(",F);

/ 42            144         \          /   66             12         \
|--- cos(2 x) + --- sin(2 x)| exp(x) + |- --- cos(2 x) - --- sin(2 x)| x exp(x)
\625            625         /          \  125            125         /

/                 12          \  2
+ |9/25 cos(2 x) - ---- sin(2 x)| x  exp(x)
\                 25          /

3
+ (1/5 cos(2 x) + 2/5 sin(2 x)) x  exp(x)
> int(f, x);

/     3         2    66      42\
|1/5 x  + 9/25 x  - --- x + ---| exp(x) cos(2 x)
\                   125     625/

/       3    12   2    12     144\
- |- 2/5 x  + ---- x  + --- x - ---| exp(x) sin(2 x)
\            25       125     625/
# The Risch-Norman integration algorithm is also based on the same idea, using
# a fundamental principle of Liouville that restricts what functions can appear
# in an elementary antiderivative of a function.  It would be difficult to
# present the complete principle here, so we'll look at one special case,
# dealing with functions of the form g exp(h), where g and h are rational
# functions and h is not
# constant.  For these, Liouville's principle says any elementary
# antiderivative must be of the form H exp(h) where H is a rational function.
# Any root of the denominator of H must be a root of the denominator of g
# and/or of h.  In particular, if g and h are polynomials so is H.  Now that
# still leaves lots of possibilities, but we can find restrictions on the
# degrees of the numerator and the denominator.
# For example:
> f:= (2*x^2+2*x+1)*exp(x^2);

2                 2
f := (2 x  + 2 x + 1) exp(x )
> F:= H(x)*exp(x^2);

2
F := H(x) exp(x )
> dF:= diff(F, x);

/  d      \      2                  2
dF := |---- H(x)| exp(x ) + 2 H(x) x exp(x )
\ dx      /
# Since g = 2x^2 + 2 x + 1 and h = x^ are polynomials, H must be a polynomial.
#  What degree can it have?  If the highest power of x in H is x^n, then the
# highest power in df is x^(n+1).  We don't want any terms higher than x^2, so
# we conclude n=1.
> H:= x -> a[1]*x + a[0];

H := x -> a[1] x + a[0]
> normal((dF - f)/exp(x^2));

2                      2
a[1] + 2 x  a[1] + 2 x a[0] - 2 x  - 2 x - 1
> seq(coeff(",x,jj)=0, jj=0 .. 2);

a[1] - 1 = 0, 2 a[0] - 2 = 0, 2 a[1] - 2 = 0
> solve({"}, {a[0], a[1]});

{a[0] = 1, a[1] = 1}
> subs(", F);

2
(x + 1) exp(x )
# We were "lucky": there were three equations in two unknowns, but they had a
# solution.  If there would not have been a solution, we would have concluded
# that there was no elementary antiderivative.
# The best-known example of an elementary function with no elementary
# antiderivative is exp(-x^2).
> H:= 'H': f:= exp(-x^2): F:= H(x)*exp(-x^2): dF:=diff(F,x);

/  d      \        2                    2
dF := |---- H(x)| exp(- x ) - 2 H(x) x exp(- x )
\ dx      /
# H must be a polynomial.  If the highest power of x in H(x) is x^n, then the
# highest power in dF is x^(n+1).  But the highest power must be 0, which is
# impossible.
#
# Now a harder example:
# For what values of b and c does (x^2 + b x + c)/(x^2 (x-1)^2) exp(1/x) have
# an elementary antiderivative?
> H:= 'H': f:= (x^2+b*x+c)/(x^2 *(x-1)^2)*exp(1/x); F:= H(x)*exp(1/x):\
dF:= diff(F, x);

2
(x  + b x + c) exp(1/x)
f := -----------------------
2        2
x  (x - 1)

/  d      \            H(x) exp(1/x)
dF := |---- H(x)| exp(1/x) - -------------
\ dx      /                   2
x
# The only possible roots of the denominator of H are 0 and 1.  If that
# denominator has an x^n:
> eval(subs(H(x)=p(x)/x^n, dF));

/  d               \
|---- p(x)         |
| dx         p(x) n|            p(x) exp(1/x)
|--------- - ------| exp(1/x) - -------------
|     n        n   |                 n  2
\    x        x  x /                x  x
# The denominator of f would have x^(n+2).  It has x^2, so n=0, i.e.
# denominator of H contains no power of x.  What about (x-1)^n?
> eval(subs(H(x)=p(x)/(x-1)^n, dF));

/  d                         \
|---- p(x)                   |
| dx              p(x) n     |            p(x) exp(1/x)
|--------- - ----------------| exp(1/x) - -------------
|        n          n        |                    n  2
\ (x - 1)    (x - 1)  (x - 1)/             (x - 1)  x
# Unless n=0, the denominator has (x-1)^(n+1).  It actually has (x-1)^2, so n
# = 1.
> eval(subs(H(x)=p(x)/(x-1), dF));

/  d                 \
|---- p(x)           |
| dx           p(x)  |            p(x) exp(1/x)
|--------- - --------| exp(1/x) - -------------
|  x - 1            2|                       2
\            (x - 1) /              (x - 1) x
> dFp:= normal(");

/ 3 /  d      \    2 /  d      \    2                     \
exp(1/x) |x  |---- p(x)| - x  |---- p(x)| - x  p(x) - p(x) x + p(x)|
\   \ dx      /      \ dx      /                          /
dFp := --------------------------------------------------------------------
2        2
x  (x - 1)
# If the  highest power of x in p(x) is x^m, then what's the highest power in
# the numerator here?
> eval(subs(p(x)= x^m, dFp));

2  m        m      2  m    m      m
exp(1/x) (x  x  m - x x  m - x  x  - x  x + x )
-----------------------------------------------
2        2
x  (x - 1)
# There's (m-1) x^(m+2) here.  Unless m=1, the highest power is x^(m+2).  If
# m=1 ...
> subs(m=1,");

2
exp(1/x) (- 2 x  + x)
---------------------
2        2
x  (x - 1)
# ... the highest power is x^2.  We want the highest power to be x^2, so m = 1.
> H:= x -> (a[1]*x + a[0])/(x-1);

a[1] x + a[0]
H := x -> -------------
x - 1
> normal((dF-f)*(x-1)^2*x^2/exp(1/x));

2         2                                  2
- 2 x  a[1] - x  a[0] + x a[1] - x a[0] + a[0] - x  - b x - c
> eqns:= { seq(coeff(", x, jj), jj = 0 .. 2)} ;

eqns := {a[1] - a[0] - b, - 2 a[1] - a[0] - 1, a[0] - c}
> solve(",{a[0], a[1],b});

{a[0] = c, b = - 3/2 c - 1/2, a[1] = - 1/2 c - 1/2}
# Conclusion: c is arbitrary, with b = -(3 c + 1)/2.  Compare to the result of
# "int".
> int(f,x);

exp(1/x)                           /  exp(1/x)                          \
-------- + exp(1) Ei(1, 1 - 1/x) - |- -------- - 2 exp(1) Ei(1, 1 - 1/x)| b
1/x - 1                           \   1/x - 1                          /

/           exp(1/x)                          \
- |exp(1/x) - -------- - 3 exp(1) Ei(1, 1 - 1/x)| c
\            1/x - 1                          /
# The "Ei" function is not elementary.  The condition b = -(3 c+1)/2 is just
# what we need to make the Ei terms cancel.
> normal(subs(b=-(3*c+1)/2, "));

exp(1/x) (x + x c - 2 c)
- 1/2 ------------------------
x - 1
>

```