$\begingroup$

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:

$u(x_1,0)=x_1^2+x_1\, , \tag{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.