$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?

Answer 1

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.

Answer 2

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:

1. 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.

2. 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}

3. 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}$$

4. Find a solution with your matrix solver.

Answer 3

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.