2

Motivation. It is well known that generating function for Catalan numbers satisfies quadratic equation. I would like to have first several coefficients of a function, implicitly defined by an algebraic equation (not necessarily a quadratic one!).

Example. 

import sympy as sp
sp.init_printing() # math as latex
from IPython.display import display
z = sp.Symbol('z')
F = sp.Function('F')(z)
equation = 1 + z * F**2 - F
display(equation)

z F^{2}{\left (z \right )} - F{\left (z \right )} + 1

solution = sp.solve(equation, F)[0]
display(solution)

\frac{1}{2 z} \left(- \sqrt{- 4 z + 1} + 1\right)

display(sp.series(solution))

enter image description here

Question. The approach where we explicitly solve the equation and then expand it as power series, works only for low-degree equations. How to obtain first coefficients of formal power series for more complicated algebraic equations?

Related.

Since algebraic and differential framework may behave differently, I posted another question.

Sympy: how to solve differential equation in formal power series?

Sergey Dovgal
  • 614
  • 6
  • 21

1 Answers1

2

I don't know a built-in way, but plugging in a polynomial for F and equating the coefficients works well enough. Although one should not try to find all coefficients at once from a large nonlinear system; those will give SymPy trouble. I take iterative approach, first equating the free term to zero and solving for c0, then equating 2nd and solving for c1, etc.

This assumes a regular algebraic equation, in which the coefficient of z**k in the equation involves the k-th Taylor coefficient of F, and does not involve higher-order coefficients. 

from sympy import *
z = Symbol('z')
d = 10                                 # how many coefficients to find
c = list(symbols('c:{}'.format(d)))    # undetermined coefficients
for k in range(d):
    F = sum([c[n]*z**n for n in range(k+1)])  # up to z**k inclusive
    equation = 1 + z * F**2 - F
    coeff_eqn = Poly(series(equation, z, n=k+1).removeO(), z).coeff_monomial(z**k)
    c[k] = solve(coeff_eqn, c[k])[0]
sol = sum([c[n]*z**n for n in range(d)])  # solution
print(series(sol + z**d, z, n=d))         # add z**d to get SymPy to print as a series

This prints 

1 + z + 2*z**2 + 5*z**3 + 14*z**4 + 42*z**5 + 132*z**6 + 429*z**7 + 1430*z**8 + 4862*z**9 + O(z**10)
  • 1
    Great post! This approach is quite universal and also works on e.g. Ricatti equation `L = z^2 + z^2 L + z^4 L'` representing divergent formal power series. It has to be slightly adjusted for alternating permutations case `F' = F^2 + 1` to incorporate initial conditions, but this is not hard to modify. I really like your idea. – Sergey Dovgal Sep 26 '17 at 21:27
  • Maybe I should include differential case in the current question and remove another one? – Sergey Dovgal Sep 26 '17 at 21:45
  • I think they are fine as is; they are sufficiently different to be separate questions. And you already have links between them, so they won't be lost. Feel free to post a solution with your adjustments for `F'=F^2+1` to that other question. –  Sep 26 '17 at 21:50