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