Following this discussion about using asymptotic expansions (i.e. polynomial power series) for numerically solving partial differential and algebraic equations (PDAE), I couldn't find any implementation of the method. So I'm thinking of implementing a SymPy function similar to Mathematica's AsymptoticDSolveValue (here) but for PDEs (here or here). So far I'm able to generate a symbolic multivariate polynomial given a list of non-negative integers D=[d1,...,dm] (here, here and here).


Now I'm able to use the symbolic multivariate polynomial to numerically solve a PDE. For example (from here) given the PDE:

$$\frac{\partial^2 u}{\partial x_1^2} - \frac{\partial^2 u}{\partial x_2^2}=0 \, ,\tag{1}$$

and the boundary conditions:

  1. $u(x_1,0)=x_1^2+x_1\, , \tag{2}$
  2. $u_{x_2}(x_1,0)=2x_1+1 \, , \tag{3}$

I can generate a 2D symbolic polynomial in Sympy:

from itertools import product
from sympy import IndexedBase, symbols, Poly
D = (5, 5) # 5 and 5 are just some random integers, any non-negative integer should do
d = len(D)
indices=list(product(*map(range, D)))
a = IndexedBase('a')
coeffs = {i: a[i] for i in indices}
vars = symbols(f'x1:{d+1}')
u = Poly(coeffs, *vars)

Equation 1 can be implemented as:



$$j(j+1)a_{i,j+1}=(i+1)(i+2)a_{i+2,j} \tag{4}$$

The first boundary condition:



$$a_{0,0}=0 \, , \, a_{1,0}=1 \, , \, a_{2,0}=1 \, , \,a_{i,0}=0 \, \forall \,(2<i) \tag{5}$$

and the second boundary condition



$$a_{0,1}=1 \, , \, a_{1,1}=2 \, , \, a_{i,1}=0 \, \forall \, (2<i) \tag{6}$$

Now from this point it is just a system of nonlinear algebraic equation of $a$s (in this specific case linear). Which should be solvable with other analytical/numerical methods.


I want to automate the process above. I want to have a function:


Where eqns is the set of symbolic partial differential expressions, fs are the set of functions we want to solve, vars are the variables and D is the dimension of our multivariate polynomials.

I would appreciate if you could help me know what is the best algorithm for this process. I will use an existing symbolic library/software like SymPy, so anything already existing doesn't need to be reimplemented.

P.S. Axiom-FriCAS has a seriesSolve function which source code can be found here. Also since I posted this question, Nicolas CELLIER was so kind to implement an early version in sympy which can be seen in this Jupyter notebook.

  • 1
    $\begingroup$ @Evil I'm not asking for a software or a specific programing language. I' asking for the algorithm and used sympy just as an example. $\endgroup$
    – Foad
    Aug 2, 2018 at 12:37
  • 2
    $\begingroup$ Welcome to CS.SE! It seems like you already have described an algorithm (and it's a nice one!). It seems to generalize in a straightforward manner. I'm wondering what kind of answer you are looking for and what is the step you don't know how to complete on your own. Am I missing some difficulty? $\endgroup$
    – D.W.
    Aug 2, 2018 at 16:30
  • $\begingroup$ @D.W. thanks a lot for the warm and with no downvote welcome :) I indeed can solve specific problems manually as I demonstrated in the post, however my intention is to make all these automatic. Although it might seem easy, there is some logic behind which i can' generalize. What I expect is to have a step by step algorithm in the form of pseudocode preferably or flowcharts which I can implement later myself in SymPy or any of other libraries I have listed here. $\endgroup$
    – Foad
    Aug 2, 2018 at 18:19
  • 1
    $\begingroup$ OK. Which part are you having difficulty generalizing? On first glance it looks like each of the steps generalizes in a natural way, but I haven't thought about it as much as you undoubtedly have, so I'm probably missing something. $\endgroup$
    – D.W.
    Aug 2, 2018 at 18:44
  • $\begingroup$ @D.W. for start, I developed equation 4 intuitively by looking into the long list of coefficients generated from the command before. I don't know how a program can find these form of patterns, given the list. secondly equations 4 or 5 should be used to calculate the coefficient with higher indices. I don't know how to algorithmically do that. $\endgroup$
    – Foad
    Aug 2, 2018 at 18:55

1 Answer 1


Instead of building a polynomial and trying to deduce the equation on the coefficients, you can just derive it directly from a general term coefficient

>>> i, j, x1, x2 = symbols('i j x1 x2')
>>> a = IndexedBase('a')
>>> u = Function('u')
>>> pde = u(x1, x2).diff(x1, 2) - u(x1, x2).diff(x2, 2)
>>> term = a[i, j]*x1**i*x2**j
>>> pdeterm = powsimp(pde.replace(u, Lambda((x1, x2), term)).doit())
>>> pdeterm
i*x1**(i - 2)*x2**j*(i - 1)*a[i, j] - j*x1**i*x2**(j - 2)*(j - 1)*a[i, j]

The PDE being identically zero means each coefficient of each x1, x2 term is 0.

>>> pdeterm.coeff(x1**i*x2**(j - 2))
-j*(j - 1)*a[i, j]
>>> pdeterm.coeff(x1**(i - 2)*x2**j)
i*(i - 1)*a[i, j]

This is the same as what you had above if you shift i and j

>>> pdeterm.coeff(x1**i*x2**(j - 2)).subs({j: j + 1})
-j*(j + 1)*a[i, j + 1]
>>> pdeterm.coeff(x1**(i - 2)*x2**j).subs({i: i + 2})
(i + 1)*(i + 2)*a[i + 2, j]

I'm not sure if there's an advantage to having them one way or another. For a general PDE, you'd want to write a function that does this automatically.

The boundary conditions can be managed in a similar way using subs.

You can then solve the system analytically with solve or numerically with nsolve.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.