%% Below, I've used Maple, treating polynomials as /expressions/. I
%% could, instead, have treated polynomials as functions.
{{{Maple}}} Wed Feb 10 21:56:13 EST 2010
|\^/| Maple 10 (SUN SPARC SOLARIS)
._|\| |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2005
\ MAPLE / All rights reserved. Maple is a trademark of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
%% Here is our numerator.
> Numxm :=
10 9 8 7 6 5 4 3 2
29 + x + x + 7 x + 2 x + 13 x - 17 x + 5 x - 34 x + 14 x - 12 x ;
%% And here is our denominator.
8 7 6 5 4 3 2
> Dxm := 8 + x + x + 4 x - 2 x + x - 11 x + 2 x - 4 x ;
%% Let's compute the quotient and remainder polynomials:
> Quox := quo(Numxm, Dxm, x, 'Remx'); Remx := Remx ;
2
Quox := x + 3
7 3
Remx := x + 3 x + 5
%% Since I know how to integrate a poly, Quox, we'll now deal with
%% integrating the rational fnc
%% Remx / Dxm .
%% The form I want to write Remx/Dxm in is:
> parfracexp := A1/(x-1) + A2/(x-1)^2 + ((B1*x)+C1)/(x^2 + x + 2) + ((B2*x)+C2)/(x^2 + x + 2)^2 + ((B3*x)+C3)/(x^2 + x + 2)^3 ;
A1 A2 B1 x + C1 B2 x + C2 B3 x + C3
parfracexp := ----- + -------- + ---------- + ------------- + -------------
x - 1 2 2 2 2 2 3
(x - 1) x + x + 2 (x + x + 2) (x + x + 2)
%% So we have a set, {A1,A2, B1,C1, B2,C2, B3,C3}, of eight unknown numbers.
%%
%% In order to solve for them, note that the difference
%% [parfracexp * Dxm] - Remx
%% must be the zero-poly; which is what I called "zip" in class and in the
%% "Primer on Polynomials".
%% So we will simply set each coeff of [parfracexp*Dxm] - Remx equal to
%% zero, then solve the resulting system of eight eqns; this will give us
%% values for {A1,A2, B1,C1, B2,C2, B3,C3}.
> Diffx := collect(simplify((parfracexp * Dxm) - Remx), x);
7 6 5
Diffx := (B1 + A1 - 1) x + (A2 + 2 A1 + C1) x + (2 B1 + 6 A1 + 3 A2 + B2) x
4
+ (4 A1 - 4 B1 + 2 C1 + C2 + 9 A2 - B2) x
3
+ (-4 C1 + 13 A2 + B3 - 3 + B1 + 5 A1 - C2 + B2) x
2
+ (-6 A1 + C3 + 18 A2 - 4 B1 + C1 - 3 B2 - 2 B3 + C2) x
+ (-3 C2 - 4 C1 + 12 A2 + 2 B2 - 2 C3 + B3 - 4 A1 + 4 B1) x - 8 A1 + 4 C1
+ 8 A2 - 5 + C3 + 2 C2
%% The following step is a mathematically-unnecessary step... but, for some
%% reason, the Maple `collect' cmd does not collect the constant term as a
%% multiplier of x^0, So, below, I multiply Diffx by x, and collect, so that
%% NOW, the eight coeffs of x^8, x^7, ..., x^1 are the coeffs I want.
> DiffxTimesx := collect(simplify(x*Diffx), x) ;
8 7
DiffxTimesx := (B1 + A1 - 1) x + (A2 + 2 A1 + C1) x
6 5
+ (2 B1 + 6 A1 + 3 A2 + B2) x + (4 A1 - 4 B1 + 2 C1 + C2 + 9 A2 - B2) x
4
+ (-4 C1 + 13 A2 + B3 - 3 + B1 + 5 A1 - C2 + B2) x
3
+ (-6 A1 + C3 + 18 A2 - 4 B1 + C1 - 3 B2 - 2 B3 + C2) x
2
+ (-3 C2 - 4 C1 + 12 A2 + 2 B2 - 2 C3 + B3 - 4 A1 + 4 B1) x
+ (-8 A1 + 4 C1 + 8 A2 - 5 + C3 + 2 C2) x
%% So there, above, are the eight coeffs I want to set to zero and then solve
%% for the eight unknows.
%% The following Maple snippet will extract the coeffs, create eqns setting
%% each coeff equal to zero, and put them in a Maple "set", which is indicated by
%% the braces.
> Eqns := { seq( op(1, op(n,DiffxTimesx)) = 0, n=1..8 ) } ;
Eqns := {B1 + A1 - 1 = 0, A2 + 2 A1 + C1 = 0, 2 B1 + 6 A1 + 3 A2 + B2 = 0,
4 A1 - 4 B1 + 2 C1 + C2 + 9 A2 - B2 = 0,
-4 C1 + 13 A2 + B3 - 3 + B1 + 5 A1 - C2 + B2 = 0,
-6 A1 + C3 + 18 A2 - 4 B1 + C1 - 3 B2 - 2 B3 + C2 = 0,
-3 C2 - 4 C1 + 12 A2 + 2 B2 - 2 C3 + B3 - 4 A1 + 4 B1 = 0,
-8 A1 + 4 C1 + 8 A2 - 5 + C3 + 2 C2 = 0}
%% Now we solve this system of linear eqns. We have eight equations in eight
%% unknowns, and it will have a unique soln.
> solve( Eqns, {A1,A2, B1,C1, B2,C2, B3,C3} );
273 59 -69 -1
{ B1 = ---, A2 = 9/64, B3 = --, C3 = 9/8, C2 = 9/8, B2 = ---, C1 = ---,
256 16 32 128
-17
A1 = --- }
256
%% As a partial check, let's ask Maple to directly compute the partial-fraction
%% decomposition of Remx/Dxm.
> convert( Remx / Dxm , parfrac, x) ;
18 + 59 x -2 + 273 x 9 17
---------------- + ---------------- + ----------- - -----------
2 3 2 2 256 (x - 1)
16 (x + x + 2) 256 (x + x + 2) 64 (x - 1)
36 - 69 x
+ ----------------
2 2
32 (x + x + 2)
%%%%%%% Let's see if our `convert' produced the same answer %%%%%%%%%%%
The immediately-above expression indeed gives the same values as what
solving "by hand" produced. F'irinstance, the expression
17
- -----------
256 (x - 1)
from the convert(...) can be rewritten as
-17 / 256
---------- .
(x - 1)
And indeed, the A1 that we solved for, up above, is
-17
A1 = --- .
256
%% Pretty Nifty, eh? %%