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 
# much simpler answer.
> 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)
>