# 1.16 Lesson 16

```# Lesson 16.  Integration of Rational Functions
# ---------------------------------------------
#
# How does Maple find antiderivatives?
#
# In contrast to the completely mechanical procedure for
# differentiation, Calculus texts tend to have a hodgepodge of
# different techniques for integration, and often it is not
# clear which technique will work for a particular integral.  If
# we can't manage to make these techniques work, there's always
# the suspicion that if we were cleverer we could find an answer.
#
# There actually are algorithms that will find an elementary
# antiderivative if it exists, or else conclude that it does not
# exist.
# Some of the basic ingredients date back to Liouville in 1833,
# but the problem was not completely solved until quite
# recently: the last case of it only in 1987.  There still isn't
# a complete theory of integrals involving special functions
# (such as erf and dilog).  Parts of the theory involves some
# very high-powered mathematics, but I'll try to give you a
# taste of some of the simpler parts.  We'll start by looking at
# how to integrate a rational function.
> restart;

# Here is a rational function, i.e. a quotient of polynomials.
> p:= x^5 + 3*x^4 + 5:
> q:= x^4 -6*x^2-8*x-3:
> f:= p/q;

5      4
x  + 3 x  + 5
f := -------------------
4      2
x  - 6 x  - 8 x - 3
# What is its integral?
> int(f, x);

2         491                  7           21       107
1/2 x  + 3 x + --- ln(x - 3) + ---------- - ---------- - --- ln(x + 1)
64                      2   16 (x + 1)    64
8 (x + 1)
# How does the typical Calculus text tell you to integrate f?
# The first step is division: write p/q = s + r/q where s and r
# are
# polynomials, with degree(r) < degree(q).
> s:= quo(p,q,x); r:= rem(p,q,x);

s := x + 3

3       2
r := 14 + 6 x  + 26 x  + 27 x
> p/q = s + r/q;

5      4                          3       2
x  + 3 x  + 5              14 + 6 x  + 26 x  + 27 x
------------------- = x + 3 + ------------------------
4      2                         4      2
x  - 6 x  - 8 x - 3              x  - 6 x  - 8 x - 3
> normal(lhs(")-rhs("));

0
# Since integrating the polynomial s is easy, we're left with
# the
# problem of integrating r/q.
# The next step is to factor the denominator q.
> factor(q);

3
(x - 3) (x + 1)
# Now r/q is supposed to be decomposed into partial fractions: a
# sum of the following form:
> eqn:= r/q = a/(x-3)+ b/(x+1) + c/(x+1)^2 + d/(x+1)^3;

3       2
14 + 6 x  + 26 x  + 27 x     a       b         c          d
eqn := ------------------------ = ----- + ----- + -------- + --------
4      2               x - 3   x + 1          2          3
x  - 6 x  - 8 x - 3                     (x + 1)    (x + 1)
> normal(rhs(eqn)-lhs(eqn));

3      3      3      2        2      2       2
(- 6 x  + b x  + a x  + c x  + 3 a x  - b x  - 26 x  - 5 b x + d x - 27 x

/   4      2
+ 3 a x - 2 c x - 3 b - 3 c - 3 d - 14 + a)  /  (x  - 6 x  - 8 x - 3)
/
# The coefficient of each power of x in the numerator gives us
# an
# equation in the unknowns a,b,c and d, which we solve.
> eqns:= { seq(coeff(numer("),x,jj), jj=0 .. 3) };

eqns := {- 3 b - 3 c - 3 d - 14 + a, - 5 b + d - 27 + 3 a - 2 c, - 6 + b + a,

c + 3 a - b - 26}
> solve(eqns, {a,b,c,d});

491                 21         107
{a = ---, d = -7/4, c = ----, b = - ---}
64                 16          64
> subs(", rhs(eqn));

491          107           21            7
---------- - ---------- + ----------- - ----------
64 (x - 3)   64 (x + 1)             2            3
16 (x + 1)    4 (x + 1)
# And of course we know how to integrate each of these terms.
# The result (added to the integral of the quotient s) is our
# antiderivative.
> int(s + ", x);

2         491                  7           21       107
1/2 x  + 3 x + --- ln(x - 3) + ---------- - ---------- - --- ln(x + 1)
64                      2   16 (x + 1)    64
8 (x + 1)
# The main complication in general is that we might not be able
# to factor the denominator.  Let's see an example where the
# denominator doesn't come factored.  I'll keep the same
# numerator p.
#
> q:= x^7 + x^5 + x+1;

7    5
q := x  + x  + x + 1
> J:= int(p/q, x);

-----
\               140802562320802750657958612847777027   6
J :=    )    _R ln(x - ------------------------------------ _R
/                99160323514247198901748344164757701
-----
_R = %1

281486806421556593632406541713136531   5
+ ------------------------------------ _R
99160323514247198901748344164757701

102810458385112832371035624411126180   4
+ ------------------------------------ _R
99160323514247198901748344164757701

107090858261128471692239849657585760   3
- ------------------------------------ _R
99160323514247198901748344164757701

323125137891580778927623424433020655   2
- ------------------------------------ _R
99160323514247198901748344164757701

101549110751312422022027258214689296
+ ------------------------------------ _R
99160323514247198901748344164757701

94139296068596134296454605704199462
- -----------------------------------)
99160323514247198901748344164757701

7             5             4             3             2
%1 := RootOf(1651331 _Z  - 2239231 _Z  - 3108282 _Z  - 1549701 _Z  - 1422918 _Z

- 303508 _Z - 101227)
> ro:= %1;

7             5             4             3             2
ro := RootOf(1651331 _Z  - 2239231 _Z  - 3108282 _Z  - 1549701 _Z  - 1422918 _Z

- 303508 _Z - 101227)
> xfac:=op(op(1,J)/_R);

140802562320802750657958612847777027   6
xfac := x - ------------------------------------ _R
99160323514247198901748344164757701

281486806421556593632406541713136531   5
+ ------------------------------------ _R
99160323514247198901748344164757701

102810458385112832371035624411126180   4
+ ------------------------------------ _R
99160323514247198901748344164757701

107090858261128471692239849657585760   3
- ------------------------------------ _R
99160323514247198901748344164757701

323125137891580778927623424433020655   2
- ------------------------------------ _R
99160323514247198901748344164757701

101549110751312422022027258214689296
+ ------------------------------------ _R
99160323514247198901748344164757701

94139296068596134296454605704199462
- -----------------------------------
99160323514247198901748344164757701
# Curious.  Why the RootOf this strange-looking polynomial,
# rather than of q?  Well, let's at least look at numerical
# values.
> rts:= [allvalues(ro)];

rts := [ - .8015786748 - .5844355395 I,  - .8015786748 + .5844355395 I,

- .1149509834 - .2971155128 I,  - .1149509834 + .2971155128 I,

.05615476610 - .5945813620 I, .05615476610 + .5945813620 I, 1.720749784]
> map(r -> subs(_R=r, xfac), rts);

[x - .0842506825 - 1.236616987 I, x - .0842506825 + 1.236616987 I,

x - .8625699449 - .5517570073 I, x - .8625699449 + .5517570073 I,

x + .5889198574 - .7214431676 I, x + .5889198574 + .7214431676 I,

x + .7158015164]
# The way we might have done things in Math 101 is as follows
# (well, almost):

> qroots:= [allvalues(RootOf(q))];

qroots := [-.7158015389,  - .5889198579 - .7214431674 I,

- .5889198579 + .7214431674 I, .08425068236 - 1.236616986 I,

.08425068236 + 1.236616986 I, .8625699450 - .5517570070 I,

.8625699450 + .5517570070 I]
# There is one real root together with three pairs of complex
# roots.
# Note that (apart from roundoff error) the logarithmic terms in
# Maple's J are ln(x - qroots[i]).
# In Math 101 we would pair the linear factors x - qroots[i]
# corresponding to complex roots to get real quadratic roots.
# I won't do that here.
> decomp:= sum(a[i]/(x - qroots[i]), i = 1..7);

a[1]                      a[2]
decomp := --------------- + -------------------------------
x + .7158015389   x + .5889198579 + .7214431674 I

a[3]                              a[4]
+ ------------------------------- + --------------------------------
x + .5889198579 - .7214431674 I   x - .08425068236 + 1.236616986 I

a[5]                               a[6]
+ -------------------------------- + -------------------------------
x - .08425068236 - 1.236616986 I   x - .8625699450 + .5517570070 I

a[7]
+ -------------------------------
x - .8625699450 - .5517570070 I
> normal(decomp);

3                     4                     5
(.1263905265 a[2] x  + .8263463540 a[2] x  - .5889198570 a[2] x

2         6
+ .9098758523 a[2] x - .8655636316 a[2] x  + a[2] x  - .8932042823 I a[4] x

3                       5                       4
+ .6281112000 I a[4] x  + .7214431670 I a[3] x  - .8497444137 I a[3] x

3                       5
- .2978270646 I a[3] x + 1.096593291 I a[3] x  - .7214431674 I a[2] x

4                       3                       5
+ .8497444154 I a[2] x  + 1.615348609 I a[7] x  + .5517570070 I a[7] x

4                                              5
+ .9518580220 I a[7] x  + 1.392150280 I a[7] x - .5517570070 I a[6] x

4                       5                       4
- .9518580238 I a[6] x  + 1.236616986 I a[5] x  + .2083716491 I a[5] x

2                                              5
- .4259630752 I a[5] x  + .8932042828 I a[5] x - 1.236616986 I a[4] x

4                       2
- .2083716495 I a[4] x  + .4259630753 I a[4] x

-9         4                 -10         2
- .1000000000*10   I a[1] x  + .7000000000*10    I a[1] x

-10                  6                       2
- .9100000000*10    I a[1] x + a[3] x  - .5546219847 I a[3] x

2                     3
+ .9098758526 a[3] x - .8655636309 a[3] x  + .1263905263 a[3] x

4                       2                       3
+ .8263463540 a[3] x  + .5546219861 I a[2] x  - .6281112010 I a[5] x

3                       2
- 1.615348605 I a[6] x  - 1.788714683 I a[6] x  + .5900521839 a[4] x

2                     3                     4
+ .7513174834 a[4] x  - .3016651726 a[4] x  - .5221233932 a[4] x

5
+ .08425068300 a[4] x  - 1.392150280 I a[6] x + .2978270653 I a[2] x

6                       2
+ .8318199883 I a[3] + a[4] x  + 1.788714681 I a[7] x  - .5262533467 I a[6]

6
+ .8049216243 I a[5] - .8049216242 I a[4] + a[5] x  + .5900521845 a[5] x

2                     3                     4
+ .7513174846 a[5] x  - .3016651727 a[5] x  - .5221233934 a[5] x

5                              6
+ .08425068300 a[5] x  - .8318199886 I a[2] + a[6] x  - 1.222591894 a[6] x

2                     3                     4
- .2732022293 a[6] x  + .7165536965 a[6] x  + 1.439591115 a[6] x

5         6                                          2
+ .8625699460 a[6] x  + a[7] x  - 1.222591894 a[7] x - .2732022313 a[7] x

3                     4                     5
+ .7165536965 a[7] x  + 1.439591115 a[7] x  + .8625699460 a[7] x

3                     5
+ .5262533467 I a[7] - 1.096593294 I a[2] x  - .5889198570 a[3] x

- .8226996931 a[7] + .6790213445 a[2] + .6790213441 a[3]

- .05483928876 a[4] - .05483928876 a[5] - .8226996931 a[6]

6                                          2
+ 1.397035276 a[1] + a[1] x  - .5546722855 a[1] x + .7748967486 a[1] x

3                     4                     5    /
- 1.082558093 a[1] x  + 1.512371842 a[1] x  - .7158015380 a[1] x )  /  (
/

(x + .7158015389) (x + .5889198579 + .7214431674 I)

(x + .5889198579 - .7214431674 I) (x - .08425068236 + 1.236616986 I)

(x - .08425068236 - 1.236616986 I) (x - .8625699450 + .5517570070 I)

(x - .8625699450 - .5517570070 I))
> eqns:= { seq(coeff(numer(")-p,x,jj)=0, jj=0..6)};

eqns := { - .8655636316 a[2] - .4259630752 I a[5] + .4259630753 I a[4]

-10
+ .7000000000*10    I a[1] - .5546219847 I a[3] - .8655636309 a[3]

+ .5546219861 I a[2] - 1.788714683 I a[6] + .7513174834 a[4]

+ 1.788714681 I a[7] + .7513174846 a[5] - .2732022293 a[6]

- .2732022313 a[7] + .7748967486 a[1] = 0,

.1263905265 a[2] + .6281112000 I a[4] + 1.096593291 I a[3] + 1.615348609 I a[7]

+ .1263905263 a[3] - .6281112010 I a[5] - 1.615348605 I a[6]

- .3016651726 a[4] - .3016651727 a[5] + .7165536965 a[6] + .7165536965 a[7]

- 1.096593294 I a[2] - 1.082558093 a[1] = 0,

- .5889198570 a[2] - 1 + .7214431670 I a[3] - .7214431674 I a[2]

+ .5517570070 I a[7] - .5517570070 I a[6] + 1.236616986 I a[5]

- 1.236616986 I a[4] + .08425068300 a[4] + .08425068300 a[5]

+ .8625699460 a[6] + .8625699460 a[7] - .5889198570 a[3] - .7158015380 a[1]

= 0,

.8263463540 a[2] - 3 - .8497444137 I a[3] + .8497444154 I a[2]

+ .9518580220 I a[7] - .9518580238 I a[6] + .2083716491 I a[5]

-9
- .2083716495 I a[4] - .1000000000*10   I a[1] + .8263463540 a[3]

- .5221233932 a[4] - .5221233934 a[5] + 1.439591115 a[6] + 1.439591115 a[7]

+ 1.512371842 a[1] = 0,

a[2] + a[3] + a[4] + a[5] + a[6] + a[7] + a[1] = 0,

- 5 + .8318199883 I a[3] - .5262533467 I a[6] + .8049216243 I a[5]

- .8049216242 I a[4] - .8318199886 I a[2] + .5262533467 I a[7]

- .8226996931 a[7] + .6790213445 a[2] + .6790213441 a[3]

- .05483928876 a[4] - .05483928876 a[5] - .8226996931 a[6]

+ 1.397035276 a[1] = 0,

.9098758523 a[2] - .8932042823 I a[4] - .2978270646 I a[3] + 1.392150280 I a[7]

-10
+ .8932042828 I a[5] - .9100000000*10    I a[1] + .9098758526 a[3]

+ .5900521839 a[4] - 1.392150280 I a[6] + .2978270653 I a[2]

+ .5900521845 a[5] - 1.222591894 a[6] - 1.222591894 a[7] - .5546722855 a[1]

= 0

}
> solve(eqns);

{a[2] = .05615476706 + .5945813612 I, a[5] =  - .8015786752 - .5844355385 I,

a[4] =  - .8015786749 + .5844355409 I, a[3] = .0561547662 - .5945813626 I,

-9
a[6] =  - .1149509833 + .2971155127 I, a[1] = 1.720749783 - .20*10   I,

a[7] =  - .114950983 - .2971155135 I}
> subs(", decomp);

-9
1.720749783 - .20*10   I     .05615476706 + .5945813612 I
------------------------ + -------------------------------
x + .7158015389       x + .5889198579 + .7214431674 I

.0561547662 - .5945813626 I       - .8015786749 + .5844355409 I
+ ------------------------------- + --------------------------------
x + .5889198579 - .7214431674 I   x - .08425068236 + 1.236616986 I

- .8015786752 - .5844355385 I      - .1149509833 + .2971155127 I
+ -------------------------------- + -------------------------------
x - .08425068236 - 1.236616986 I   x - .8625699450 + .5517570070 I

- .114950983 - .2971155135 I
+ -------------------------------
x - .8625699450 - .5517570070 I
# Note that the coefficients here are (up to roundoff error) the
# same ones Maple found.
# There are some better ways to get the partial fraction
# decomposition than solving a large system of equations.
# Suppose a polynomial q(x) of degree n with leading coefficient
# 1 has no repeated roots.
# Then the coefficient a[i] is p(qroots[i])/q'(qroots[i]).
> pdq:= unapply(p/diff(q,x), x);

5      4
x  + 3 x  + 5
pdq := x -> ---------------
6      4
7 x  + 5 x  + 1
> seq(pdq(qroots[ii]), ii=1..7);

1.720749784, .05615476612 + .5945813617 I, .05615476612 - .5945813617 I,

- .8015786745 + .5844355392 I,  - .8015786745 - .5844355392 I,

- .1149509834 + .2971155127 I,  - .1149509834 - .2971155127 I
# The Rothstein-Trager method uses the fact that the integral
# can be expressed as a sum of terms b[j] ln(v[j]), where the
# b[j] are
# the roots of  resultant(p - z q', q, x) (as a polynomial in
# z), and v[j] is the greatest common divisor (gcd) of the
# polynomials p - b[j] q' and q.  This v[j] sometimes contains
# more than one factor x - x[i].
# So that's where the strange polynomial came from.
# Here's one where the Rothstein-Trager method does result in a
> resultant(p - z*diff(q,x), q, x);

2            3            4            5
101227 + 303508 z + 1422918 z  + 1549701 z  + 3108282 z  + 2239231 z

7
- 1651331 z
> p:= 4*x^4 - 33*x^2 - 12*x;
> q:= x^6 - 4*x^4 + 3*x^3 + 22*x^2 - 28;

4       2
p := 4 x  - 33 x  - 12 x

6      3      4       2
q := x  + 3 x  - 4 x  + 22 x  - 28
> resultant(p - z*diff(q,x), q, x);

6                  4                  2
- 470730907024 z  + 1412192721072 z  - 1412192721072 z  + 470730907024
> solve(");

1, -1, 1, -1, 1, -1
> p - diff(q,x);

4       2             5       3
4 x  - 42 x  - 56 x - 6 x  + 16 x
> gcd(", q);

3      2
x  - 2 x  + 7
> int(p/q, x);

3      2            3      2
ln(x  - 2 x  + 7) - ln(x  + 2 x  - 4)
>

```