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