Ask Your Question
1

find correct coefficients in $n$

asked 2020-02-25 15:32:04 -0600

brettpim gravatar image

updated 2020-02-26 08:19:39 -0600

def desc_prod_recip(x,y):
    t0 = 1
    t1 = y*(y-1)/(2*x)
    t2 = (y*(y+1)*(y-1)*(3*y-2))/(24*x^2)
    t3 = (y^2*(y-1)^2*(y+1)*(y+2))/(48*x^3)
    t4 = y*(y-1)*(y+1)*(y+2)*(y+3)*(15*y^3-15*y^2-10*y+8)/(5760*x^4)
    return (x^(-y))*(t0 + t1 + t2 + t3 + t4)
n,t = var('n t')
E1 = desc_prod_recip(n,t).expand().simplify()
E1.coefficients(n)

gives the answer:

 [[1/96*t^7 - 1/576*t^6, -t - 4],
 [-1/16*t^4, -t - 3],
 [1/(n^t), 0],
 [1/384*t^8 - 1/30*t^5 - 5/1152*t^4 + 1/32*t^3 + 1/288*t^2 - 1/120*t, -t- 4],
 [1/48*t^6 + 1/48*t^5 - 1/48*t^3 + 1/24*t^2, -t - 3],
 [1/8*t^4 - 1/12*t^3 - 1/8*t^2 + 1/12*t, -t - 2],
 [1/2*t^2 - 1/2*t, -t - 1]]

It seems to detect all the powers of $n$ properly except $n^{-t}$. How can I process this symbolic expression so that sagemath computes the coefficients of $n$ properly?

edit retag flag offensive close merge delete

Comments

I just realised that it also has two separate coefficients for $n^{-t-4}$ and $n^{-t-3}$, so there are other issues with .coefficients()

brettpim gravatar imagebrettpim ( 2020-02-26 17:59:22 -0600 )edit

1 answer

Sort by ยป oldest newest most voted
1

answered 2020-04-13 20:29:45 -0600

dan_fulea gravatar image

I have no idea why do we want to really do the above, and really in this way (so we lose control on the powers of $n$ while constructing an expression in $n$, then want to get back these powers...) but here is an answer. It is designed to work with the above example. Again, i would never do this. It is a better way to use or isolate the total factor $n^{-t}$.

First of all, let us initialize the E1 without a full simplification.

sage: def desc_prod_recip(x,y): 
....:     t0 = 1 
....:     t1 = y*(y-1)/(2*x) 
....:     t2 = (y*(y+1)*(y-1)*(3*y-2))/(24*x^2) 
....:     t3 = (y^2*(y-1)^2*(y+1)*(y+2))/(48*x^3) 
....:     t4 = y*(y-1)*(y+1)*(y+2)*(y+3)*(15*y^3-15*y^2-10*y+8)/(5760*x^4) 
....:     return (x^(-y))*(t0 + t1 + t2 + t3 + t4) 
....: n,t = var('n t'); 
....: E1 = desc_prod_recip(n,t)                                                                                                                                                                       

sage: E1 = E1.expand()                                                                                                                                                                                

sage: for operand in E1.operands(): 
....:     print(operand) 
....:                                                                                                                                                                                                 
1/384*t^8/(n^4*n^t)
1/48*t^6/(n^3*n^t)
1/96*t^7/(n^4*n^t)
1/8*t^4/(n^2*n^t)
1/48*t^5/(n^3*n^t)
-1/576*t^6/(n^4*n^t)
1/2*t^2/(n*n^t)
-1/12*t^3/(n^2*n^t)
-1/16*t^4/(n^3*n^t)
-1/30*t^5/(n^4*n^t)
1/(n^t)
-1/2*t/(n*n^t)
-1/8*t^2/(n^2*n^t)
-1/48*t^3/(n^3*n^t)
-5/1152*t^4/(n^4*n^t)
1/12*t/(n^2*n^t)
1/24*t^2/(n^3*n^t)
1/32*t^3/(n^4*n^t)
1/288*t^2/(n^4*n^t)
-1/120*t/(n^4*n^t)

Now we can put the hands on each operand. The code has the last one in its hands, and we can use for instance...

sage: operand.factor_list()                                                                                                                                                                           
[(n, -t - 4), (t, 1), (-1/120, 1)]

Now we can parse the "factor list" above with entries of the shape (the_factor, its_power), ask if the first component is n, if yes, record its_power, and then build the product of all not considered factors to the corresponding powers. This may be easily done, but unfortunately...

sage: 1/(n^t)                                                                                                                                                                                         
1/(n^t)
sage: _.factor_list()                                                                                                                                                                                 
[(n^t, -1)]

so the implementation needs to go recursively deeper. A first approximation but not the solution would be...

sage: ops = E1.operands()                                                                                                                                                                             
sage: n_powers = [ power for op in ops for base, power in op.factor_list() if base == n ]                                                                                                             
sage: n_powers = list(set(n_powers))                                                                                                                                                                  
sage: n_powers.sort()                                                                                                                                                                                 
sage: n_powers                                                                                                                                                                                        
[-t - 4, -t - 3, -t - 2, -t - 1]
sage: # so -t is missing here already, since 1/n^t is not correctly parsed...

And now we can build a dictionary with keys the powers of $n$, values their "coefficients".

sage: dic = {}                                                                                                                                                                                        
sage: for fixed_power in n_powers: 
....:     dic[fixed_power] = sum([ prod([ base^power  
....:                                     for base, power in op.factor_list() 
....:                                     if (base, power) != (n, fixed_power) ]) 
....:                              for op in ops 
....:                              if  (n, fixed_power) in op.factor_list()]) 
....:                                                                                                                                                                                                 
sage: dic                                                                                                                                                                                             
{-t - 4: 1/384*t^8 + 1/96*t^7 - 1/576*t^6 - 1/30*t^5 - 5/1152*t^4 + 1/32*t^3 + 1/288*t^2 - 1/120*t,
 -t - 3: 1/48*t^6 + 1/48*t^5 - 1/16*t^4 - 1/48*t^3 + 1/24*t^2,
 -t - 2: 1/8*t^4 - 1/12*t^3 - 1/8*t^2 + 1/12*t,
 -t - 1: 1/2*t^2 - 1/2*t}
sage:

This is only the half of the road, but we can still ask for the difference...

sage: diff = E1 - sum([ n^key*val for key, val in dic.items() ])                                                                                                                                      
sage: diff.simplify_full()                                                                                                                                                                            
1/(n^t)
sage:

As i said, i would never program such an adventure, instead change the method to attack the problem, the situation.

edit flag offensive delete link more

Comments

thank you @dan_fulea. I am definitely open to trying to do my calculation a different way. My example was a simplification of what I am really trying to do. Essentially I have a number of different expressions of forms similar to the above and I have to take some sums and products of them and then afterwards I need to extract the coefficients of particular powers of n. I was getting incorrect answers in some small examples and identified the issue that way. What is the method you would recommend for doing what I describe? Thank you for your time and expertise.

brettpim gravatar imagebrettpim ( 2020-04-20 17:54:21 -0600 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2020-02-25 15:32:04 -0600

Seen: 141 times

Last updated: Apr 13