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