$1/(x^6 + 1)$ partial fraction decomposition, with computer?

Question

This is a re-post from StackOverflow, I was advised to post it here.

https://stackoverflow.com/questions/64101194/partial-fraction-decomposition

How do I find the constants A,B,C,D,K,S such that

$$ \frac{1}{x^6+1} = \frac{Ax+B}{x^2+1} + \frac{Cx+D}{x^2-\sqrt 3 x+1} + \frac{Kx+S}{x^2+\sqrt 3 x+1} $$

is true for every real x.

I need some sympy code maybe, not sure. Or... any other Python lib which could help here.

I tried by hand but it's not easy at all... and after 1 hour of calculating, I found that I have probably made some mistake.

I tried partial fraction decomposition in SymPy but it does not go that far.

I tried Wolfram Alpha too, but it also does not decompose to that level of detail, it seems.

WA attempt

See the alternate forms which WA gives below.

EDIT: I did a second try entirely by hand and I got these:

\begin{align}A &= 0,\\ B &= \frac13,\\ C &= -\frac1{2\sqrt3},\\ D &= \frac13,\\ K &= \frac1{2\sqrt3},\\ S &= \frac13. \end{align}

Could someone verify if these are correct?
And in general... how can I automate this task via SymPy or WA?

multiply through by the product of the denominators (you can be smarter and not multiply by common factors). Now you have a linear equation in the space of polynomials; use some matrix solve() — Calvin Khor, Sep 28 '20 at 12:51
@CalvinKhor That's the thing, I know how to do this manually but not by code. So how do I do this by code? What if I need to decompose something slightly more complicated e.g. $(x^3-1)/(x^{10} + 1)$? — peter.petrov, Sep 28 '20 at 12:53
This particular rational function is handled here. In general it is probably easier to handle $1/(x^m+1)$ using complex number and roots of unity. You need to combine complex conjugate terms to get real coefficients (needed for calculus style indefinite integrals and such). — Jyrki Lahtonen, Sep 28 '20 at 12:54
@peter.petrov I believe my comment is perfectly servicable pseudocode? Haven't bothered to code it up so maybe I'm dumb — Calvin Khor, Sep 28 '20 at 12:54
@JyrkiLahtonen It's handed but noone gave a decent final answer, it seems :) — peter.petrov, Sep 28 '20 at 12:55
Sorry about jumping the gun a bit there. Anyway, the method in Paul Garret's answer is how I would approach this (and combine the complex conjugate pairs). — Jyrki Lahtonen, Sep 28 '20 at 12:57
In your WA attempt, Apart supports Factor with Extension, for example: Partial fraction expansion with complex numbers. And btw, instead of writing the WA Mathematica code into WA, you can write it into an online wolfram mathematica notebook and then click run (shift + enter). — Vepir, Oct 02 '20 at 16:35

Alessio K · Answer 1 · 2020-09-28T12:59:58.793

What you got is indeed correct since we have $$ \frac{1}{3}\frac{1}{x^2+1} + \frac{-\frac{x}{2\sqrt{3}}+\frac{1}{3}}{x^2-\sqrt 3 x+1} + \frac{\frac{x}{2\sqrt{3}}+\frac{1}{3}}{x^2+\sqrt 3 x+1} $$ $$=\frac{1}{3}\frac{1}{x^2+1}+\frac{2-x^2}{3(x^2-\sqrt{3}x+1)(x^2+\sqrt{3}+1)}$$ $$=\frac{1}{3}\frac{1}{x^2+1}+\frac{2-x^2}{3(x^4-x^2+1)}$$ $$=\frac{1}{3}\big[\frac{x^4-x^2+1+(2-x^2)(x^2+1)}{(x^2+1)(x^4-x^2+1)}\big]$$ $$=\frac{1}{3}\frac{3}{x^6+1}=\frac{1}{x^6+1}.$$

WA also agrees.

score 1 · Answer 2 · answered Sep 28 '20 at 14:01

You're trying to solve for coefficients $a_0,\dots a_{2M-1}$ in $$ \frac1P = \sum_{i=0}^M \frac{a_{2i} x + a_{2i+1}}{Q_i}$$

where $Q_i$ are 2nd order. and $P$ is some high order $K$ polynomial. Presumably $Q_i$ are factors of $P$. A naive algorithm:

let let $S_i = P/Q_i = \sum_{k=0}^{K-2} s_{i,k} x^k$. Also choose $s_{i,-1}:=0$ and $s_{i,K-2+l}:=0$ for $l>0$, and maybe other things zero, you'll find out in debug.
Multiply by $P$ to get \begin{align} 1 &= \sum_{i=0}^{M} S_i (a_i x + b_i) \\ &= \sum_{i=0}^M \sum_{k=0}^{K-2} a_{2i}s_{i,k}x^{k+1} + a_{2i+1}s_{i,k}x^{k} \\ &= \sum_{k=0}^{K-1}\sum_{i=0}^M (a_{2i}s_{i,k-1} + a_{2i+1}s_{i,k})x^k\\ &= \sum_{i=0}^M (a_{2i}s_{i,-1} + a_{2i+1}s_{i,0}) + \sum_{i=0}^M(a_{2i}s_{i,0} + a_{2i+1}s_{i,1})x + \sum_{i=0}^M(a_{2i}s_{i,1} + a_{2i+1}s_{i,2})x^2 + \dots \end{align}
Recognize this as the matrix equation $$ \begin{bmatrix}1\\0\\ \vdots \\ 0\end{bmatrix} = \begin{bmatrix} s_{0,-1} & s_{0,0} & s_{1,-1} & s_{1,0} & s_{2,-1} & s_{2,0} & \dots \\ s_{0,0} & s_{0,1} & s_{1,0}& s_{1,1}& s_{2,0} & s_{2,1}& \dots \\ s_{0,1} & s_{0,2} & s_{1,1}& s_{1,2}& s_{2,1} & s_{2,2}& \dots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix}a_0\\ a_1\\ \vdots \\ a_{2M-1}\end{bmatrix}$$
Find a solution with your matrix solver.

score 1 · Answer 3 · answered Oct 02 '20 at 15:35

A very simple approach would be to just substitute in 6 values for $x$ to get a linear system of equations. For example, substituting in $x\in\{-2\sqrt3,-\sqrt3,0,\sqrt3,2\sqrt3,3\sqrt3\}$ gives the following equations

$$\begin{bmatrix}-2/13&1/13&-2/19&1/19&-2/7&1/7\\-1/4&1/4&-1/7&1/7&-1&1\\0&1&0&1&0&1\\1/4&1/4&1&1&1/7&1/7\\2/13&1/13&2/7&1/7&2/19&1/19\\3/28&1/28&3/19&1/19&3/37&1/37\end{bmatrix}\begin{bmatrix}A\sqrt3\\B\\C\sqrt3\\D\\K\sqrt3\\S\end{bmatrix}=\begin{bmatrix}1/1729\\1/28\\1\\1/28\\1/1729\\1/19683\end{bmatrix}$$

which is very easy to solve using most things.

Note: This assumes the given form is correct, so it must work on any 6 chosen points. Just because it works for 6 points does not prove it will work for all points. In some cases you may also notice the matrix on the LHS to be singular. If there are no solutions, the chosen form is wrong. If there are multiple solutions, different points should be used.

Yeah, basically that's how I solved it by hand. Thanks. – peter.petrov Oct 03 '20 at 10:49 — peter.petrov, Oct 03 '20 at 10:49

$1/(x^6 + 1)$ partial fraction decomposition, with computer?

3 Answers3