2

I need to find the coefficient of a term in a rather long, nasty expansion. I have a polynomial, say f(x) = (x+x^2)/2 and then a function that is defined recursively: g_k(x,y) = y*f(g_{k-1}(x,y)) with g_0(x,y)=yx.

I want to know, say, the coefficient of x^2y^4 in g_10(x,y)

I've coded this up as

import sympy
x, y = sympy.symbols('x y')

def f(x):
      return (x+x**2)/2
def g(x,y,k):
     if k==0:
          return y*x
     else:
          return y*f(g(x,y,k-1))

fxn = g(x,y,2)     
fxn.expand().coeff(x**2).coeff(y**4)
> 1/4

So far so good.

But now I want to find a coefficient for k = 10. Now fxn = g(x,y,10) and then fxn.expand() is very slow. Obviously there are a lot of steps going on, so it's not a surprise. But my knowledge of sympy is rudimentary - I've only started using it specifically because I need to be able to find these coefficients. I could imagine that there may be a way to get sympy to recognize that everything is a polynomial and so it can more quickly find a particular coefficient, but I haven't been able to find examples doing that.

Is there another approach through sympy to get this coefficient, or anything I can do to speed it up?

Joel
  • 22,598
  • 6
  • 69
  • 93
  • The answer for `k=10` is `0`. Looking closely at your problem you will see that the lowest order of `y` is `k+1`. Hence for `k=1` there will be no coefficient `y**4`. Or rather the coefficient will be `0`. – laolux Aug 07 '17 at 18:28
  • For that particular coefficient, yes. I'll be looking for other coefficients as well. Your answer looks to be well on the way to what I'm after. Will wait to see if other suggestions come up before accepting. – Joel Aug 07 '17 at 20:09

1 Answers1

1

I assume you are only interested in the coefficients given and not the whole polynomial g(x,y,10). So you can redefine your function g to get rid of higher orders in every step of the recursion. This will significantly speed up your calculation.

def g(x,y,k):
    if k==0:
        return y*x
    else:
        temp = y*f(g(x,y,k-1)) + sympy.O(y**5) + sympy.O(x**3)
        return temp.expand().removeO()

Works as follows: First everything of the order O(y**5), O(x**3) (and higher) will be grouped and then discarded. Keep in mind you loose lots of information!

Also have a look here: Sympy: Drop higher order terms in polynomial

laolux
  • 1,445
  • 1
  • 17
  • 28
  • This answer addresses much of my situation. An issue that comes up for me is that my `f(x)` will potentially have a constant term. So it could be `f(x) = (1+x+x^2)/3`. The `symp.O(y**5)` type approach would still work if I know i'm after a power of `y` that is small enough, but the constant term would prevent me from being able to throw out an `x**3` term since `f(x)**3` would have terms of order 0, 1, and 2. Any suggestions for this? – Joel Aug 07 '17 at 20:12
  • I don't get your problem with the constant term. If you are after `x**2` and your function does not include anything along `x**(-1)` then you can safely throw away coefficients the order of `O(x**3)`. The benefits will be less, but it's still valid. – laolux Aug 08 '17 at 12:48
  • Assuming f=(1+x+x^2)/3 and zeta_0 = yx gives: zeta_1 =yf(x) = y(1+x+x^2)/3 = y/3 + yx/3 + yx^2/3. Then zeta_2 = y/3 + yf(x)/3 +y f(x)^2/3, and you can see that there are O(1) and O(x) terms coming from the f(x)^2 term which is (1+x+x^2)^2/9. This will continue for any O(x^3) term which has (1+x+x^2)^3/27. – Joel Aug 08 '17 at 20:35
  • I am not sure what you try to tell me with zeta_n. Where does it come from? `g(x,y,1)=y*f(g(x,y,0))=y*f(y*x)` so clearly zeta_1 is not g(x,y,1). Anyways, it's not important for my argument. I am well aware of the fact that with every step your functions will create new expressions of the order of `O(x**3)`. But my point is that these terms (`O(x**3)`) will NEVER generate any terms of order `O(1), O(x), O(x**2)`. Thus `O(x**3)` can be safely discarded in any step. – laolux Aug 09 '17 at 07:54
  • There are many examples where I need to do this sort of thing. I gave one [mcve]. An O(X**3) will generate an O(1) term if we plug in a function X which itself has O(1) terms. Let's do an example with straight up iteration f(x) = (1+x+x^2), and I only want to know f(f(x)) to 1st order in x. Then f(f(x)) = 1+(1+x+x^2)+(1+x+x^2)^2 = 3 + 3x+ O(x^3). If I truncate f to hat{f} = 1+x, I get hat{f}(hat{f}(x)) = 2+x. – Joel Aug 09 '17 at 20:37
  • Ok, now I get your confusion. I don't ask you to truncate your function f, as this will clearly not work. Instead I ask you to truncate your results. Let's do an example. f(x) = (1+x+x^2). Then first iteration gives f1 = f(x). We now truncate f1 -> hat{f1}. Now we do the second iteration f2 = f(hat{f1}) = x^2 + 3*x + 3. Now truncate f2 -> hat{f2}. Then f3 = f(hat{f2}) = 9x^2 + 21*x + 13. For comparison: f(f(f(x))) = 13 + 21*x + 37*x^2 + O(x^3). So my suggestion does indeed work if properly applied to the intermediate results and not the function itself. This is exactly what my code does. – laolux Aug 10 '17 at 07:39