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