[FIXED] Partial fraction decomposition using SymPy / Python

Issue

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

1/(x**6+1) = (A*x+B)/(x**2+1) + (C*x+D)/(x**2-sqrt(3)*x+1) + (K*x+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: 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:

A = 0
B = 1/3
C = -1/(2*sqrt(3))
D = 1/3
K = 1/(2*sqrt(3))
S = 1/3

How can I verify if these are correct?

Edit 2

The main point of my question is: how to do this with some nice/reusable Python code?

Solution

You can do this using apart in sympy but apart will look for a rational factorisation by default so you have to tell it to work in Q(sqrt(3)):

In [37]: apart(1/(x**6+1))                                                                                                                     
Out[37]: 
        2                     
       x  - 2           1     
- ─────────────── + ──────────
    ⎛ 4    2    ⎞     ⎛ 2    ⎞
  3⋅⎝x  - x  + 1⎠   3⋅⎝x  + 1⎠

In [36]: apart(1/(x**6+1), extension=sqrt(3))                                                                                                  
Out[36]: 
       √3⋅x - 2            √3⋅x + 2           1     
- ───────────────── + ───────────────── + ──────────
    ⎛ 2           ⎞     ⎛ 2           ⎞     ⎛ 2    ⎞
  6⋅⎝x  - √3⋅x + 1⎠   6⋅⎝x  + √3⋅x + 1⎠   3⋅⎝x  + 1⎠

EDIT: The comments are asking for a way to find this in more general cases without needing to know that sqrt(3) generates the extension.

We can use .apart(full=True) to compute the full PFE over linear factors:

In [29]: e = apart(1/(x**6+1), full=True).doit()

In [30]: e
Out[30]: 
       √3   ⅈ           √3   ⅈ          √3   ⅈ           √3   ⅈ                            
     - ── - ─         - ── + ─          ── - ─           ── + ─                            
       2    2           2    2          2    2           2    2           ⅈ           ⅈ    
- ────────────── - ────────────── - ────────────── - ────────────── + ───────── - ─────────
    ⎛    √3   ⅈ⎞     ⎛    √3   ⅈ⎞     ⎛    √3   ⅈ⎞     ⎛    √3   ⅈ⎞   6⋅(x + ⅈ)   6⋅(x - ⅈ)
  6⋅⎜x + ── + ─⎟   6⋅⎜x + ── - ─⎟   6⋅⎜x - ── + ─⎟   6⋅⎜x - ── - ─⎟                        
    ⎝    2    2⎠     ⎝    2    2⎠     ⎝    2    2⎠     ⎝    2    2⎠ 

This is not what is wanted because you want to have real quadratic denominators rather than introducing complex numbers. The terms here come in complex conjugate pairs though so we can combine the pairs:

In [46]: terms1 = []

In [47]: for a in terms:
    ...:     if conjugate(a) not in terms1:
    ...:         terms1.append(a)
    ...: 

In [49]: terms_real = [(t+t.conjugate()) for t in terms1]

In [51]: Add(*(factor(cancel(t)) for t in terms_real))
Out[51]: 
      √3⋅x - 2            √3⋅x + 2            1     
- ───────────────── + ───────────────── + ──────────
    ⎛ 2           ⎞     ⎛ 2           ⎞     ⎛ 2    ⎞
  6⋅⎝x  - √3⋅x + 1⎠   6⋅⎝x  + √3⋅x + 1⎠   3⋅⎝x  + 1⎠

Note that in general there might not be any simple expressions for the roots (Abel-Ruffini) so this kind of expression for the partial fraction expansion will not succeed for all possible denominator polynomials. This is why .apart by default computes an expansion over irreducible denominators (something that can always succeed).

Answered By - Oscar Benjamin