The computation that we just carried out, gives the desired result, but it has a serious flaw: The equations were manually assembled. That is inconvenient and error-prone, and it does not lead to a useful generalisation. In this section we explore a better algorithm to set up the equations.
r: 1/(x^2*(x^2 + 1));
1 ----------- 2 2 x (x + 1)
Now, we enter the three elementary fractions that may occur in the partial fraction decomposition of that expression:
p1: a/x;
a
-
x
p2: b/x^2;
b
--
2
x
p3: (c*x + d)/(x^2 + 1);
c x + d
-------
2
x + 1
The partial fraction decomposition is the sum of these three fractions:
p1 + p2 + p3;
c x + d a b
------- + - + --
2 x 2
x + 1 x
To compare the unknown numerator coefficients with the numerator of the given expression, we have to rewrite that sum on a common denominator:
ratsimp(%);
3 2
(c + a) x + (d + b) x + a x + b
---------------------------------
4 2
x + x
Now, we access the numerator of this fraction:
n: num(%);
3 2
(c + a) x + (d + b) x + a x + b
This is a sum, but we need a list of its terms. To rewrite a sum as a list of terms, we have to use the function maplist:
terms: maplist (lambda ([item ], item), %);
3 2 [(c + a) x , (d + b) x , a x, b]
The following expression transforms the list of terms into a list of equations. For each term in terms, the term coefficient is equated to the coefficient of the corresponding term in r
map(lambda ([term], coeff(term, x, hipow(term, x)) = coeff(r, x, hipow(term, x))),
terms);
[c + a = 0, d + b = 0, a = 0, b = 1]
This ist is a system of linear equation. We can solve it with solve:
solve(%, [a, b, c, d]);
[[a = 0, b = 1, C = 0, d = - 1]]
at(p1 + p2 + p3, first(%));
1 1
-- - ------
2 2
x x + 1