Algorithm for using power series to numerically solve a partial differential equation given a boundary condition?

Question

Motivation:

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).

Example:

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:

pde=u.diff(0,0).add(-u.diff(1,1))

Implying

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

The first boundary condition:

u.eval(1,0)

Implying

$$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

u.diff(1).eval(1,0)

giving

$$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.

Question:

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

AsymptoticPdeSolve(eqns,fs,vars,D)

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.

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.